#include "TFile.h" #include "TH1F.h" #include "TH2F.h" #include "TCanvas.h" #include "TStopwatch.h" #include "TROOT.h" #include "TSystem.h" #include "TTree.h" #include "TString.h" /* #include "RhoSelector/TPidSelector.h" #include "RhoBase/TCandList.h" #include "RhoBase/TCandListIterator.h" #include "RhoBase/TCandidate.h" #include "RhoBase/TFactory.h" */ void config_histo(TH1 *h, TString tx, TString ty) { h->SetLineWidth(2); //h->SetLineColor(1); //h->SetFillColor(1); h->GetXaxis()->SetTitleOffset(1.2); h->GetXaxis()->SetTitleColor(1); h->GetXaxis()->SetLabelSize(0.05); h->GetXaxis()->SetTitleSize(0.05); h->GetXaxis()->SetNdivisions(505); h->GetYaxis()->SetTitleOffset(1.6); h->GetYaxis()->SetTitleFont(42); h->GetYaxis()->SetLabelSize(0.05); h->GetYaxis()->SetTitleSize(0.05); h->SetMinimum(0); h->SetFillStyle(1000); h->SetFillColor(3); h->SetLineWidth(1); h->SetXTitle(tx); h->SetYTitle(ty); } void printCand(TLorentzVector l, TVector3 p) { cout <<"vtx("<LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C"); basiclibs(); // **** the library is needed in order to use the TCandidate and TCandList etc gSystem->Load("libRho"); TCanvas *c1=new TCanvas("c1","c1",1000,700); c1->Divide(3,2); // **** access the datastructure holding the particle lists // TFile* f = new TFile(fname.Data()); TTree *t=(TTree*)f->Get("cbmsim") ; // **** for every event, a TCLonesArray with the candidates is stored in cbmsim // TClonesArray *fChrgCands=new TClonesArray("TCandidate"); TClonesArray *fNeutCands=new TClonesArray("TCandidate"); TClonesArray *fMcCands=new TClonesArray("TCandidate"); t->SetBranchAddress("PndChargedCandidates",&fChrgCands); t->SetBranchAddress("PndNeutralCandidates",&fNeutCands); t->SetBranchAddress("PndMcTracks",&fMcCands); TCandidate *tc; // **** create and setup some histos for QA plots // TH1F *phimass = new TH1F("phimass","phi cands",100,0.98,1.1); TH1F *pi0mass = new TH1F("pi0mass","pi0 cands",100,0.130-0.05,0.130+0.05); TH1F *dsmass = new TH1F("dsmass","Ds cands",100,1.968-0.06,1.968+0.06); TH1F *ds0mass = new TH1F("ds0mass","Ds0 cands",100,2.317-0.1,2.317+0.1); TH1F *ppmass = new TH1F("ppmass","pbarp cands",100,4.291-0.1,4.291+0.1); TH1F *pp2mass = new TH1F("pp2mass","sum m_{miss}+m_{Ds}",100,4.2855-0.02,4.2855+0.02); TH1F *nmult=new TH1F("nmult","# neutrals",15,0,15); TH1F *pdiff=new TH1F("pdiff","momentum difference",100,0,2); TH2F *pdiff2=new TH2F("pdiff2","gamma corr",50,0,2.5,50,0,2.5); config_histo(phimass,"m(K^{+} K^{-}) [GeV/c^{2}]","entries"); config_histo(pi0mass,"m(#gamma#gamma) [GeV/c^{2}]","entries"); config_histo(dsmass,"m(#phi#pi^{#pm}) [GeV/c^{2}]","entries"); config_histo(ds0mass,"m(D_{s}#pi^{0}) [GeV/c^{2}]","entries"); config_histo(ppmass,"m(D_{s} D_{s0}*) [GeV/c^{2}]","entries"); config_histo(pp2mass,"m_{miss}+m_{D_{s}} [GeV/c^{2}]","entries"); /*phimass->SetMinimum(0); pi0mass->SetMinimum(0); dsmass->SetMinimum(0); ds0mass->SetMinimum(0); ppmass->SetMinimum(0); */ if (num==0) num= t->GetEntriesFast(); cout <<"\n####### Processing "<SetCriterion("veryLoose"); TPidSimplePionSelector *piSel = new TPidSimplePionSelector(); piSel->SetCriterion("veryLoose"); TH2F *hpro=new TH2F("hpro","",300,-1,1,300,-1,1); int i1,i2; // **** loop over all _events_ // for (Int_t j=0; j< num;j++) { if ((j%100)==0) cout <<"evt "<GetEntry(j); TFactory::Instance()->Reset(); mcCands.Cleanup(); chargedCands.Cleanup(); chargedCands2.Cleanup(); neutralCands.Cleanup(); // **** loop over all Candidates and add them to the list allCands // for (i1=0; i1GetEntriesFast(); i1++){ tc = (TCandidate *)fChrgCands->At(i1); TLorentzVector l=tc->P4(); TVector3 p=tc->Pos(); if (pr!=0) { cout <<"1 --> "; cout <<"c("<Charge()<<") vtx("< "; cout <<"c("<Charge()<<") vtx("<"<<*tc<Fill(p.Mag()); chargedCands2.Add(*tc); } for (i1=0; i1GetEntriesFast(); i1++){ tc = (TCandidate *)fNeutCands->At(i1); double minang=10; for (i2=0;i2P3(); TVector3 chrg=chargedCands[i2].Pos(); double ang=gam.Angle(chrg); //pdiff2->Fill(gam.Angle(chrg), tc->P4().E()); if (minang>ang) minang=ang; } pdiff2->Fill(minang, tc->P4().E()); pdiff->Fill(tc->P4().E()); //tc->SetEnergy(tc->P4().E()*1.07); if (minang>0.3) neutralCands.Add(*tc); } for (i1=0; i1GetEntriesFast(); i1++){ tc = (TCandidate *)fMcCands->At(i1); mcCands.Add(*tc); } // **** select all the basic list // cout <<"c:"<Fill(neutralCands.GetLength()); plusCands.Select(chargedCands2 ,plusSel); minusCands.Select(chargedCands2 ,minusSel); // **** pid selection // kpCands.Select(plusCands ,kSel); kmCands.Select(minusCands ,kSel); piCands.Select(chargedCands2 ,piSel); cout <<"chrg:"<Fill(tc->M()); phiCands.Select(phiMSel); cout <<"phi:"<Fill(tc->M()); TLorentzVector lds=tc->P4(); TLorentzVector ini(0,0,8.824,9.812); TLorentzVector miss=ini-lds; //ds0mass->Fill(miss.M()); pp2mass->Fill(miss.M()+lds.M()); } dsCands.Select(dsMSel); pi0Cands.Combine(neutralCands,neutralCands); TCandListIterator iterPi0(pi0Cands); while (tc=iterPi0.Next()) pi0mass->Fill(tc->M()); pi0Cands.Select(pi0MSel); ds0Cands.Combine(dsCands,pi0Cands); TCandListIterator iterDs0(ds0Cands); while (tc=iterDs0.Next()) ds0mass->Fill(tc->M()); ppCands.Combine(ds0Cands,dsCands); //cout << ds0Cands.GetLength()<<" "<Fill(tc->M()); //printTree(tc); } TCandListIterator iterChrg(chargedCands); while(tc=iterChrg.Next()) { int idx=tc->GetMcIdx(); //if (idx!=-1) pdiff->Fill(tc->P()-mcCands[idx]->P()); } } // **** plot all that stuff // c1->cd(1); phimass->Draw(); c1->cd(2); pi0mass->Draw(); c1->cd(3); dsmass->Draw(); c1->cd(4); ds0mass->Draw(); c1->cd(5); ppmass->Draw(); //pdiff->Draw(); c1->cd(6); pp2mass->Draw(); //pdiff2->Draw("colz"); //hpro->Draw(); c1->cd(); //timer.Stop(); //Double_t rtime = timer.RealTime(); //Double_t ctime = timer.CpuTime(); //printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime); } /*void printTree(TCandidate *tc, int level=0) { int i=0; int nd=tc->NDaughters(); if (nd==0) return; cout <Uid()<<"("< "; for (i=0;iDaughter(i)->Uid()<<" "; cout <Daughter(i),level+1); if (level==0) cout <