// reconstruction of pi0 & eta //last changed on 2013-08-30 class RhoCandList; class RhoCandidate; class PndAnaPidSelector; class PndAnaPidCombiner; class PndAnalysis; void ana_pi0(int nevts=0) { TStopwatch timer; timer.Start(); // *** some variables int i=0,j=0, k=0; TString OutFile="run_outfile.root"; // *** the files coming from the simulation TString inPidFile = "pid_complete.root"; // this file contains the PndPidCandidates and McTruth TString inParFile = "simparams.root"; gStyle->SetOptFit(1011); gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C"); FairLogger::GetLogger()->SetLogToFile(kFALSE); // *** initialization FairRunAna* fRun = new FairRunAna(); FairRuntimeDb* rtdb = fRun->GetRuntimeDb(); fRun->SetInputFile(inPidFile); FairParRootFileIo* parIO = new FairParRootFileIo(); parIO->open(inParFile); rtdb->setFirstInput(parIO); rtdb->setOutput(parIO); fRun->SetOutputFile(OutFile); fRun->Init(); // *** create an output file for all histograms TFile *out = TFile::Open("1K_output_ana_hist.root","RECREATE"); // *** create some histograms //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ pi0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ TH1F *hpizm_mcf_bf = new TH1F("hpizm_mcf_bf","#pi_{0} mass distribution (before mass constraint fit); GeV/c^{2}; counts",200,0,0.28); TH1F *hpizm_chi2_mf = new TH1F("hpizm_chi2_mf","Mass constraint fit #chi^{2} of #pi_{0}; #chi^{2}; counts",200,0,0.15); TH1F *hpizm_prob = new TH1F("hpizm_prob","Chi2 probability distribution of #pi_{0} mass constraint fit; prob; counts",200,0,1); TH1F *hpizm_mcf = new TH1F("hpizm_mcf","#pi_{0} mass distribution (after mass constraint fit); GeV/c^{2}; counts",200,0,0.28); TH1F *hpizm_bargam = new TH1F("hpizm_bargam","#pi_{0} mass distribution (#gamma#gamma @barrel); GeV/c^{2}; counts",200,0,0.28); TH1F *hpizm_fwdgam = new TH1F("hpizm_fwdgam","#pi_{0} mass distribution (#gamma#gamma @forward end-cap); GeV/c^{2}; counts",200,0,0.28); TH1F *hpizm_bwdgam = new TH1F("hpizm_bwdgam","#pi_{0} mass distribution (#gamma#gamma @backward end-cap); GeV/c^{2}; counts",200,0,0.28); TH1F *hpizm_bfgam = new TH1F("hpizm_bfgam","#pi_{0} mass distribution (#gamma_{1} @barrel & #gamma_{2} @forward end-cap); GeV/c^{2}; counts",200,0,0.28); TH1F *hpizm_bwbgam = new TH1F("hpizm_bwbgam","#pi_{0} mass distribution (#gamma_{1} @backward end-cap & #gamma_{2} @barrel); GeV/c^{2}; counts",200,0,0.28); TH1F *hpizm_bwfgam = new TH1F("hpizm_bwfgam","#pi_{0} mass distribution (#gamma_{1} @backward end-cap & #gamma_{2} @forward end-cap); GeV/c^{2}; counts",200,0,0.28); //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~gamma~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ TH1F *hrawgam_mass = new TH1F("hrawgam_mass","Photon Candidate Mass Distribution; M_{#gamma} (GeV/c^{2}); counts",400,-0.1,2); TH2F *hrawgame_theta_sim = new TH2F("hrawgame_theta_sim","Photon energy distribution vs. #theta;#theta degree;E_{#gamma} GeV",200,-1,181,200,-1,5); TH2F *hgame_theta_sim_mc = new TH2F("hgame_theta_sim_mc","McTruth photon energy distribution vs. #theta ;#theta degree;E_{#gamma} GeV",200,-1,181,200,-1,5); TH1F *hgame_sim_all = new TH1F("hgame_sim_all","Photon energy distribution (all); E_{#gamma} (GeV); counts",200,0,1); TH1F *hgame_sim_bar = new TH1F("hgame_sim_bar","Photon energy distribution (barrel EMC); E_{#gamma} (GeV); counts",200,0,0.5); TH1F *hgame_sim_fwd = new TH1F("hgame_sim_fwd","Photon energy distribution (forward end-cap EMC); E_{#gamma} (GeV); counts",200,0,1.5); TH1F *hgame_sim_bwd = new TH1F("hgame_sim_bwd","Photon energy distribution (backward end-cap EMC); E_{#gamma} (GeV); counts",200,0,0.4); TH1F *hgame_sim_all_mc = new TH1F("hgame_sim_all_mc","McTruth photon energy distribution (all); E_{#gamma} (GeV); counts",200,0,2); TH1F *hgame_sim_bar_mc = new TH1F("hgame_sim_bar_mc","McTruth photon energy distribution (barrel EMC); E_{#gamma} (GeV); counts",200,0,0.8); TH1F *hgame_sim_fwd_mc = new TH1F("hgame_sim_fwd_mc","McTruth photon energy distribution (forward end-cap EMC); E_{#gamma} (GeV); counts",200,0,3); TH1F *hgame_sim_bwd_mc = new TH1F("hgame_sim_bwd_mc","McTruth photon energy distribution (backward end-cap EMC); E_{#gamma} (GeV); counts",200,0,0.4); // *** the data reader object PndAnalysis* theAnalysis = new PndAnalysis(); if (nevts==0) nevts= theAnalysis->GetEntries(); // *** TCandLists for the analysis RhoCandList mctruth, pp, eplus, pipbase, pimbase, pip, pim,rawgambase, rawgam, bar_gam, fwd_gam, bwd_gam; RhoCandList pi0_bargam, pi0_fwdgam, pi0_bwdgam, pi0_bfgam, pi0_bwfgam, pi0_bwbgam, pi0_raw, pi0_reco; RhoCandList twopi, twopi_reco, eta_raw, eta_reco; //define TLorentzVector Double_t wmms,mms,eds,eeta,eeplus,enu; TLorentzVector pbeam, pds,pkplus,pkminus,ppids,peta,peplus,pnu,ppipe,ppime,ppize; TLorentzVector Pkppi,Pkmpi,Petae,Pnue,Ppmpz,Ppppz,pbp; //pbarpSystem TLorentzVector pt(0,0,0,0.938272); TLorentzVector pbeam(0,0,8.0,8.0548); // PDG mass (GeV/c^2) double pdgmass_eta = 0.548; double pdgmass_pi0 = 0.135; double pdgmass_gamma = 0.00001; // *** Mass selector for cands RhoMassParticleSelector *etaMassSel=new RhoMassParticleSelector("etaSelctor",pdgmass_eta,pdgmass_eta*0.5); RhoMassParticleSelector *pi0MassSel=new RhoMassParticleSelector("pi0Selctor",pdgmass_pi0,pdgmass_pi0*2); RhoMassParticleSelector *pi0_fineMassSel=new RhoMassParticleSelector("pi0_fSelctor",pdgmass_pi0,pdgmass_pi0*0.5); RhoMassParticleSelector *gammaMassSel=new RhoMassParticleSelector("gSelctor",pdgmass_gamma,pdgmass_gamma); // *** // the event loop // *** while (theAnalysis->GetEvent() && i++FillList(mctruth,"McTruth"); theAnalysis->FillList(pipbase, "PionLoosePlus"); theAnalysis->FillList(pimbase, "PionLooseMinus"); theAnalysis->FillList(rawgambase,"Neutral"); pip.Cleanup(); pim.Cleanup(); //rawgambase.Select(gammaMassSel); rawgam.Cleanup(); // use Monte Carlo PID information to clean particle lists for (j=0; jPdgCode() == 211)pip.Add(pipbase[j]); } for (j=0; jPdgCode() == -211) pim.Add(pimbase[j]); } for (j=0; jPdgCode() == 22) rawgam.Add(rawgambase[j]); rawgam.Add(rawgambase[j]); hrawgam_mass->Fill(rawgambase[j]->M()); TLorentzVector rawGamma = rawgambase[j]->P4(); double degree_raw = 180*rawGamma.Theta()/3.141592654; hrawgame_theta_sim->Fill(degree_raw ,rawGamma.E()); } for (j=0;jPdgCode(); RhoCandidate *mcmother = mctruth[j]->TheMother(); if (mcmother) int pdg_mother = mcmother->PdgCode(); // eta from Ds+ if(mcmother && pdg_mother==431 && pdg_mc==221 ) { RhoCandidate *trueeta = mctruth[j]; double trueetaM = trueeta->M(); TLorentzVector tlveta = trueeta->P4(); TVector3 trueetaMom = trueeta->P(); TVector3 trueetaVtx = trueeta->Pos(); } // gamma from pi0 if(mcmother && pdg_mother==111 && pdg_mc==22 ) { RhoCandidate *gam_mc = mctruth[j]; TLorentzVector mcGamma = mctruth[j]->P4(); hgame_sim_all_mc->Fill(mctruth[j]->E()); double gam_theta_mc = 180*mcGamma.Theta()/3.141592654; hgame_theta_sim_mc->Fill(gam_theta_mc,mcGamma.E()); if (gam_theta_mc>0.35 && gam_theta_mc<2.4) // @ EMC barrel { hgame_sim_bar_mc->Fill(mctruth[j]->E()); } if (gam_theta_mc<0.35) // @ EMC forward end-cup { hgame_sim_fwd_mc->Fill(mctruth[j]->E()); } if (gam_theta_mc>2.4) // @ EMC backward end-cup { hgame_sim_bwd_mc->Fill(mctruth[j]->E()); } } } //-------------------------------- reconstruction ---------------------------------------------- ///////////////////////////////////// pi0 //////////////////////////////////////////////////////////// //`````````````````````````````````````select gammas``````````````````````` bar_gam.Cleanup(); fwd_gam.Cleanup(); for (j=0;jP4(); double gam_theta = tlvGamma.Theta(); double gam_e = tlvGamma.E(); hgame_sim_all->Fill(gam_e); if (gam_theta>0.35 && gam_theta<2.4) // @ EMC barrel { hgame_sim_bar->Fill(gam_e); if (gam_e>0.1)bar_gam.Add( rawgam[j] ); } if (gam_theta<0.35) // @ EMC forward end-cup { hgame_sim_fwd->Fill(gam_e); if (gam_e>0.3)fwd_gam.Add( rawgam[j] ); } if (gam_theta>2.4) // @ EMC backward end-cup { hgame_sim_bwd->Fill(gam_e); if (gam_e>0.1)bwd_gam.Add( rawgam[j] ); } }//`````````````````````````````````````````````````````````````````` // pi0 -> gamma gamma pi0_bargam.Cleanup(); pi0_fwdgam.Cleanup(); pi0_bwdgam.Cleanup(); pi0_bfgam.Cleanup(); pi0_bwbgam.Cleanup(); pi0_bwfgam.Cleanup(); pi0_bargam.Combine(bar_gam,bar_gam); pi0_fwdgam.Combine(fwd_gam,fwd_gam); pi0_bwdgam.Combine(bwd_gam,bwd_gam); pi0_bfgam.Combine(bar_gam,fwd_gam); pi0_bwfgam.Combine(bwd_gam,fwd_gam); pi0_bwbgam.Combine(bar_gam,bwd_gam); pi0_bargam.Select(pi0MassSel); pi0_fwdgam.Select(pi0MassSel); pi0_bwdgam.Select(pi0MassSel); pi0_bfgam.Select(pi0MassSel); pi0_bwfgam.Select(pi0MassSel); pi0_bwbgam.Select(pi0MassSel); for (j=0;jFill(pi0_bargam[j]->M()); } for (j=0;jFill(pi0_fwdgam[j]->M()); } for (j=0;jFill(pi0_bwdgam[j]->M()); } for (j=0;jFill(pi0_bfgam[j]->M()); } for (j=0;jFill(pi0_bwfgam[j]->M()); } for (j=0;jFill(pi0_bwbgam[j]->M()); } // sum up the three components to raw pi0 pi0_raw.Cleanup(); pi0_bargam.Select(pi0_fineMassSel); pi0_fwdgam.Select(pi0_fineMassSel); pi0_bwdgam.Select(pi0_fineMassSel); pi0_bfgam.Select(pi0_fineMassSel); pi0_bwfgam.Select(pi0_fineMassSel); pi0_bwbgam.Select(pi0_fineMassSel); pi0_raw.Append(pi0_bargam); pi0_raw.Append(pi0_fwdgam); pi0_raw.Append(pi0_bfgam); pi0_raw.Append(pi0_bwdgam); pi0_raw.Append(pi0_bwbgam); pi0_raw.Append(pi0_bwfgam); // do mass constraint fitting for pi0 pi0_reco.Cleanup(); for (j=0;jFill(pi0_raw[j]->M()); // mass of raw pi0 PndKinFitter mfitter(pi0_raw[j]); mfitter.AddMassConstraint(pdgmass_pi0); mfitter.Fit(); double chi2_pi_m = mfitter.GetChi2(); hpizm_chi2_mf->Fill(chi2_pi_m); int ndf_m = 1; double prob = TMath::Prob(chi2_pi_m,ndf_m); hpizm_prob->Fill(prob); //if (prob>0.6 && prob<0.8) if (chi2_pi_m>0 && chi2_pi_m<0.0001) { hpizm_mcf->Fill(pi0_raw[j]->M()); pi0_reco.Add( pi0_raw[j] ); } } //pi0_reco.SetType(111); }//end event loop // *** write out all the histos out->cd(); hpizm_mcf_bf->Write(); hpizm_chi2_mf->Write(); hpizm_prob->Write(); hpizm_mcf->Write(); hpizm_bargam->Write(); hpizm_fwdgam->Write(); hpizm_bwdgam->Write(); hpizm_bfgam->Write(); hpizm_bwfgam->Write(); hpizm_bwbgam->Write(); hgame_sim_all->Write(); hgame_sim_bar->Write(); hgame_sim_fwd->Write(); hgame_sim_bwd->Write(); hgame_sim_all_mc->Write(); hgame_sim_bar_mc->Write(); hgame_sim_fwd_mc->Write(); hgame_sim_bwd_mc->Write(); hrawgam_mass->Write(); hrawgame_theta_sim->Write(); hgame_theta_sim_mc->Write(); out->Save(); timer.Stop(); printf("\n ************** Job finished successfully! ^-^ *************\n"); Double_t rtime = timer.RealTime(); Double_t ctime = timer.CpuTime(); printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime); }//end of macro