// reconstruction of Ds- //created in 2013-04-10 (/private/workspace/0418_subjobs/2_2K) //modified in 2013-08-02 to adapt new version of pandaroot // in /private/mymacro/dsminus //last changed on 2013-09-13 class RhoCandList; class RhoCandidate; class PndAnaPidSelector; class PndAnaPidCombiner; class PndAnalysis; // *** plot ctau histogarm void ctauhistogramm(TH1F* histo, TString savefile) { histo->SetLineColor(kBlack); histo->Draw(""); TF1 *fExpo = new TF1("fExpo","expo",0,100); fExpo->SetLineColor(kRed); fExpo->SetLineWidth(2); histo->Fit("fExpo","RQ"); printf("Slope %f", fExpo->GetParameter("Slope")); fExpo->GetParameter("Slope"); TPaveText *pt = new TPaveText(44.5781,59.32019,104.1472,206.744,""); pt->UseCurrentStyle(); pt->SetFillColor(kWhite); pt->SetBorderSize(1); TString myString; std::stringstream virtualString; virtualString<< "c#tau = "; virtualString<< TMath::Abs(1/fExpo->GetParameter("Slope")); virtualString<< "#pm"; virtualString<GetParameter("Slope"))*TMath::Abs(1/fExpo->GetParameter("Slope"))*fExpo->GetParError(fExpo->GetParNumber("Slope")); virtualString<< " cm"; TString bTagisGreaterThanSTRING = virtualString.str(); pt->AddText(bTagisGreaterThanSTRING.Data() ); pt->Draw(); gPad->SetLogy(); //gPad->Print(savefile,"pdf"); // nicht fuer powerpoint geeignet gPad->Print(savefile,"eps"); // gPad->Print("result_pictures/fMeanLife.gif","gif"); // ganz gut, aber keine Vektorgrafik // gPad->Print("result_pictures/fMeanLife.png","png"); // ganz gut, aber keine Vektorgrafik // gPad->Print("result_pictures/fMeanLife.jpg","jpg"); // schlechte qualitaet } void mytask(int nevts= 0) { TStopwatch timer; timer.Start(); // *** some variables int i=0,j=0, k=0, l=0; TString OutFile="output.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(); /*TFile mcfile(inPidFile, "READ"); TTree* mctree = (TTree*)mcfile.Get("cbmsim"); TClonesArray* mcarray = new TClonesArray("PndMCTrack"); mctree->SetBranchAddress("MCTrack", &mcarray); */ // *** create an output file for all histograms TFile *out = TFile::Open("1K_output_ana_hist.root","RECREATE"); // *** create some histograms //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Ds- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ TH1F *hdsmm_mctrk = new TH1F("hdsmm_mctrk","Ds^{-} mass distribution (MC track); GeV/c^{2}; counts",200,1.6,2.4); TH1F *hdsmm_before = new TH1F("hdsmm_before","Ds^{-} mass distribution (Before Fitting); GeV/c^{2}; counts",200,1.6,2.4); //-------------vertex fit------------- TH1F *hdsmm_vf = new TH1F("hdsmm_vf","Ds^{-} mass distribution after vertex fit; GeV/c^{2}; counts",200,1.6,2.4); TH1F *hdsm_vx_vf = new TH1F("hdsm_vx_vf","Ds^{-} vertex distribution (After Vertex Fit); X_{reco} cm; counts",200,-0.2,0.2); TH1F *hdsm_vtx_prob = new TH1F("hdsm_vtx_prob","Chi2 probability distribution of Ds^{-} vtx fit; prob; counts",200,0,1); TH1F *hdsm_chi2_vf = new TH1F("hdsm_chi2_vf","Vertex fit #chi^{2} of Ds^{-}; #chi^{2}; counts",200,0,20); //TH1F *hdsm_chi2_vf_mc = new TH1F("hdsm_chi2_vf_mc","Vertex fit #chi^{2} of MC true Ds^{-}; #chi^{2}; counts",200,0,80); // TCanvas *cvdsm_chi2_vf = new TCanvas("cvdsm_chi2_vf","Vertex fit #chi^{2} of Ds^{-}; #chi^{2}; counts",600,400); //--------------mass constraint fit--------- TH1F *hdsmm_mcf_bef = new TH1F("hdsmm_mcf_bef","Ds^{-} mass distribution (before mass constraint fit); GeV/c^{2}; counts",200,1.6,2.4); TH1F *hdsm_chi2_mf = new TH1F("hdsm_chi2_mf","Mass constraint fit #chi^{2} of Ds^{-}; #chi^{2}; counts",200,0,2); TH2F *hdsm_chi2vsmass = new TH2F("hdsm_chi2vsmass","Ds^{-} mass vs #chi^{2} of mass constraint fit;#chi^{2};mass GeV/c^{2}",200,-1,2,200,1,3); //TH1F *hdsm_chi2_mf_mc = new TH1F("hdsm_chi2_mf_mc","Mass constraint fit #chi^{2} of MC true Ds^{-}; #chi^{2}; counts",200,0,80); //TCanvas *cvdsm_chi2_mf = new TCanvas("cvdsm_chi2_mf","Mass constraint fit #chi^{2} of Ds^{-}; #chi^{2}; counts",600,400); TH1F *hdsmm_mcf_rej = new TH1F("hdsmm_mcf_rej","Ds^{-} rejected mass (mass constraint fit); GeV/c^{2}; counts",200,1.6,2.4); TH1F *hdsm_mcf_prob = new TH1F("hdsm_mcf_prob","Chi2 probability distribution of Ds^{-} mass constraint fit; prob; counts",200,0,1); TH1F *hdsmm_reco = new TH1F("hdsmm_reco","Ds^{-} mass distribution; GeV/c^{2}; counts",200,1.6,2.4); TH2F *hdvpos = new TH2F("hdvpos","(x,y) projection of fitted decay vertex location distribution of Ds^{-}-> K^{+} K^{-} #pi^{-};cm;cm",200,-1,1,200,-1,1); TH1F *hdsm_ptmc = new TH1F("hdsm_ptmc","Ds^{-} Pt in MC; Pt_{MC} GeV/c; counts",200,0,1); TH1F *hdsm_pzmc = new TH1F("hdsm_pzmc","Ds^{-} Pz in MC; Pz_{MC} GeV/c; counts",200,0,8); TH1F *hdsm_pt = new TH1F("hdsm_pt","Ds^{-} Pt distribution; Pt_{reco} GeV/c; counts",200,0,2); TH1F *hdsm_pz = new TH1F("hdsm_pz","Ds^{-} Pz distribution; Pt_{reco} GeV/c; counts",200,0,8); TH1F *hdsm_ptre_rltv = new TH1F("hdsm_ptre_rltv","Ds^{-} Pt relative resolution; #Delta Pt/Pt_{MC}; counts",200,-1,1); TH1F *hdsm_pzre_rltv = new TH1F("hdsm_pzre_rltv","Ds^{-} Pz relative resolution; #Delta Pz/Pz_{MC}; counts",200,-1,1); TH1F *hdsm_mre = new TH1F("hdsm_mre","Ds^{-} mass resolution; #delta M (M_{reco}-M_{MC}) GeV/c^{2}; counts",200,-1,1); TH1F *hdsm_vxre = new TH1F("hdsm_vxre","Ds^{-} Vertex_X location distribution; #Delta X (X_{reco}-X_{MC}) cm; counts",200,-0.2,0.2); TH1F *hdsm_vyre = new TH1F("hdsm_vyre","Ds^{-} Vertex_Y location distribution; #Delta Y (Y_{reco}-Y_{MC}) cm; counts",200,-0.2,0.2); TH1F *hdsm_vzre = new TH1F("hdsm_vzre","Ds^{-} Vertex_Z location distribution; #Delta Z (Z_{reco}-Z_{MC}) cm; counts",200,-0.2,0.2); TH1F *hdsm_vxmc = new TH1F("hdsm_vxmc","Ds^{-} MC Vertex_X location; X_{MC} cm; counts",200,-0.2,0.2); TH1F *hdsm_vymc = new TH1F("hdsm_vymc","Ds^{-} MC Vertex_Y location; Y_{MC} cm; counts",200,-0.2,0.2); TH1F *hdsm_vzmc = new TH1F("hdsm_vzmc","Ds^{-} MC Vertex_Z location; Z_{MC} cm; counts",200,-0.2,0.2); TH1F *hdsm_rmp = new TH1F("hdsm_rmp","Decay length of Ds^{-}; RM/P cm; counts",200,0,2); TH1F *hdsm_rmp_mct = new TH1F("hdsm_rmp_mct","Decay length of Ds^{-} (MCTruth matched); RM/P cm; counts",200,0,2); // *** the data reader object PndAnalysis* theAnalysis = new PndAnalysis(); if (nevts==0) nevts= theAnalysis->GetEntries(); // *** RhoCandLists for the analysis RhoCandList mctrk, pp, dsp, dsm, eta, neutrino_e, eplus, kam, kap, pimds, kambase, kapbase, pimdsbase; RhoCandList tm_dsm, vf_dsm, reco_dsm, mcf_pi0, mcf_eta; //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); // ds mass double pdgmass_ds = 1.96849; // *** Mass selector for the Ds-/Ds+ cands RhoMassParticleSelector *dsMassSel=new RhoMassParticleSelector("dsSelector",pdgmass_ds,0.5); // GeV RhoMassParticleSelector *ds_fineMassSel=new RhoMassParticleSelector("ds_fineSelector",pdgmass_ds,0.5); //int sum=0; // *** // the event loop // *** while (theAnalysis->GetEvent() && i++GetEntry(i); // *** the MC Truth objects //theAnalysis->FillList(mctrk,"McTruth"); // *** Select with no PID info ('All'); type and mass are set theAnalysis->FillList(kambase, "KaonLooseMinus"); theAnalysis->FillList(kapbase, "KaonLoosePlus"); theAnalysis->FillList(pimdsbase, "PionLooseMinus"); kam.Cleanup(); kap.Cleanup(); pimds.Cleanup(); // select pdg code for (j=0; jPdgCode() == -321)kam.Add(kambase[j]); } for (j=0; jPdgCode() == 321) kap.Add(kapbase[j]); } for (j=0; jGetMcTruth(); // RhoCandidate *mother = pimdsbase[j]->TheMother(); // int mopdg=-1; // if (mother) mopdg = mother->PdgCode();&& mopdg == -431 if (pimdsbase[j]->PdgCode() == -211 ) pimds.Add(pimdsbase[j]); } } // end event loop //cout << "sum = " << sum <cd(); // plot decay length of ds- timer.Stop(); Double_t rtime = timer.RealTime(); Double_t ctime = timer.CpuTime(); printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime); } // end of macro