- Compute parameter covariances including absorber dispersion effects
[u/mrichter/AliRoot.git] / MUON / MUONmassPlot_ESD.C
CommitLineData
5473f16a 1#if !defined(__CINT__) || defined(__MAKECINT__)
2// ROOT includes
48d6b312 3#include "TTree.h"
5473f16a 4#include "TBranch.h"
5#include "TClonesArray.h"
6#include "TLorentzVector.h"
7#include "TFile.h"
8#include "TH1.h"
ac3c5325 9#include "TH2.h"
5473f16a 10#include "TParticle.h"
11#include "TTree.h"
12#include <Riostream.h>
8cde4af5 13#include <TGeoManager.h>
21939432 14#include <TROOT.h>
5473f16a 15
16// STEER includes
17#include "AliRun.h"
0ff94351 18#include "AliLog.h"
5473f16a 19#include "AliRunLoader.h"
20#include "AliHeader.h"
21#include "AliLoader.h"
22#include "AliStack.h"
0ff94351 23#include "AliMagFMaps.h"
1c2bdf00 24#include "AliESDEvent.h"
995a61aa 25#include "AliESDVertex.h"
0ff94351 26#include "AliTracker.h"
5473f16a 27
28// MUON includes
211c52eb 29#include "AliMUONTrackParam.h"
37827b29 30#include "AliMUONTrackExtrap.h"
5473f16a 31#include "AliESDMuonTrack.h"
32#endif
33//
34// Macro MUONmassPlot.C for ESD
35// Ch. Finck, Subatech, April. 2004
36//
37
38// macro to make invariant mass plots
39// for combinations of 2 muons with opposite charges,
40// from root file "MUON.tracks.root" containing the result of track reconstruction.
41// Histograms are stored on the "MUONmassPlot.root" file.
42// introducing TLorentzVector for parameter calculations (Pt, P,rap,etc...)
43// using Invariant Mass for rapidity.
44
45// Arguments:
8cde4af5 46// ExtrapToVertex (default -1)
47// <0: no extrapolation;
48// =0: extrapolation to (0,0,0);
49// >0: extrapolation to ESDVertex if available, else to (0,0,0)
5473f16a 50// FirstEvent (default 0)
51// LastEvent (default 0)
52// ResType (default 553)
53// 553 for Upsilon, anything else for J/Psi
54// Chi2Cut (default 100)
55// to keep only tracks with chi2 per d.o.f. < Chi2Cut
56// PtCutMin (default 1)
57// to keep only tracks with transverse momentum > PtCutMin
58// PtCutMax (default 10000)
59// to keep only tracks with transverse momentum < PtCutMax
60// massMin (default 9.17 for Upsilon)
61// & massMax (default 9.77 for Upsilon)
62// to calculate the reconstruction efficiency for resonances with invariant mass
63// massMin < mass < massMax.
64
65// Add parameters and histograms for analysis
66
690d2205 67Bool_t MUONmassPlot(char* filename = "generated/galice.root", Int_t ExtrapToVertex = -1, char* geoFilename = "geometry.root",
8cde4af5 68 Int_t FirstEvent = 0, Int_t LastEvent = 10000, char* esdFileName = "AliESDs.root", Int_t ResType = 553,
5473f16a 69 Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,
70 Float_t massMin = 9.17,Float_t massMax = 9.77)
71{
72 cout << "MUONmassPlot " << endl;
73 cout << "FirstEvent " << FirstEvent << endl;
74 cout << "LastEvent " << LastEvent << endl;
75 cout << "ResType " << ResType << endl;
76 cout << "Chi2Cut " << Chi2Cut << endl;
77 cout << "PtCutMin " << PtCutMin << endl;
78 cout << "PtCutMax " << PtCutMax << endl;
79 cout << "massMin " << massMin << endl;
80 cout << "massMax " << massMax << endl;
81
82
83 //Reset ROOT and connect tree file
84 gROOT->Reset();
85
5473f16a 86 // File for histograms and histogram booking
87 TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
88 TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
ac3c5325 89 TH1F *hPtMuonPlus = new TH1F("hPtMuonPlus", "Muon+ Pt (GeV/c)", 100, 0., 20.);
90 TH1F *hPtMuonMinus = new TH1F("hPtMuonMinus", "Muon- Pt (GeV/c)", 100, 0., 20.);
5473f16a 91 TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
92 TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
93 TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
ac3c5325 94 TH1F *hInvMassBg = new TH1F("hInvMassBg", "Mu+Mu- invariant mass BG(GeV/c2)", 480, 0., 12.);
f57d136a 95 TH2F *hInvMassAll_vs_Pt = new TH2F("hInvMassAll_vs_Pt","hInvMassAll_vs_Pt",480,0.,12.,80,0.,20.);
96 TH2F *hInvMassBgk_vs_Pt = new TH2F("hInvMassBgk_vs_Pt","hInvMassBgk_vs_Pt",480,0.,12.,80,0.,20.);
97 TH1F *hInvMassRes;
6ab68e5f 98 TH1F *hPrimaryVertex = new TH1F("hPrimaryVertex","SPD reconstructed Z vertex",150,-15,15);
5473f16a 99
100 if (ResType == 553) {
101 hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
102 } else {
103 hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.);
104 }
105
106 TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5);
107 TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5,-2);
108 TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
109 TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
ac3c5325 110 TH2F *hThetaPhiPlus = new TH2F("hThetaPhiPlus", "Theta vs Phi +", 760, -190., 190., 400, 160., 180.);
111 TH2F *hThetaPhiMinus = new TH2F("hThetaPhiMinus", "Theta vs Phi -", 760, -190., 190., 400, 160., 180.);
5473f16a 112
113
114 // settings
115 Int_t EventInMass = 0;
6678cd54 116 Int_t EventInMassMatch = 0;
117 Int_t NbTrigger = 0;
118
5473f16a 119 Float_t muonMass = 0.105658389;
120// Float_t UpsilonMass = 9.46037;
121// Float_t JPsiMass = 3.097;
122
22ccc301 123 Int_t fCharge1, fCharge2;
5473f16a 124 Double_t fPxRec1, fPyRec1, fPzRec1, fE1;
125 Double_t fPxRec2, fPyRec2, fPzRec2, fE2;
5473f16a 126
127 Int_t ntrackhits, nevents;
128 Double_t fitfmin;
e3c6ae4c 129 Double_t fZVertex=0;
211c52eb 130 Double_t fYVertex=0;
131 Double_t fXVertex=0;
690d2205 132 Double_t errXVtx=0;
133 Double_t errYVtx=0;
5473f16a 134
135 TLorentzVector fV1, fV2, fVtot;
cb75342e 136
8cde4af5 137 // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
138 if (!gGeoManager) {
139 TGeoManager::Import(geoFilename);
140 if (!gGeoManager) {
141 Error("MUONmass_ESD", "getting geometry from file %s failed", filename);
142 return kFALSE;
143 }
144 }
145
0ff94351 146 // set mag field
147 // waiting for mag field in CDB
148 printf("Loading field map...\n");
b97b210c 149 AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
0ff94351 150 AliTracker::SetFieldMap(field, kFALSE);
cb75342e 151
5473f16a 152 // open run loader and load gAlice, kinematics and header
153 AliRunLoader* runLoader = AliRunLoader::Open(filename);
154 if (!runLoader) {
8cde4af5 155 Error("MUONmass_ESD", "getting run loader from file %s failed", filename);
5473f16a 156 return kFALSE;
157 }
6b092dfc 158/*
159 runLoader->LoadgAlice();
5473f16a 160 if (!gAlice) {
161 Error("MUONmass_ESD", "no galice object found");
162 return kFALSE;
163 }
6b092dfc 164*/
5473f16a 165
166 // open the ESD file
167 TFile* esdFile = TFile::Open(esdFileName);
168 if (!esdFile || !esdFile->IsOpen()) {
169 Error("MUONmass_ESD", "opening ESD file %s failed", esdFileName);
170 return kFALSE;
171 }
172
1c2bdf00 173 AliESDEvent* esd = new AliESDEvent();
48d6b312 174 TTree* tree = (TTree*) esdFile->Get("esdTree");
175 if (!tree) {
176 Error("CheckESD", "no ESD tree found");
177 return kFALSE;
178 }
8b4a9b89 179// tree->SetBranchAddress("ESD", &esd);
180 esd->ReadFromTree(tree);
f9ebb3bd 181
48d6b312 182
5473f16a 183 runLoader->LoadHeader();
184 nevents = runLoader->GetNumberOfEvents();
211c52eb 185
186 AliMUONTrackParam trackParam;
0ff94351 187
5473f16a 188 // Loop over events
189 for (Int_t iEvent = FirstEvent; iEvent <= TMath::Min(LastEvent, nevents - 1); iEvent++) {
190
191 // get current event
192 runLoader->GetEvent(iEvent);
48d6b312 193
5473f16a 194 // get the event summary data
48d6b312 195 tree->GetEvent(iEvent);
5473f16a 196 if (!esd) {
48d6b312 197 Error("CheckESD", "no ESD object found for event %d", iEvent);
5473f16a 198 return kFALSE;
199 }
200
f9ebb3bd 201 // get the SPD reconstructed vertex (vertexer) and fill the histogram
43939bd8 202 AliESDVertex* Vertex = (AliESDVertex*) esd->GetVertex();
43939bd8 203 if (Vertex->GetNContributors()) {
211c52eb 204 fZVertex = Vertex->GetZv();
205 fYVertex = Vertex->GetYv();
206 fXVertex = Vertex->GetXv();
690d2205 207 errXVtx = Vertex->GetXRes();
208 errYVtx = Vertex->GetYRes();
211c52eb 209 }
f9ebb3bd 210 hPrimaryVertex->Fill(fZVertex);
211
48d6b312 212 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
5473f16a 213
214 // printf("\n Nb of events analysed: %d\r",iEvent);
48d6b312 215 // cout << " number of tracks: " << nTracks <<endl;
5473f16a 216
37827b29 217 // set the magnetic field for track extrapolations
218 AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
5473f16a 219 // loop over all reconstructed tracks (also first track of combination)
220 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
221
8cde4af5 222 AliESDMuonTrack* muonTrack = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack)));
5473f16a 223
8cde4af5 224 // extrapolate to vertex if required and available
225 if (ExtrapToVertex > 0 && Vertex->GetNContributors()) {
22ccc301 226 trackParam.GetParamFromUncorrected(*muonTrack);
690d2205 227 AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex, errXVtx, errYVtx);
8cde4af5 228 trackParam.SetParamFor(*muonTrack); // put the new parameters in this copy of AliESDMuonTrack
229 } else if ((ExtrapToVertex > 0 && !Vertex->GetNContributors()) || ExtrapToVertex == 0){
22ccc301 230 trackParam.GetParamFromUncorrected(*muonTrack);
690d2205 231 AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0., 0., 0.);
8cde4af5 232 trackParam.SetParamFor(*muonTrack); // put the new parameters in this copy of AliESDMuonTrack
211c52eb 233 }
8cde4af5 234
22ccc301 235 fCharge1 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
236
237 muonTrack->LorentzP(fV1);
8cde4af5 238
5473f16a 239 ntrackhits = muonTrack->GetNHit();
240 fitfmin = muonTrack->GetChi2();
241
242 // transverse momentum
243 Float_t pt1 = fV1.Pt();
244
245 // total momentum
246 Float_t p1 = fV1.P();
247
248 // Rapidity
249 Float_t rapMuon1 = fV1.Rapidity();
250
251 // chi2 per d.o.f.
252 Float_t ch1 = fitfmin / (2.0 * ntrackhits - 5);
253// printf(" px %f py %f pz %f NHits %d Norm.chi2 %f charge %d\n",
22ccc301 254// fPxRec1, fPyRec1, fPzRec1, ntrackhits, ch1, fCharge1);
5473f16a 255
256 // condition for good track (Chi2Cut and PtCut)
257
258 if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
259
260 // fill histos hPtMuon and hChi2PerDof
261 hPtMuon->Fill(pt1);
262 hPMuon->Fill(p1);
263 hChi2PerDof->Fill(ch1);
264 hRapMuon->Fill(rapMuon1);
22ccc301 265 if (fCharge1 > 0) {
ac3c5325 266 hPtMuonPlus->Fill(pt1);
22ccc301 267 hThetaPhiPlus->Fill(fV1.Phi()*180./TMath::Pi(),fV1.Theta()*180./TMath::Pi());
ac3c5325 268 } else {
269 hPtMuonMinus->Fill(pt1);
22ccc301 270 hThetaPhiMinus->Fill(fV1.Phi()*180./TMath::Pi(),fV1.Theta()*180./TMath::Pi());
ac3c5325 271 }
5473f16a 272 // loop over second track of combination
273 for (Int_t iTrack2 = iTrack + 1; iTrack2 < nTracks; iTrack2++) {
274
8cde4af5 275 AliESDMuonTrack* muonTrack2 = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack2)));
276
277 // extrapolate to vertex if required and available
278 if (ExtrapToVertex > 0 && Vertex->GetNContributors()) {
22ccc301 279 trackParam.GetParamFromUncorrected(*muonTrack2);
690d2205 280 AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex, errXVtx, errYVtx);
8cde4af5 281 trackParam.SetParamFor(*muonTrack2); // put the new parameters in this copy of AliESDMuonTrack
282 } else if ((ExtrapToVertex > 0 && !Vertex->GetNContributors()) || ExtrapToVertex == 0){
22ccc301 283 trackParam.GetParamFromUncorrected(*muonTrack2);
690d2205 284 AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0., 0., 0.);
8cde4af5 285 trackParam.SetParamFor(*muonTrack2); // put the new parameters in this copy of AliESDMuonTrack
211c52eb 286 }
8cde4af5 287
8cde4af5 288 fCharge2 = Int_t(TMath::Sign(1.,muonTrack2->GetInverseBendingMomentum()));
5473f16a 289
22ccc301 290 muonTrack2->LorentzP(fV2);
5473f16a 291
8cde4af5 292 ntrackhits = muonTrack2->GetNHit();
293 fitfmin = muonTrack2->GetChi2();
5473f16a 294
295 // transverse momentum
296 Float_t pt2 = fV2.Pt();
297
298 // chi2 per d.o.f.
299 Float_t ch2 = fitfmin / (2.0 * ntrackhits - 5);
300
301 // condition for good track (Chi2Cut and PtCut)
302 if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax)) {
303
304 // condition for opposite charges
22ccc301 305 if ((fCharge1 * fCharge2) == -1) {
5473f16a 306
307 // invariant mass
308 fVtot = fV1 + fV2;
309 Float_t invMass = fVtot.M();
310
311 // fill histos hInvMassAll and hInvMassRes
312 hInvMassAll->Fill(invMass);
313 hInvMassRes->Fill(invMass);
ac3c5325 314 hInvMassAll_vs_Pt->Fill(invMass,fVtot.Pt());
6678cd54 315 Int_t ptTrig;
316 if (ResType == 553)
8d4fefab 317 ptTrig = 0x20;// mask for Hpt unlike sign pair
6678cd54 318 else
8d4fefab 319 ptTrig = 0x10;// mask for Lpt unlike sign pair
6678cd54 320
f57d136a 321 if (esd->GetTriggerMask() & ptTrig) NbTrigger++;
5473f16a 322 if (invMass > massMin && invMass < massMax) {
323 EventInMass++;
8cde4af5 324 if (muonTrack2->GetMatchTrigger() && (esd->GetTriggerMask() & ptTrig))// match with trigger
6678cd54 325 EventInMassMatch++;
326
5473f16a 327 hRapResonance->Fill(fVtot.Rapidity());
328 hPtResonance->Fill(fVtot.Pt());
329 }
330
22ccc301 331 } // if (fCharge1 * fCharge2) == -1)
5473f16a 332 } // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
8cde4af5 333 delete muonTrack2;
5473f16a 334 } // for (Int_t iTrack2 = iTrack + 1; iTrack2 < iTrack; iTrack2++)
335 } // if (ch1 < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
8cde4af5 336 delete muonTrack;
5473f16a 337 } // for (Int_t iTrack = 0; iTrack < nrectracks; iTrack++)
338
339 hNumberOfTrack->Fill(nTracks);
48d6b312 340 // esdFile->Delete();
5473f16a 341 } // for (Int_t iEvent = FirstEvent;
342
ac3c5325 343// Loop over events for bg event
344
345 Double_t thetaPlus, phiPlus;
346 Double_t thetaMinus, phiMinus;
347 Float_t PtMinus, PtPlus;
348
349 for (Int_t iEvent = 0; iEvent < hInvMassAll->Integral(); iEvent++) {
350
351 hThetaPhiPlus->GetRandom2(phiPlus, thetaPlus);
352 hThetaPhiMinus->GetRandom2(phiMinus,thetaMinus);
353 PtPlus = hPtMuonPlus->GetRandom();
354 PtMinus = hPtMuonMinus->GetRandom();
355
356 fPxRec1 = PtPlus * TMath::Cos(TMath::Pi()/180.*phiPlus);
357 fPyRec1 = PtPlus * TMath::Sin(TMath::Pi()/180.*phiPlus);
358 fPzRec1 = PtPlus / TMath::Tan(TMath::Pi()/180.*thetaPlus);
359
360 fE1 = TMath::Sqrt(muonMass * muonMass + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1);
361 fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1);
362
363 fPxRec2 = PtMinus * TMath::Cos(TMath::Pi()/180.*phiMinus);
364 fPyRec2 = PtMinus * TMath::Sin(TMath::Pi()/180.*phiMinus);
365 fPzRec2 = PtMinus / TMath::Tan(TMath::Pi()/180.*thetaMinus);
366
367 fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2);
368 fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2);
369
370 // invariant mass
371 fVtot = fV1 + fV2;
372
373 // fill histos hInvMassAll and hInvMassRes
374 hInvMassBg->Fill(fVtot.M());
375 hInvMassBgk_vs_Pt->Fill( fVtot.M(), fVtot.Pt() );
376 }
377
5473f16a 378 histoFile->Write();
379 histoFile->Close();
380
6678cd54 381 cout << endl;
5473f16a 382 cout << "EventInMass " << EventInMass << endl;
6678cd54 383 cout << "NbTrigger " << NbTrigger << endl;
384 cout << "EventInMass match with trigger " << EventInMassMatch << endl;
5473f16a 385
386 return kTRUE;
387}
388