class RhoCandList; class RhoCandidate; class PndAnaPidSelector; class PndAnaPidCombiner; class PndAnalysis; // *** routine to only keep PID matched candidates in list int SelectTruePid(PndAnalysis *ana, RhoCandList &l) { int removed = 0; for (int ii=l.GetLength()-1;ii>=0;--ii) { if ( !(ana->McTruthMatch(l[ii])) ) { l.Remove(l[ii]); removed++; } } return removed; } void printCand(RhoCandidate *c) { TLorentzVector lv=c->P4(); cout <PdgCode()<<" ("<P4() - l[j]->P4(); bool chkmc = (l[i]->GetMcTruth()==l[j]->GetMcTruth()); bool chktrk = (fabs(dl.X())AddParticle("pbarpSystem","pbarpSystem",1.9,kFALSE,0.1,0,"",88888); TStopwatch fTimer; // *** some variables int i=0,j=0, k=0, l=0; gStyle->SetOptFit(1011); // *** the output file for FairRunAna TString OutFile="output.root"; // *** the files coming from the simulation TString inPidFile = "dpm_qa_pid.root"; // this file contains the PndPidCandidates and McTruth TString inParFile = "dpm_qa_par.root"; // *** PID table with selection thresholds; can be modified by the user TString pidParFile = TString(gSystem->Getenv("VMCWORKDIR"))+"/macro/params/all.par"; // *** initialization FairLogger::GetLogger()->SetLogToFile(kFALSE); FairRunAna* fRun = new FairRunAna(); FairRuntimeDb* rtdb = fRun->GetRuntimeDb(); fRun->SetInputFile(inPidFile); // *** setup parameter database FairParRootFileIo* parIO = new FairParRootFileIo(); parIO->open(inParFile); FairParAsciiFileIo* parIOPid = new FairParAsciiFileIo(); parIOPid->open(pidParFile.Data(),"in"); rtdb->setFirstInput(parIO); rtdb->setSecondInput(parIOPid); rtdb->setOutput(parIO); fRun->SetOutputFile(OutFile); fRun->Init(); // *** create an output file for all histograms TFile *out = TFile::Open("output_ana.root","RECREATE"); // *** create some histograms TH1F *hmomtrk = new TH1F("hmomtrk","track momentum (all)",200,0,5); TH1F *hthttrk = new TH1F("hthttrk","track theta (all)",200,0,3.1415); TH1F *hjpsim_all = new TH1F("hjpsim_all","J/#psi mass (all)",200,0,4.5); TH1F *hpsim_all = new TH1F("hpsim_all","#psi(2S) mass (all)",200,0,5); TH1F *hjpsim_lpid = new TH1F("hjpsim_lpid","J/#psi mass (loose pid)",200,0,4.5); TH1F *hpsim_lpid = new TH1F("hpsim_lpid","#psi(2S) mass (loose pid)",200,0,5); TH1F *hjpsim_tpid = new TH1F("hjpsim_tpid","J/#psi mass (tight pid)",200,0,4.5); TH1F *hpsim_tpid = new TH1F("hpsim_tpid","#psi(2S) mass (tight pid)",200,0,5); TH1F *hjpsim_trpid = new TH1F("hjpsim_trpid","J/#psi mass (true pid)",200,0,4.5); TH1F *hpsim_trpid = new TH1F("hpsim_trpid","#psi(2S) mass (true pid)",200,0,5); TH1F *hjpsim_ftm = new TH1F("hjpsim_ftm","J/#psi mass (full truth match)",200,0,4.5); TH1F *hpsim_ftm = new TH1F("hpsim_ftm","#psi(2S) mass (full truth match)",200,0,5); TH1F *hjpsim_nm = new TH1F("hjpsim_nm","J/#psi mass (no truth match)",200,0,4.5); TH1F *hpsim_nm = new TH1F("hpsim_nm","#psi(2S) mass (no truth match)",200,0,5); TH1F *hjpsim_diff = new TH1F("hjpsim_diff","J/#psi mass diff to truth",100,-2,2); TH1F *hpsim_diff = new TH1F("hpsim_diff","#psi(2S) mass diff to truth",100,-2,2); TH1F *hjpsim_vf = new TH1F("hjpsim_vf","J/#psi mass (vertex fit)",200,0,4.5); TH1F *hjpsim_4cf = new TH1F("hjpsim_4cf","J/#psi mass (4C fit)",200,0,4.5); TH1F *hjpsim_mcf = new TH1F("hjpsim_mcf","J/#psi mass (mass constraint fit)",200,0,4.5); TH1F *hjpsi_chi2_vf = new TH1F("hjpsi_chi2_vf", "J/#psi: #chi^{2} vertex fit",100,0,10); TH1F *hpsi_chi2_4c = new TH1F("hpsi_chi2_4c", "#psi(2S): #chi^{2} 4C fit",100,0,250); TH1F *hjpsi_chi2_mf = new TH1F("hjpsi_chi2_mf", "J/#psi: #chi^{2} mass fit",100,0,10); TH1F *hjpsi_prob_vf = new TH1F("hjpsi_prob_vf", "J/#psi: Prob vertex fit",100,0,1); TH1F *hpsi_prob_4c = new TH1F("hpsi_prob_4c", "#psi(2S): Prob 4C fit",100,0,1); TH1F *hjpsi_prob_mf = new TH1F("hjpsi_prob_mf", "J/#psi: Prob mass fit",100,0,1); TH2F *hvpos = new TH2F("hvpos","(x,y) projection of fitted decay vertex",100,-2,2,100,-2,2); // // Now the analysis stuff comes... // // *** the data reader object PndAnalysis* theAnalysis = new PndAnalysis(); if (nevts==0) nevts= theAnalysis->GetEntries(); // *** RhoCandLists for the analysis RhoCandList chrg, muplus, muminus, piplus, piminus, jpsi, psi2s; // *** Mass selector for the jpsi cands double m0_jpsi = TDatabasePDG::Instance()->GetParticle("J/psi")->Mass(); // Get nominal PDG mass of the J/psi RhoMassParticleSelector *jpsiMassSel=new RhoMassParticleSelector("jpsi",m0_jpsi,1.0); // *** the lorentz vector of the initial psi(2S) TLorentzVector ini(0, 0, 6.231552, 7.240065); // *** // the event loop // *** int cntdbltrk=0, cntdblmc=0, cntdblboth=0, cnttrk=0, cnt_dbl_jpsi=0, cnt_dbl_psip=0; while (theAnalysis->GetEvent() && i++FillList(chrg, "Charged"); theAnalysis->FillList(muplus, "MuonAllPlus"); theAnalysis->FillList(muminus, "MuonAllMinus"); theAnalysis->FillList(piplus, "PionAllPlus"); theAnalysis->FillList(piminus, "PionAllMinus"); // *** momentum and theta histograms for (j=0;jFill(muplus[j]->P()); hthttrk->Fill(muplus[j]->P4().Theta()); } for (j=0;jFill(muminus[j]->P()); hthttrk->Fill(muminus[j]->P4().Theta()); } cnttrk += chrg.GetLength(); int n1, n2, n3; countDoubles(chrg,n1,n2,n3); cntdbltrk += n1; cntdblmc += n2; cntdblboth += n3; // *** combinatorics for J/psi -> mu+ mu- jpsi.Combine(muplus, muminus); // *** // *** do the TRUTH MATCH for jpsi // *** jpsi.SetType(443); int nm = 0; for (j=0;jFill( jpsi[j]->M() ); if (theAnalysis->McTruthMatch(jpsi[j])) { nm++; hjpsim_ftm->Fill( jpsi[j]->M() ); hjpsim_diff->Fill( jpsi[j]->GetMcTruth()->M() - jpsi[j]->M() ); } else hjpsim_nm->Fill( jpsi[j]->M() ); } if (nm>1) cnt_dbl_jpsi++; // *** // *** do VERTEX FIT (J/psi) // *** for (j=0;jFill(chi2_vtx); hjpsi_prob_vf->Fill(prob_vtx); if ( prob_vtx > 0.01 ) // when good enough, fill some histos { RhoCandidate *jfit = jpsi[j]->GetFit(); // access the fitted cand TVector3 jVtx=jfit->Pos(); // and the decay vertex position hjpsim_vf->Fill(jfit->M()); hvpos->Fill(jVtx.X(),jVtx.Y()); } } // *** some rough mass selection jpsi.Select(jpsiMassSel); // *** combinatorics for psi(2S) -> J/psi pi+ pi- psi2s.Combine(jpsi, piplus, piminus); // *** // *** do the TRUTH MATCH for psi(2S) // *** psi2s.SetType(88888); nm = 0; for (j=0;jFill( psi2s[j]->M() ); if (theAnalysis->McTruthMatch(psi2s[j])) { nm++; hpsim_ftm->Fill( psi2s[j]->M() ); hpsim_diff->Fill( psi2s[j]->GetMcTruth()->M() - psi2s[j]->M() ); } else hpsim_nm->Fill( psi2s[j]->M() ); } if (nm>1) cnt_dbl_psip++; // *** // *** do 4C FIT (initial psi(2S) system) // *** for (j=0;jFill(chi2_4c); hpsi_prob_4c->Fill(prob_4c); if ( prob_4c > 0.01 ) // when good enough, fill some histo { RhoCandidate *jfit = psi2s[j]->Daughter(0)->GetFit(); // get fitted J/psi hjpsim_4cf->Fill(jfit->M()); } } // *** // *** do MASS CONSTRAINT FIT (J/psi) // *** for (j=0;jFill(chi2_m); hjpsi_prob_mf->Fill(prob_m); if ( prob_m > 0.01 ) // when good enough, fill some histo { RhoCandidate *jfit = jpsi[j]->GetFit(); // access the fitted cand hjpsim_mcf->Fill(jfit->M()); } } // *** // *** TRUE PID combinatorics // *** // *** do MC truth match for PID type SelectTruePid(theAnalysis, muplus); SelectTruePid(theAnalysis, muminus); SelectTruePid(theAnalysis, piplus); SelectTruePid(theAnalysis, piminus); // *** all combinatorics again with true PID jpsi.Combine(muplus, muminus); for (j=0;jFill( jpsi[j]->M() ); jpsi.Select(jpsiMassSel); psi2s.Combine(jpsi, piplus, piminus); for (j=0;jFill( psi2s[j]->M() ); // *** // *** LOOSE PID combinatorics // *** // *** and again with PidAlgoMvd;PidAlgoStt;PidAlgoDrc and loose selection theAnalysis->FillList(muplus, "MuonLoosePlus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc;PidAlgoDisc;PidAlgoMdtHardCuts"); theAnalysis->FillList(muminus, "MuonLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc;PidAlgoDisc;PidAlgoMdtHardCuts"); theAnalysis->FillList(piplus, "PionLoosePlus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc;PidAlgoDisc"); theAnalysis->FillList(piminus, "PionLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc;PidAlgoDisc"); jpsi.Combine(muplus, muminus); for (j=0;jFill( jpsi[j]->M() ); jpsi.Select(jpsiMassSel); psi2s.Combine(jpsi, piplus, piminus); for (j=0;jFill( psi2s[j]->M() ); // *** // *** TIGHT PID combinatorics // *** // *** and again with PidAlgoMvd;PidAlgoStt and tight selection theAnalysis->FillList(muplus, "MuonTightPlus", "PidAlgoMdtHardCuts"); theAnalysis->FillList(muminus, "MuonTightMinus", "PidAlgoMdtHardCuts"); theAnalysis->FillList(piplus, "PionLoosePlus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc"); theAnalysis->FillList(piminus, "PionLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc"); jpsi.Combine(muplus, muminus); for (j=0;jFill( jpsi[j]->M() ); jpsi.Select(jpsiMassSel); psi2s.Combine(jpsi, piplus, piminus); for (j=0;jFill( psi2s[j]->M() ); } // *** write out all the histos out->cd(); hmomtrk->Write(); hthttrk->Write(); hjpsim_all->Write(); hpsim_all->Write(); hjpsim_lpid->Write(); hpsim_lpid->Write(); hjpsim_tpid->Write(); hpsim_tpid->Write(); hjpsim_trpid->Write(); hpsim_trpid->Write(); hjpsim_ftm->Write(); hpsim_ftm->Write(); hjpsim_nm->Write(); hpsim_nm->Write(); hpsim_diff->Write(); hjpsim_diff->Write(); hjpsim_vf->Write(); hjpsim_4cf->Write(); hjpsim_mcf->Write(); hjpsi_chi2_vf->Write(); hpsi_chi2_4c->Write(); hjpsi_chi2_mf->Write(); hjpsi_prob_vf->Write(); hpsi_prob_4c->Write(); hjpsi_prob_mf->Write(); hvpos->Write(); out->Save(); // Extract the maximal used memory an add is as Dart measurement // This line is filtered by CTest and the value send to CDash FairSystemInfo sysInfo; Float_t maxMemory=sysInfo.GetMaxMemory(); cout << ""; cout << maxMemory; cout << "" << endl; fTimer.Stop(); Double_t rtime = fTimer.RealTime(); Double_t ctime = fTimer.CpuTime(); Float_t cpuUsage=ctime/rtime; cout << ""; cout << cpuUsage; cout << "" << endl; cout << endl; cout << "Real time " << rtime << " s, CPU time " << ctime << "s" << endl; cout << "CPU usage " << cpuUsage*100. << "%" << endl; cout << "Max Memory " << maxMemory << " MB" << endl; cout << "Macro finished successfully." << endl; return 0; }