diff --git a/.gitignore b/.gitignore index 9719c91..9dcb86b 100644 --- a/.gitignore +++ b/.gitignore @@ -12,6 +12,7 @@ pxi-time-order evt2root evt2hist ev22txt +EventBuilder *.so *.d diff --git a/Monitor_raw.C b/Monitor_raw.C index cede9ba..34ebf03 100644 --- a/Monitor_raw.C +++ b/Monitor_raw.C @@ -1,31 +1,31 @@ #define Monitor_raw_cxx -#include "Monitor_raw.h" + #include #include #include + #include #include #include -#include + //############################################ User setting int rawEnergyRange[2] = {500, 6000}; // in ch, {min, max} -int energyRange[3] = {2, 40, 2000}; // in keV, {resol, min, max} +int energyRange[3] = {1, 40, 2000}; // in keV, {resol, min, max} TString e_corr = "correction_e.dat"; Long64_t maxEvent = -1; +bool isBuildEventRoot = true; // is true, no histogram + + //############################################ end of user setting -ULong64_t NumEntries = 0; -ULong64_t ProcessedEntries = 0; -Float_t Frac = 0.1; ///Progress bar -TStopwatch StpWatch; +#include "Monitor_raw.h" -vector> eCorr; //############################################ histogram declaration @@ -37,39 +37,50 @@ TH2F * heCalvID; //############################################ BEGIN void Monitor_raw::Begin(TTree * tree){ TString option = GetOption(); - - NumEntries = tree->GetEntries(); - - printf("======================== Creating histograms\n"); - 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]); - heCal[i] = new TH1F (Form("heCal%02d", i), Form("e%02d (Cali.)", i), (energyRange[2]-energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + totnumEntry = tree->GetEntries(); + + printf( "=========================================================================== \n"); + printf( "========================== Monitor_raw.C/h ================================ \n"); + printf( "====== total Entry : %lld \n", totnumEntry); + printf( "=========================================================================== \n"); + + + if( isBuildEvent == false && isBuildEventRoot == false) { + printf("======================== Creating histograms\n"); + + 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]); + heCal[i] = new TH1F (Form("heCal%02d", i), Form("e%02d (Cali.)", i), (energyRange[2]-energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + + switch (i%4){ + case 0: he[i]->SetLineColor(2);break; + case 1: he[i]->SetLineColor(kYellow+3);break; + case 2: he[i]->SetLineColor(kGreen+2);break; + case 3: he[i]->SetLineColor(4);break; + } + switch (i%4){ + case 0: heCal[i]->SetLineColor(2);break; + case 1: heCal[i]->SetLineColor(kYellow+3);break; + case 2: heCal[i]->SetLineColor(kGreen+2);break; + case 3: heCal[i]->SetLineColor(4);break; + } + } - switch (i%4){ - case 0: he[i]->SetLineColor(2);break; - case 1: he[i]->SetLineColor(kYellow+3);break; - case 2: he[i]->SetLineColor(kGreen+2);break; - case 3: he[i]->SetLineColor(4);break; - } - switch (i%4){ - case 0: heCal[i]->SetLineColor(2);break; - case 1: heCal[i]->SetLineColor(kYellow+3);break; - case 2: heCal[i]->SetLineColor(kGreen+2);break; - case 3: heCal[i]->SetLineColor(4);break; - } - } - - heCalvID = new TH2F("heCalvID", "ID vs Energy (Cali.); ID; Energy", NCRYSTAL, 0, NCRYSTAL, (energyRange[2]-energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + heCalvID = new TH2F("heCalvID", "ID vs Energy (Cali.); ID; Energy", NCRYSTAL, 0, NCRYSTAL, (energyRange[2]-energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); - - printf("======================== end of histograms creation.\n"); + + printf("------------------------ end of histograms creation.\n"); + } printf("======================== Load parameters.\n"); eCorr = LoadCorrectionParameters(e_corr); - StpWatch.Start(); - printf("======================== Start processing....\n"); + time0 = 0; + timeDiff = 0; + + ClearTreeData(); + } @@ -79,11 +90,11 @@ Bool_t Monitor_raw::Process(Long64_t entry){ ProcessedEntries++; /*********** Progress Bar ******************************************/ - if (ProcessedEntries>NumEntries*Frac-1) { - TString msg; msg.Form("%llu", NumEntries/1000); + if (ProcessedEntries>totnumEntry*Frac-1) { + TString msg; msg.Form("%llu", totnumEntry/1000); int len = msg.Sizeof(); printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n", - Frac*100, len, ProcessedEntries/1000,NumEntries/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac); + Frac*100, len, ProcessedEntries/1000,totnumEntry/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac); StpWatch.Start(kFALSE); Frac+=0.1; } @@ -93,22 +104,25 @@ Bool_t Monitor_raw::Process(Long64_t entry){ Terminate(); } + if( isBuildEventRoot || isBuildEvent ) entry = index[entry]; + b_ID->GetEntry(entry); b_energy->GetEntry(entry); - b_timestamp->GetEntry(entry); + b_energy_timestamp->GetEntry(entry); - int detID = mapping[ID]; - - if( 0 <= detID && detID < 100 ){ - he[detID]->Fill(e); - - //double eCal = ApplyCorrection(eCorr, detID, e); - double eCal = eCorr[detID][0]+eCorr[detID][1]*e; - heCal[detID]->Fill(eCal); - - heCalvID->Fill(detID, eCal); + if ( isBuildEventRoot || isBuildEvent ) { + BuildEvent(); + }else{ + if( detID < 100 ){ + he[detID]->Fill(energy); + + ///double eCal = ApplyCorrection(eCorr, detID, e); //slow, why? + double eCal = eCorr[detID][0] + eCorr[detID][1]*energy; + heCal[detID]->Fill(eCal); + + heCalvID->Fill(detID, eCal); + } } - return kTRUE; } @@ -116,38 +130,47 @@ Bool_t Monitor_raw::Process(Long64_t entry){ //############################################ TERMINATE void Monitor_raw::Terminate(){ - printf("============================== finishing.\n"); + + printf("============================== finished.\n"); gROOT->cd(); - int nCrystalPerClover = 4; - int nClover = NCRYSTAL / nCrystalPerClover; - - TCanvas * cc = new TCanvas("cc", "cc", 2000, 2000); - if( cc->GetShowEventStatus() == 0 ) cc->ToggleEventStatus(); - cc->Divide(1, 9, 0); + if ( isBuildEventRoot || isBuildEvent ) { + saveFile->cd(); + newtree->Write(); + saveFile->Close(); + }else{ - for (Int_t i = 0; i < nClover; i++) { - int canvasID = i + 1; - cc->cd(canvasID); - cc->cd(canvasID)->SetGrid(); - cc->cd(canvasID)->SetTickx(2); - cc->cd(canvasID)->SetTicky(2); - cc->cd(canvasID)->SetBottomMargin(0.06); - cc->cd(canvasID)->SetLogy(); + int nCrystalPerClover = 4; + int nClover = NCRYSTAL / nCrystalPerClover; - for( Int_t j = 0; j < nCrystalPerClover; j++){ - int hID = nCrystalPerClover*i+ j ; - heCal[hID]->Draw("same"); - //he[hID]->Draw("same"); + TCanvas * cc = new TCanvas("cc", "cc", 2000, 2000); + if( cc->GetShowEventStatus() == 0 ) cc->ToggleEventStatus(); + cc->Divide(1, 9, 0); + + for (Int_t i = 0; i < nClover; i++) { + int canvasID = i + 1; + cc->cd(canvasID); + cc->cd(canvasID)->SetGrid(); + cc->cd(canvasID)->SetTickx(2); + cc->cd(canvasID)->SetTicky(2); + cc->cd(canvasID)->SetBottomMargin(0.06); + cc->cd(canvasID)->SetLogy(); + + for( Int_t j = 0; j < nCrystalPerClover; j++){ + int hID = nCrystalPerClover*i+ j ; + heCal[hID]->Draw("same"); + //he[hID]->Draw("same"); + } } + cc->SetCrosshair(1); + + TCanvas * c1 = new TCanvas("c1", "c1", 1000, 1000); + if( c1->GetShowEventStatus() == 0 ) c1->ToggleEventStatus(); + c1->SetLogz(); + c1->SetGridx(); + heCalvID->SetNdivisions(-409, "X"); + heCalvID->Draw("colz"); + } - cc->SetCrosshair(1); - - TCanvas * c1 = new TCanvas("c1", "c1", 1000, 1000); - if( c1->GetShowEventStatus() == 0 ) c1->ToggleEventStatus(); - c1->SetLogz(); - c1->SetGridx(); - heCalvID->SetNdivisions(-409, "X"); - heCalvID->Draw("colz"); } diff --git a/Monitor_raw.h b/Monitor_raw.h index 9bac8a4..46cc0c6 100644 --- a/Monitor_raw.h +++ b/Monitor_raw.h @@ -3,8 +3,12 @@ #include #include +#include #include +#include #include +#include +#include #include "mapping.h" #include "armory/AnalysisLibrary.h" @@ -19,17 +23,20 @@ public : // Declaration of leaf types Long64_t evID; - UShort_t ID; - UShort_t e; - ULong64_t t; + UShort_t detID; + UShort_t energy; + ULong64_t energy_t; // List of branches TBranch *b_data_ID; //! TBranch *b_ID; //! TBranch *b_energy; //! - TBranch *b_timestamp; //! + TBranch *b_energy_timestamp; //! - Monitor_raw(TTree * /*tree*/ =0) : fChain(0) { } + Monitor_raw(TTree * /*tree*/ =0) : fChain(0) { isBuildEvent = false; + timeWindow = 100; + saveFileName = "test.root"; + } virtual ~Monitor_raw() { } virtual Int_t Version() const { return 2; } virtual void Begin(TTree *tree); @@ -44,8 +51,44 @@ public : virtual TList *GetOutputList() const { return fOutput; } virtual void SlaveTerminate(); virtual void Terminate(); + + void SetBuildEvent(bool buildOption = false){ isBuildEvent = buildOption;} ClassDef(Monitor_raw,0); + + //=================== progress + ULong64_t ProcessedEntries = 0; + Float_t Frac = 0.1; ///Progress bar + TStopwatch StpWatch; + + //=================== correction parameters + vector> eCorr; + + //=================== output tree; + bool isBuildEvent; + int timeWindow; + ULong64_t time0; //time-0 for each event + int timeDiff; + Long64_t * index; //! + + TFile * saveFile; //! + TTree * newtree; //! + TString saveFileName; + Long64_t totnumEntry; /// of original root + + //tree + void ClearTreeData(); + void BuildEvent(); + + Int_t eventID; + double e[NCRYSTAL]; + ULong64_t e_t[NCRYSTAL]; + double bgo[NBGO]; + ULong64_t bgo_t[NBGO]; + Short_t other[NOTHER]; + Short_t multi; + + }; #endif @@ -55,15 +98,187 @@ void Monitor_raw::Init(TTree *tree) { // Set branch addresses and branch pointers if (!tree) return; - fChain = tree; + + fChain = (TChain *) tree; fChain->SetMakeClass(1); - fChain->SetBranchAddress("evID", &evID, &b_data_ID); - fChain->SetBranchAddress("id", &ID, &b_ID); - fChain->SetBranchAddress("e", &e, &b_energy); - fChain->SetBranchAddress("t", &t, &b_timestamp); + fChain->SetBranchAddress("evID", &evID, &b_data_ID); + fChain->SetBranchAddress("id", &detID, &b_ID); + fChain->SetBranchAddress("e", &energy, &b_energy); + fChain->SetBranchAddress("e_t", &energy_t, &b_energy_timestamp); + + if( isBuildEventRoot || isBuildEvent){ + + printf("======================== Buidling Index using the timestamp\n"); + tree->BuildIndex("e_t"); + TTreeIndex *in = (TTreeIndex*) tree->GetTreeIndex(); + index = in->GetIndex(); + + //for(int i = 0; i < 100; i++){ + // printf(" %3d | %lld \n", i, index[i]); + //} + + + TString option = GetOption(); + + //printf("======================== Formation of output file name\n"); + //int numFile = fChain->GetListOfFiles()->GetLast() + 1; //need input of TChain + // + //printf(".......... number of files : %d \n", numFile); + // + //if( numFile > 0 ) { + // int oldRunNum = -100; + // bool contFlag = false; // is runNumber continue; + // for( int i = 0; i < numFile ; i++){ + // TString name = ((TChain *) fChain)->GetListOfFiles()->At(i)->GetTitle(); + // int found = name.Last('/'); + // name.Remove(0, found + 1 ); // this should give "XXX_run0XX.root" + // TString prefix = name; + // found = name.Last('.'); + // name.Remove(found); // this should give "XXX_run0XX" + // found = name.Last('_'); + // int runNum = name.Remove(0, found+4).Atoi(); // this should give the 3 digit run number + // + // if( i == 0 ) { + // found = prefix.First("_"); + // prefix.Remove(found); + // saveFileName = prefix + "_run"; + // } + // + // if( runNum == oldRunNum + 1 ){ + // int kk = saveFileName.Sizeof(); + // if( contFlag == false ){ + // saveFileName.Remove(kk-2); //remove the "-" + // saveFileName += "-"; + // }else{ + // saveFileName.Remove(kk-5); //remove the runNum and "-" + // } + // contFlag = true; + // } + // if( runNum > oldRunNum + 1) contFlag = false; + // + // saveFileName += Form("%03d_", runNum); + // oldRunNum = runNum; + // } + // int kk = saveFileName.Sizeof(); + // saveFileName.Remove(kk-2); // remove the last "-" + // saveFileName += ".root"; + //}else{ + // gROOT->ProcessLine(".q"); + //} + + if( option != "" ) saveFileName = option; + + printf("save file name : %s \n", saveFileName.Data()); + printf("---------------------------------------------\n"); + + //TODO; + saveFile = new TFile(saveFileName, "recreate"); + saveFile->cd(); + newtree = new TTree("tree", "tree"); + + printf("======================== Create output tree\n"); + + newtree->Branch("evID", &eventID, "event_ID/l"); + newtree->Branch("e", e, Form("e[%d]/D", NCRYSTAL)); + newtree->Branch("e_t", e_t, Form("e_timestamp[%d]/l", NCRYSTAL)); + //newtree->Branch("p", pileup, Form("pile_up_flag[%d]/s", NCRYSTAL)); + //newtree->Branch("hit", hit, Form("hit[%d]/s", NCRYSTAL)); + + newtree->Branch("bgo", bgo, Form("BGO_e[%d]/D", NBGO)); + newtree->Branch("bgo_t", bgo_t, Form("BGO_timestamp[%d]/l", NBGO)); + + newtree->Branch("other", other, Form("other_e[%d]/D", NOTHER)); + + newtree->Branch("multi", &multi, "multiplicity_crystal/I"); + + + if( eCorr.size() > 0 ) { + TMacro energyCorr(e_corr); + energyCorr.Write("energyCorr"); + } + + } + + + printf("======================== Start processing....\n"); + StpWatch.Start(); + eventID = 0; + } + +void Monitor_raw::ClearTreeData(){ + + for( int i = 0; i < NCRYSTAL; i++){ + e[i] = TMath::QuietNaN(); + e_t[i] = 0; + //pileup[i] = 0; + //hit[i] = 0; + } + for( int i = 0; i < NBGO; i++) { + bgo[i] = TMath::QuietNaN(); + bgo_t[i] = 0 ; + } + for( int i = 0; i < NOTHER; i++) { + other[i] = TMath::QuietNaN(); + } + multi = 0; +} + +void Monitor_raw::BuildEvent(){ + + if( time0 == 0) time0 = energy_t; + + timeDiff = (int) (energy_t - time0); + + if( timeDiff < timeWindow ) { + + if ( detID < NCRYSTAL ){ + e[detID] = energy; + e_t[detID] = energy_t; + multi++; + } + if ( 100 <= detID && detID < 100 + NBGO ){ + bgo[detID-100] = energy; + bgo_t[detID-100] = energy_t; + } + if ( 200 <= detID && detID < 200 + NOTHER){ + other[detID-200] = energy; + } + + //printf("%d | %3d %6d %10llu, %3d\n", multi, detID, energy, energy_t, timeDiff); + + }else{ + //---- end of event + eventID ++; + + saveFile->cd(); + newtree->Fill(); + + ClearTreeData(); + + /// fill 1st data of an event + time0 = energy_t; + timeDiff = 0; + + if ( detID < NCRYSTAL ){ + e[detID] = energy; + e_t[detID] = energy_t; + multi = 1; + } + if ( 100 <= detID && detID < 100 + NBGO ){ + bgo[detID-100] = energy; + bgo_t[detID-100] = energy_t; + } + if ( 200 <= detID && detID < 200 + NOTHER){ + other[detID-200] = energy; + } + + } +} + + Bool_t Monitor_raw::Notify() { return kTRUE;