diff --git a/Analyzer.C b/Analyzer.C index 322d601..60d023f 100644 --- a/Analyzer.C +++ b/Analyzer.C @@ -12,18 +12,20 @@ //############################################ User setting int rawEnergyRange[2] = {0, 6000}; // in ch -int energyRange[3] = {1, 0, 6000}; // keV {resol, min, max} +int energyRange[3] = {1, 50, 2000}; // keV {resol, min, max} double BGO_threshold = 0; // in ch -TString e_corr = "";//"correction_e.dat"; +TString e_corr = "correction_e.dat"; + +bool save_ev2 = true; //############################################ end of user setting //############################################ histogram declaration -TH2F * hevID; +TH2F * heVID; TH1F * he[NCRYSTAL]; TH2F * hgg[NCRYSTAL][NCRYSTAL]; @@ -31,7 +33,7 @@ TH2F * hgg[NCRYSTAL][NCRYSTAL]; TH2F * hcoin; ///----- after calibration and BGO veto -TH2F * heCalvID; +TH2F * heCalVID; TH1F * heCal[NCRYSTAL]; TH2F * hcoinBGO; @@ -51,11 +53,11 @@ void Analyzer::Begin(TTree * tree){ printf("======================== Histograms declaration\n"); - hevID = new TH2F("hevID", "e vs ID; det ID; e [ch]", NCRYSTAL, 0, NCRYSTAL, rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]); - heCalvID = new TH2F("heCalvID", Form("eCal vs ID (BGO veto > %.1f); det ID; Energy [keV]", BGO_threshold), NCRYSTAL, 0, NCRYSTAL, (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + heVID = new TH2F("heVID", "e vs ID; det ID; e [ch]", NCRYSTAL, 0, NCRYSTAL, rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]); + heCalVID = new TH2F("heCalVID", Form("eCal vs ID (BGO veto > %.1f); det ID; Energy [keV]", BGO_threshold), NCRYSTAL, 0, NCRYSTAL, (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); - hevID->SetNdivisions(-409, "X"); - heCalvID->SetNdivisions(-409, "X"); + heVID->SetNdivisions(-409, "X"); + heCalVID->SetNdivisions(-409, "X"); for( int i = 0; i < NCRYSTAL; i ++){ he[i] = new TH1F( Form("he%02d", i), Form("e -%02d", i), rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]); @@ -75,13 +77,14 @@ void Analyzer::Begin(TTree * tree){ hcoinBGO = new TH2F("hcoinBGO", Form("detector coin. (BGO veto > %.1f); det ID; det ID", BGO_threshold), NCRYSTAL, 0, NCRYSTAL, NCRYSTAL, 0 , NCRYSTAL); hcrystalBGO = new TH2F("hcrystalBGO", Form("crystal vs BGO ; det ID; BGO ID"), NCRYSTAL, 0, NCRYSTAL, NBGO, 0 , NBGO); - printf("======================== End of histograms declaration\n"); - printf("======================== Load parameters.\n"); eCorr = LoadCorrectionParameters(e_corr); + + saveEV2 = save_ev2; } + //############################################ PROCESS Bool_t Analyzer::Process(Long64_t entry){ @@ -99,8 +102,8 @@ Bool_t Analyzer::Process(Long64_t entry){ b_energy->GetEntry(entry); b_time->GetEntry(entry); - b_pileup->GetEntry(entry); - b_hit->GetEntry(entry); + //b_pileup->GetEntry(entry); + //b_hit->GetEntry(entry); b_bgo->GetEntry(entry); b_bgoTime->GetEntry(entry); b_other->GetEntry(entry); @@ -108,10 +111,7 @@ Bool_t Analyzer::Process(Long64_t entry){ if( multi == 0 ) return kTRUE; - int numGatedData=0; - vector gatedID; - gatedID.clear(); - double eCal[NCRYSTAL] ={0.0}; + for( int i = 0; i < NCRYSTAL; i++) eCal[i] = TMath::QuietNaN(); ///=========== Looping Crystals for( int detID = 0; detID < NCRYSTAL ; detID ++){ @@ -121,7 +121,7 @@ Bool_t Analyzer::Process(Long64_t entry){ //if( pileup[detID] == 1 ) continue; //======== Fill raw data - hevID->Fill( detID, e[detID]); + heVID->Fill( detID, e[detID]); he[detID]->Fill(e[detID]); for( int kk = 0 ; kk < NBGO; kk++){ @@ -142,20 +142,19 @@ Bool_t Analyzer::Process(Long64_t entry){ return kTRUE; } } - - //----- for ev2 file - if( saveEV2 ){ - numGatedData ++; - gatedID.push_back(detID); - } //========= apply correction - int order = (int) eCorr[detID].size(); - for( int i = 0; i < order ; i++){ - eCal[detID] += eCorr[detID][i] * TMath::Power(e[detID], i); + //int order = (int) eCorr[detID].size(); + //for( int i = 0; i < order ; i++){ + // eCal[detID] += eCorr[detID][i] * TMath::Power(e[detID], i); + //} + if( e_corr == "" ){ + eCal[detID] = e[detID]; + }else{ + eCal[detID] = eCorr[detID][0] + eCorr[detID][1] * e[detID]; } - - heCalvID->Fill( detID, eCal[detID]); + + heCalVID->Fill( detID, eCal[detID]); heCal[detID]->Fill(eCal[detID]); for( int detJ = detID +1; detJ < NCRYSTAL; detJ++) { @@ -165,38 +164,7 @@ Bool_t Analyzer::Process(Long64_t entry){ } - if ( saveEV2){ - - short * out0 = new short[1]; - short * outa = new short[1]; - short * outb = new short[1]; - - out0[0] = numGatedData; - fwrite(out0, 1, 1, outEV2); - - for( int i = 0; i < (int) gatedID.size(); i++){ - int id = gatedID[i]; - outa[0] = id; - fwrite(outa, 1, 1, outEV2); - outb[0] = eCal[id]; - fwrite(outb, 2, 1, outEV2); - } - - fwrite(out0, 1, 1, outEV2); - - /** - int len = (int) gatedID.size(); - char out[2*len+2]; - out[0] = numGatedData; - for( int i = 0; i < (int) gatedID.size(); i++){ - int id = gatedID[i]; - out[2*i+1] = id; - out[2*i+2] = eCal[id]; - } - out[2*len+1]=numGatedData; - fwrite(out,3*len+2,sizeof(out),outEV2); - */ - } + if ( saveEV2) Save2ev2(); return kTRUE; @@ -219,11 +187,11 @@ void Analyzer::Terminate(){ cCanvas->cd(1); cCanvas->cd(1)->SetLogz(1); - hevID->Draw("colz"); + heVID->Draw("colz"); cCanvas->cd(2); cCanvas->cd(2)->SetLogz(1); - heCalvID->Draw("colz"); + heCalVID->Draw("colz"); cCanvas->cd(3); cCanvas->cd(3)->SetLogz(1); diff --git a/Analyzer.h b/Analyzer.h index 1b09294..50a9b83 100644 --- a/Analyzer.h +++ b/Analyzer.h @@ -79,12 +79,11 @@ public : Float_t Frac; ///Progress bar TStopwatch StpWatch; + void Save2ev2(); bool saveEV2; FILE * outEV2; TString outEV2Name; - - - + double eCal[NCRYSTAL]; }; #endif @@ -108,15 +107,18 @@ void Analyzer::Init(TTree *tree) fChain->SetBranchAddress("evID", &evID, &b_event_ID); fChain->SetBranchAddress("e", e, &b_energy); fChain->SetBranchAddress("e_t", e_t, &b_time); - fChain->SetBranchAddress("p", p, &b_pileup); - fChain->SetBranchAddress("hit", hit, &b_hit); + //fChain->SetBranchAddress("p", p, &b_pileup); + //fChain->SetBranchAddress("hit", hit, &b_hit); fChain->SetBranchAddress("bgo", bgo, &b_bgo); fChain->SetBranchAddress("bgo_t", bgo_t, &b_bgoTime); fChain->SetBranchAddress("other", other, &b_other); fChain->SetBranchAddress("multi", &multi, &b_multi); + TString option = GetOption(); + if ( option != "" ) outEV2Name = option; + if( saveEV2 ){ - printf("======================== Create output ev2 : %s\n", outEV2Name.Data()); + printf("======================== Create output ev2 : \033[1;31m%s\033[0m\n", outEV2Name.Data()); outEV2 = fopen(outEV2Name.Data(), "w+"); } @@ -150,4 +152,42 @@ void Analyzer::SlaveTerminate(){ } +void Analyzer::Save2ev2(){ + short * out0 = new short[1]; + short * outa = new short[1]; + short * outb = new short[1]; + + int count = 0; + for (int i = 0; i < NCRYSTAL ; i++){ + if( !TMath::IsNaN(eCal[i]) ) count++; + } + + out0[0] = count; + fwrite(out0, 1, 1, outEV2); + + for( int i = 0; i < count; i++){ + if( TMath::IsNaN(eCal[i]) ) continue; + outa[0] = i; + fwrite(outa, 1, 1, outEV2); + outb[0] = eCal[i]; + fwrite(outb, 2, 1, outEV2); + } + + fwrite(out0, 1, 1, outEV2); + + + /** + int len = (int) gatedID.size(); + char out[2*len+2]; + out[0] = numGatedData; + for( int i = 0; i < (int) gatedID.size(); i++){ + int id = gatedID[i]; + out[2*i+1] = id; + out[2*i+2] = eCal[id]; + } + out[2*len+1]=numGatedData; + fwrite(out,3*len+2,sizeof(out),outEV2); + */ +} + #endif // #ifdef Analyzer_cxx diff --git a/armory/AnalysisLibrary.h b/armory/AnalysisLibrary.h index ffb3b39..8bfe3c1 100644 --- a/armory/AnalysisLibrary.h +++ b/armory/AnalysisLibrary.h @@ -233,7 +233,7 @@ std::vector> FindMatchingPair(std::vector enX, std:: } -std::vector> LoadCorrectionParameters(TString corrFile){ +std::vector> LoadCorrectionParameters(TString corrFile, bool show=false){ printf(" load correction parameters : %s", corrFile.Data()); std::ifstream file; @@ -268,14 +268,16 @@ std::vector> LoadCorrectionParameters(TString corrFile){ file.close(); printf(".... done\n"); - printf("===== correction parameters \n"); - for( int i = 0; i < (int) corr.size(); i++){ - printf("det : %2d | ", i ); - int len = (int) corr[i].size(); - for( int j = 0; j < len - 1 ; j++){ - printf("%14.6f, ", corr[i][j]); + if( show ){ + printf("===== correction parameters \n"); + for( int i = 0; i < (int) corr.size(); i++){ + printf("det : %2d | ", i ); + int len = (int) corr[i].size(); + for( int j = 0; j < len - 1 ; j++){ + printf("%14.6f, ", corr[i][j]); + } + printf("%14.6f\n", corr[i][len-1]); } - printf("%14.6f\n", corr[i][len-1]); } }else{ diff --git a/armory/Analyzer_Utili.c b/armory/Analyzer_Utili.c index eacc73f..99a9cdd 100644 --- a/armory/Analyzer_Utili.c +++ b/armory/Analyzer_Utili.c @@ -2,7 +2,7 @@ void listDraws(void) { printf("------------------- List of Plots -------------------\n"); printf(" newCanvas() - Create a new Canvas\n"); printf("-----------------------------------------------------\n"); - printf(" rawEvID() - e vs ID\n"); + printf(" eVID() - e vs ID\n"); printf(" drawE() - e for all %d detectors\n", NCRYSTAL); //printf(" drawGG() - Gamma - Gamma Coincident for all %d detectors\n", NCRYSTAL); printf("-----------------------------------------------------\n"); @@ -18,12 +18,16 @@ void newCanvas(int sizeX = 800, int sizeY = 600, int posX = 0, int posY = 0){ cNewCanvas->cd(); } -//void rawEvID(bool cal = false){ -// TCanvas * cRawID = (TCanvas *) gROOT->FindObjectAny("cRawID"); -// if( cRawID == NULL ) cRawID = new TCanvas("cRawID", "raw ID", 1000, 800); -// cRawID->cd(1)->SetGrid(); -// cal ? heCalVID->Draw("colz") : heVID->Draw("colz"); -//} +void eVID(bool cal = false, bool logz = false){ + TCanvas * cRawID = (TCanvas *) gROOT->FindObjectAny("cRawID"); + if( cRawID == NULL ) cRawID = new TCanvas("cRawID", "raw ID", 1000, 800); + + if( cRawID->GetShowEventStatus() == 0 ) cRawID->ToggleEventStatus(); + + cRawID->cd(1)->SetGrid(); + if( logz ) cRawID->cd(1)->SetLogz(); + cal ? heCalVID->Draw("colz") : heVID->Draw("colz"); +} void drawE(int CloverID = -1, bool cali = false, bool isLogy = false, double xMin = 0, double xMax = 0){ @@ -39,10 +43,16 @@ void drawE(int CloverID = -1, bool cali = false, bool isLogy = false, double xMi TCanvas *cRawE = (TCanvas *) gROOT->FindObjectAny("cRawE"); if( cRawE == NULL ) cRawE = new TCanvas("cRawE", cali ? "Cal e" : "Raw e", size * nClover, size * nCrystalPerClover); + + if( cRawE->GetShowEventStatus() == 0 ) cRawE->ToggleEventStatus(); + cRawE->Clear(); - if( CloverID >= 0 ) nClover = 1; - cRawE->Divide(nClover, nCrystalPerClover, 0); - + if( CloverID >= 0 ) { + nClover = 1; + cRawE->Divide(nClover, 1); + }else{ + cRawE->Divide(nClover, nCrystalPerClover, 0); + } ///find max y double maxY = 0; @@ -56,6 +66,18 @@ void drawE(int CloverID = -1, bool cali = false, bool isLogy = false, double xMi ///printf("max Y : %f \n", maxY); for (Int_t i = 0; i < nClover; i++) { + + int hID = nCrystalPerClover * CloverID + i ; + if( cali ) { + heCal[hID]->SetMaximum(maxY); + heCal[hID]->Draw(""); + }else{ + he[hID]->SetMaximum(maxY); + he[hID]->Draw(""); + } + + + /* for( Int_t j = 0; j < nCrystalPerClover; j++){ int canvasID = CloverID < 0 ? nClover*j+ i + 1 : j + 1; cRawE->cd(canvasID); @@ -75,7 +97,7 @@ void drawE(int CloverID = -1, bool cali = false, bool isLogy = false, double xMi he[hID]->SetMaximum(maxY); he[hID]->Draw(""); } - } + }*/ } cRawE->SetCrosshair(1); diff --git a/process_run b/process_run index 9405d4a..94fff15 100755 --- a/process_run +++ b/process_run @@ -19,6 +19,8 @@ isMerge=1 if [ $# -gt 1 ]; then isMerge=$2; fi isBuildEvents=1 if [ $# -gt 2 ]; then isBuildEvents=$3; fi +isAnalysis=1 +if [ $# -gt 3 ]; then isAnalysis=$4; fi RED='\033[1;31m' YELLOW='\033[1;33m' @@ -28,6 +30,12 @@ BLUE='\033[0;34m' Cyan='\033[0;36m' NC='\033[0m' +if [ -f $DATA_DIR/$RunFolder/*.evt ]; then + echo -e "found evt files." +else + echo -e "cannot found any evt files. Abort." + exit +fi echo -e "$RED>>> `date` >>>>>>>>>>>>>>>>>>>>>>> Merge evt files to ${RunFolder}_raw.root $NC" @@ -86,94 +94,17 @@ fi if [ ${isBuildEvents} -eq -1 ]; then echo -e "$YELLOW forced by user $NC" ./armory/EventBuilder ${RunFolder}"_raw.root" - fi echo -e "$RED>>> `date` >>>>>>>>>>>>>>>>>>>>>>> Build Events finsihed.$NC" echo -e "$RED>>> `date` >>>>>>>>>>>>>>>>>>>>>>> Analysis $NC" +if [ ${isAnalysis} -eq 1 ]; then + root -l "process_run.c(\"${RunFolder}.root\")" +fi +if [ ${isAnalysis} -eq -1 ]; then + echo -e "$YELLOW forced by user $NC" + root -l "process_run.c(\"${RunFolder}.root\")" +fi - - - - -#cd $RunFolder -#ls -lhtr *.evt -#cd $DATA_DIR/$RunFolder -#fileList=$(ls *.evt) -#numFile=$(ls -lhtr *.evt | wc -l) -#count=0 -# -#cd $DIR - - - - - - -#for a in $fileList -#do -# count=$((${count}+1)) -# echo -e "$YELLOW>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ["$count"/"$numFile"]" -# echo -e $a -# echo -e "------------------------------------------$NC" -# -# #............ check if *.to file exist, if yes, check timestamp, if timestamp is eariler, sort time -# -# if [ $isMerge -eq 1 ]; then -# evtDateTime=`stat -c "%Z" $DATA_DIR/$RunFolder/$a | sort -rn | head -1` -# istoExist=`ls -1 $a.to 2> /dev/null | wc -l` -# if [ ${istoExist} -gt 0 ]; then -# echo "--- found $a.to" -# toDateTime=`stat -c "%Z" $a.to | sort -rn | head -1` -# if [ ${evtDateTime} -ge ${toDateTime} ]; then -# $sortTime $DATA_DIR/$RunFolder/$a -# else -# echo "$a.to is generated after $a. skip." -# fi -# else -# echo "cannot find $a.to, sort" -# $sortTime $DATA_DIR/$RunFolder/$a -# fi -# elif [ $isMerge -eq 0 ]; then -# echo "skipped time sort by user." -# else -# echo "force Time sort" -# $sortTime $DATA_DIR/$RunFolder/$a -# fi -# -# -# #............ check if *.root file exist, if yes, check timestamp, if timestamp is eariler, built -# -# if [ $isBuildEvents -eq 1 ]; then -# len=`echo $a | wc -c` -# len=$(($len-5)) -# rootFile=${a:0:$len}".root" -# isRootExist=`ls -1 $rootFile 2> /dev/null | wc -l` -# if [ ${isRootExist} -gt 0 ]; then -# echo "--- found $rootFile" -# rootDateTime=`stat -c "%Z" $rootFile | sort -rn | head -1` -# if [ ${toDateTime} -ge ${rootDateTime} ]; then -# $DIR/pixie2root $a.to -# else -# echo "$rootFile is generated after $a.to skip." -# fi -# else -# echo "cannot find $rootFile, build events" -# $DIR/pixie2root $a.to -# fi -# elif [ $isMerge -eq 0 ]; then -# echo "skipped event build by user." -# else -# echo "force event build" -# $DIR/pixie2root $a.to -# fi -# -#done -# -#fileList=$(ls *.root) - - - -cd $DIR diff --git a/process_run.c b/process_run.c new file mode 100644 index 0000000..84c30b5 --- /dev/null +++ b/process_run.c @@ -0,0 +1,15 @@ + + +void process_run(TString rootFile){ + + TFile * file = new TFile(rootFile, "read"); + + TTree * tree = (TTree *) file->Get("tree"); + + TString ev2FileName = rootFile; + ev2FileName.Remove(rootFile.First(".")); + ev2FileName.Append(".ev2"); + + tree->Process("Analyzer.C+", ev2FileName); + +}