// LastEvent (default 0)
// ResType (default 553)
// 553 for Upsilon, anything else for J/Psi
-// NSigma (default 3)
-// the number of combinations is counted around the resonance mass
-// within +/- NSigma times the nominal sigma's
-// (0.099 GeV for Upsilon, 0.0615 GeV for J/Psi)
// Chi2Cut (default 100)
// to keep only tracks with chi2 per d.o.f. < Chi2Cut
-// PtCut (default 1)
-// to keep only tracks with transverse momentum > PtCut
-
-// IMPORTANT NOTICE FOR USERS:
-// under "root" or "root.exe", execute the following commands:
-// 1. "gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER")" to get the right path at compilation time
-// 2. ".x loadlibs.C" to load the shared libraries
-// 3. ".L MUONrecoNtuple.C+"
-// 4. ".x MUONmassPlot.C()" with the right arguments according to the list above
-
-void MUONmassPlot(Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553, Float_t Nsig = 3., Float_t Chi2Cut = 100., Float_t PtCut = 1.)
+// 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 << "PtCut" << PtCut << endl;
+ 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
//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");
UInt_t fUniqueID;
UInt_t fBits;
Int_t Tracks_;
- Int_t Tracks_fCharge[5];
- Float_t Tracks_fPxRec[5];
- Float_t Tracks_fPyRec[5];
- Float_t Tracks_fPzRec[5];
- Float_t Tracks_fZRec[5];
- Float_t Tracks_fZRec1[5];
- Int_t Tracks_fNHits[5];
- Float_t Tracks_fChi2[5];
- Float_t Tracks_fPxGen[5];
- Float_t Tracks_fPyGen[5];
- Float_t Tracks_fPzGen[5];
- UInt_t Tracks_fUniqueID[5];
- UInt_t Tracks_fBits[5];
+ 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);
// 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)", 240, 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, 1., 5.);
+ 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,-4.5.,-2);
+ TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
+ 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 <<endl;
+ cout << " number of tracks: " << Tracks_ <<endl;
nbytes += MUONtrackReco->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 > 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
// 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 > 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]);
+ 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 > PtCut))
+ } // if ((Tracks_fChi2[t2] < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
} // for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++)
- } // if ((Tracks_fChi2[t1] < Chi2Cut) && (pt1 > PtCut))
+ } // 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 MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t muonMass)
{
- Float_t muonMass = 0.10566;
+ // 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;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff