#define DecayFinder_cxx #include "DecayFinder.h" #include #include #include #include #include #include #include #include //===================== Settings Bool_t debug = false; TString pidCutFileName = "PIDCuts.root"; double distThreshold = 0.2; double forwardTime = 250; //ms double backwardTime = 100; //ms //===================== histograms TH2F * hPID; TH1F ** hDecay ; TH2F * hDist; TH1F ** hvetoF; TH2F * hIonsPos; TH2F * hBetaPos; //===================== Parameters int tick2ns = 8; double tick2ms = tick2ns / 1e6; TFile * cutFile; TCutG * cut = NULL; TObjArray * cutList; int numCut = 0; int numMatch = 0; //===================== Begins void DecayFinder::Begin(TTree * /*tree*/){ TString option = GetOption(); hPID = new TH2F("hPID", "PID; TOF; dE", 400, -260, -220, 400, 0, 6500); hDist = new TH2F("hDist", "dist; dx; dy", 100, -distThreshold, distThreshold, 100, -distThreshold, distThreshold); hIonsPos = new TH2F("hIonsPos", "Ions Pos; x; y", 200, 0, 1, 200, 0, 1); hBetaPos = new TH2F("hBetaPos", "Beta Pos; x; y", 200, 0, 1, 200, 0, 1); //-------------- GetCut; cutFile = new TFile(pidCutFileName, "UPDATE"); bool listExist = cutFile->GetListOfKeys()->Contains("cutList"); if( !listExist ) { cutList = new TObjArray(); }else{ cutList = (TObjArray*) cutFile->FindObjectAny("cutList"); numCut = cutList->GetLast()+1; printf("------------ found %d cuts in %s \n", numCut, pidCutFileName.Data()); for( int k = 0; k < numCut; k++){ if( cutList->At(k) != NULL ){ printf("found a cut at %2d \n", k); }else{ printf(" No cut at %2d \n", k); } } } hDecay = new TH1F *[numCut]; hvetoF = new TH1F *[numCut]; for( int i = 0; i < numCut; i++){ hDecay[i] = new TH1F(Form("hDecay%02d",i), Form("Decay cut-%02d ; [ms]; count", i), 50, -backwardTime, forwardTime); hvetoF[i] = new TH1F(Form("hvetoF%02d",i), Form("veto-F cut-%02d ; [ms]; count", i), 100, 0, 6000); } printf("=====================================\n"); } double dist(double x1, double x2, double y1, double y2){ return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); } //===================== Process Bool_t DecayFinder::Process(Long64_t entry){ if( entry > (Long64_t)1e6 ) return kTRUE; b_TOF->GetEntry(entry); b_energy->GetEntry(entry); b_crossTime->GetEntry(entry); b_crossEnergy->GetEntry(entry); b_flag->GetEntry(entry); b_vetoFlag->GetEntry(entry); b_xIons->GetEntry(entry); b_yIons->GetEntry(entry); b_dyIonsTime->GetEntry(entry); b_veto_r->GetEntry(entry); b_veto_f->GetEntry(entry); ULong64_t ImplantEntry = 0; ULong64_t ImplantTime = 0; double ImplantX = TMath::QuietNaN(); double ImplantY = TMath::QuietNaN(); if( debug && entry % 10 == 0 ) printf("------------- %llu\n", entry); if( flag & 1 ) hPID->Fill(TOF, energy); int cutID = -1; for(int i = 0; i < numCut; i++){ cut = (TCutG*) cutList->At(i); if( cut->IsInside(TOF, energy) ) cutID = i; } if( cutID == -1 ) return kTRUE; if( (flag & 1) == 0 ) return kTRUE; /// no beam //if( (flag & 2) == 0 ) return kTRUE; /// no Ions if( (flag & 4) == 0 ) return kTRUE; /// no beta //if( veto_f > 0 ) return kTRUE; if( veto_r > 0 ) return kTRUE; if( dyIonsTime[0] == 0 ) return kTRUE; ImplantEntry = entry; ImplantTime = dyIonsTime[0]; ImplantX = xIons; ImplantY = yIons; hIonsPos->Fill( xIons, yIons); if( debug ) printf("===========%8lld, %13llu, %f, %f\n", ImplantEntry, ImplantTime, ImplantX, ImplantY); for( ULong64_t k = ImplantEntry - 10000; k < totNumEntry ; k++){ b_flag->GetEntry(k); b_vetoFlag->GetEntry(k); b_dyBetaTime->GetEntry(k); b_crossTime->GetEntry(k); b_crossEnergy->GetEntry(k); b_veto_f->GetEntry(k); b_veto_r->GetEntry(k); b_xBeta->GetEntry(k); b_yBeta->GetEntry(k); if ( k == ImplantEntry ) continue; //if( crossTime > 0 ) continue; if( crossEnergy > 0 ) continue; //if( (flag & 2) == 2 ) continue; /// has Ions if( (flag & 4) == 0 ) continue; /// no Beta if( !TMath::IsNaN(veto_f) ) continue; if( !TMath::IsNaN(veto_r) ) continue; if( dyBetaTime[0] == 0) continue; if( dist(ImplantX, xBeta, ImplantY, yBeta) > distThreshold ) continue; if( k < ImplantEntry && ( ImplantTime - dyBetaTime[0] ) * tick2ms > backwardTime) continue; if( k > ImplantEntry && ( dyBetaTime[0] - ImplantTime ) * tick2ms > forwardTime ) break; clock.Stop("timer"); Double_t time = clock.GetRealTime("timer"); clock.Start("timer"); printf(" %llu[%.2f%%], searched next %lld events | match case %d | expect : %5.2f min\r", ImplantEntry, ImplantEntry*100./totNumEntry, k - ImplantEntry, numMatch, totNumEntry*time/(ImplantEntry+1.)/60.); if( debug ) printf(" %8lld, %13llu, %f, %f\n", k, dyBetaTime[0], xBeta, yBeta); if( debug ) printf("+++++++++++%8s, %13llu | %.f ms\n", "+", dyBetaTime[0], (dyBetaTime[0] - ImplantTime)*tick2ms); hDist->Fill(ImplantX - xBeta, ImplantY - yBeta); numMatch ++; double dT = 0; if( dyBetaTime[0] > ImplantTime ){ dT = (dyBetaTime[0] - ImplantTime)*tick2ms; }else{ dT = (ImplantTime - dyBetaTime[0])*tick2ms; dT = -dT; } hDecay[cutID]->Fill(dT); hvetoF[cutID]->Fill(crossEnergy); hBetaPos->Fill(xBeta, yBeta); } return kTRUE; } void DecayFinder::Terminate(){ //============== Save the results TFile * haha = new TFile("results.root", "recreate"); haha->cd(); for( int i = 0; i < numCut; i++){ hDecay[i]->Write(); } hPID->Write(); haha->Close(); //============== Plot result gStyle->SetOptStat("neiou"); TCanvas * cDecay = new TCanvas("cDecay", "Decay", 1000, 1000); cDecay->Divide(2,2); int padID = 1; cDecay->cd(padID); hPID->Draw("colz"); for( int i = 0; i < numCut; i++){ cut = (TCutG*) cutList->At(i); cut->SetLineColor(i+2); cut->Draw("same"); } padID++; cDecay->cd(padID); //cDecay->cd(padID)->SetLogy(); double yMax = 0; for( int i = 0; i < numCut; i++){ if( hDecay[i]->GetMaximum() > yMax ) yMax = hDecay[i]->GetMaximum(); } for( int i = 0; i < numCut; i++){ hDecay[i]->SetLineColor(i+2); hDecay[i]->SetMaximum(yMax *1.2); hDecay[i]->Draw(i == 0 ? "" : "same"); } //padID++; cDecay->cd(padID); //hDist->Draw("colz"); // //padID++; cDecay->cd(padID); //for( int i = 0; i < numCut; i++){ // hvetoF[i]->SetLineColor(i+2); // hvetoF[i]->Draw(i == 0 ? "" : "same"); //} padID++; cDecay->cd(padID); hIonsPos->Draw("colz"); padID++; cDecay->cd(padID); hBetaPos->Draw("colz"); }