// macro to make invariant mass plots // for combinations of 2 muons with opposite charges, // from root file "MUONtrackReco.root" containing the result of track reconstruction, // generated by the macro "MUONrecoNtuple.C". // Histograms are stored on the "MUONmassPlot.root" file. // A model for macros using the Ntuple in the file "MUONtrackReco.root" // may be found in "MUONtrackRecoModel.C": // it has been obtained by reading with Root a file "MUONtrackReco.root", // and executing the command: // MUONtrackReco->MakeCode("MUONtrackRecoModel.C") // Arguments: // FirstEvent (default 0) // LastEvent (default 0) // ResType (default 553) // 553 for Upsilon, anything else for J/Psi // Chi2Cut (default 100) // to keep only tracks with chi2 per d.o.f. < Chi2Cut // PtCutMin (default 1) // to keep only tracks with transverse momentum > PtCutMin // PtCutMax (default 10000) // to keep only tracks with transverse momentum < PtCutMax // massMin (default 9.17 for Upsilon) // & massMax (default 9.77 for Upsilon) // to calculate the reconstruction efficiency for resonances with invariant mass // massMin < mass < massMax. // compile MUONrecoNtuple.C and load shared library in this macro // Add parameters and histograms for analysis void MUONmassPlot(Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553, Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,Float_t massMin = 9.17,Float_t massMax = 9.77) { cout << "MUONmassPlot " << endl; cout << "FirstEvent " << FirstEvent << endl; cout << "LastEvent " << LastEvent << endl; cout << "ResType " << ResType << endl; // cout << "Nsig " << Nsig << endl; cout << "Chi2Cut " << Chi2Cut << endl; cout << "PtCutMin " << PtCutMin << endl; cout << "PtCutMax " << PtCutMax << endl; cout << "massMin " << massMin << endl; cout << "massMax " << massMax << endl; ////////////////////////////////////////////////////////// // This file has been automatically generated // (Thu Sep 21 14:53:11 2000 by ROOT version2.25/02) // from TTree MUONtrackReco/MUONtrackReco // found on file: MUONtrackReco.root ////////////////////////////////////////////////////////// //Reset ROOT and connect tree file gROOT->Reset(); // Dynamically link some shared libs if (gClassTable->GetID("AliRun") < 0) { gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER"); gROOT->LoadMacro("loadlibs.C"); loadlibs(); // compile MUONrecoNtuple and load shared library gSystem->CompileMacro("$(ALICE_ROOT)/macros/MUONrecoNtuple.C","kf"); } TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("MUONtrackReco.root"); if (!f) { f = new TFile("MUONtrackReco.root"); } TTree *MUONtrackReco = (TTree*)gDirectory->Get("MUONtrackReco"); MUONtrackReco->SetMakeClass(1); //Declaration of leaves types Int_t fEvent; UInt_t fUniqueID; UInt_t fBits; Int_t Tracks_; Int_t Tracks_fCharge[500]; Float_t Tracks_fPxRec[500]; Float_t Tracks_fPyRec[500]; Float_t Tracks_fPzRec[500]; Float_t Tracks_fZRec[500]; Float_t Tracks_fZRec1[500]; Int_t Tracks_fNHits[500]; Float_t Tracks_fChi2[500]; Float_t Tracks_fPxGen[500]; Float_t Tracks_fPyGen[500]; Float_t Tracks_fPzGen[500]; UInt_t Tracks_fUniqueID[500]; UInt_t Tracks_fBits[500]; //Set branch addresses //MUONtrackReco->SetBranchAddress("Header",&Header); MUONtrackReco->SetBranchAddress("fEvent",&fEvent); MUONtrackReco->SetBranchAddress("fUniqueID",&fUniqueID); MUONtrackReco->SetBranchAddress("fBits",&fBits); // MUONtrackReco->SetBranchAddress("Tracks_",&Tracks_); MUONtrackReco->SetBranchAddress("Tracks",&Tracks_); MUONtrackReco->SetBranchAddress("Tracks.fCharge",Tracks_fCharge); MUONtrackReco->SetBranchAddress("Tracks.fPxRec",Tracks_fPxRec); MUONtrackReco->SetBranchAddress("Tracks.fPyRec",Tracks_fPyRec); MUONtrackReco->SetBranchAddress("Tracks.fPzRec",Tracks_fPzRec); MUONtrackReco->SetBranchAddress("Tracks.fZRec",Tracks_fZRec); MUONtrackReco->SetBranchAddress("Tracks.fZRec1",Tracks_fZRec1); MUONtrackReco->SetBranchAddress("Tracks.fNHits",Tracks_fNHits); MUONtrackReco->SetBranchAddress("Tracks.fChi2",Tracks_fChi2); MUONtrackReco->SetBranchAddress("Tracks.fPxGen",Tracks_fPxGen); MUONtrackReco->SetBranchAddress("Tracks.fPyGen",Tracks_fPyGen); MUONtrackReco->SetBranchAddress("Tracks.fPzGen",Tracks_fPzGen); MUONtrackReco->SetBranchAddress("Tracks.fUniqueID",Tracks_fUniqueID); MUONtrackReco->SetBranchAddress("Tracks.fBits",Tracks_fBits); // This is the loop skeleton // To read only selected branches, Insert statements like: // MUONtrackReco->SetBranchStatus("*",0); // disable all branches // TTreePlayer->SetBranchStatus("branchname",1); // activate branchname Int_t nentries = MUONtrackReco->GetEntries(); Int_t nbytes = 0; // for (Int_t i=0; iGetEntry(i); // } ///////////////////////////////////////////////////////////////// // Here comes the specialized part for MUONmassPlot ///////////////////////////////////////////////////////////////// // File for histograms and histogram booking TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE"); TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.); TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.); TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.); TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.); if (ResType == 553) { TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.); } else { TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.); } TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5); TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,2.,4.5); TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,2.,4.5); TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.); Int_t EventInMass = 0; Float_t muonMass = 0.105658389; Float_t UpsilonMass = 9.46037; Float_t JPsiMass = 3.097; // Loop over events for (Int_t event = FirstEvent; event <= TMath::Min(LastEvent, nentries - 1); event++) { // get current event Int_t NumberOfTrack = 0; cout <<" event: " << event <GetEntry(event); // loop over all reconstructed tracks (also first track of combination) for (Int_t t1 = 0; t1 < Tracks_; t1++) { // transverse momentum Float_t pt1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1]); // total momentum Float_t p1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1] + Tracks_fPzRec[t1] * Tracks_fPzRec[t1] ); // Rapidity Float_t rapMuon1 = rapParticle(Tracks_fPxRec[t1],Tracks_fPyRec[t1],Tracks_fPzRec[t1],muonMass); // chi2 per d.o.f. Float_t ch1 = Tracks_fChi2[t1] / (2.0 * Tracks_fNHits[t1] - 5); printf(" px %f py %f pz %f NHits %d Norm.chi2 %f \n",Tracks_fPxRec[t1],Tracks_fPyRec[t1],Tracks_fPzRec[t1],Tracks_fNHits[t1],ch1); // condition for good track (Chi2Cut and PtCut) if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) { // fill histos hPtMuon and hChi2PerDof NumberOfTrack++; hPtMuon->Fill(pt1); hPMuon->Fill(p1); hChi2PerDof->Fill(ch1); hRapMuon->Fill(rapMuon1); // loop over second track of combination for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++) { // transverse momentum Float_t pt2 = TMath::Sqrt(Tracks_fPxRec[t2] * Tracks_fPxRec[t2] + Tracks_fPyRec[t2] * Tracks_fPyRec[t2]); // chi2 per d.o.f. Float_t ch2 = Tracks_fChi2[t2] / (2.0 * Tracks_fNHits[t2] - 5); // condition for good track (Chi2Cut and PtCut) if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax)) { // condition for opposite charges if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1) { // invariant mass Float_t invMass = MuPlusMuMinusMass(Tracks_fPxRec[t1], Tracks_fPyRec[t1], Tracks_fPzRec[t1], Tracks_fPxRec[t2], Tracks_fPyRec[t2], Tracks_fPzRec[t2],muonMass); // fill histos hInvMassAll and hInvMassRes hInvMassAll->Fill(invMass); hInvMassRes->Fill(invMass); Float_t pxRes = Tracks_fPxRec[t1]+Tracks_fPxRec[t2]; Float_t pyRes = Tracks_fPyRec[t1]+Tracks_fPyRec[t2]; Float_t pzRes = Tracks_fPzRec[t1]+Tracks_fPzRec[t2]; if (invMass > massMin && invMass < massMax) { EventInMass++; Float_t rapRes = 0; if (ResType == 553) { rapRes = rapParticle(pxRes,pyRes,pzRes,UpsilonMass); } else { rapRes = rapParticle(pxRes,pyRes,pzRes,JPsiMass); } hRapResonance->Fill(rapRes); hPtResonance->Fill(TMath::Sqrt(pxRes*pxRes+pyRes*pyRes)); } } //if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1) } // if ((Tracks_fChi2[t2] < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax)) } // for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++) } // if ((Tracks_fChi2[t1] < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) ) } // for (Int_t t1 = 0; t1 < Tracks_; t1++) hNumberOfTrack->Fill(NumberOfTrack); } // for (Int_t event = FirstEvent; histoFile->Write(); histoFile->Close(); cout << "MUONmassPlot " << endl; cout << "FirstEvent " << FirstEvent << endl; cout << "LastEvent " << LastEvent << endl; cout << "ResType " << ResType << endl; // cout << "Nsig " << Nsig << endl; cout << "Chi2Cut " << Chi2Cut << endl; cout << "PtCutMin " << PtCutMin << endl; cout << "PtCutMax " << PtCutMax << endl; cout << "massMin " << massMin << endl; cout << "massMax " << massMax << endl; cout << "EventInMass " << EventInMass << endl; } Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t muonMass) { // return invariant mass for particle 1 & 2 whose masses are equal to muonMass Float_t e1 = TMath::Sqrt(muonMass * muonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1); Float_t e2 = TMath::Sqrt(muonMass * muonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2); return (TMath::Sqrt(2.0 * (muonMass * muonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2))); } Float_t rapParticle(Float_t Px, Float_t Py, Float_t Pz, Float_t Mass) { // return rapidity for particle // Particle energy Float_t Energy = TMath::Sqrt(Px*Px+Py*Py+Pz*Pz+Mass*Mass); // Rapidity Float_t Rapidity = 0.5*TMath::Log((Energy+Pz)/(Energy-Pz)); return Rapidity; }