changes in the MagF constructor
[u/mrichter/AliRoot.git] / MUON / DIMUONFakes.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 // ROOT includes
3 #include <TTree.h>
4 #include <TFile.h>
5 #include <TH1.h>
6 #include <TCanvas.h>
7 #include <Riostream.h>
8 #include <TROOT.h>
9 #include <TClonesArray.h>
10 #include <TGeoGlobalMagField.h>
11 #include <TLorentzVector.h>
12
13 // STEER includes
14 #include "AliLog.h"
15 #include "AliMagF.h"
16 #include "AliESDEvent.h"
17 #include "AliESDMuonTrack.h"
18 #include "AliESDMuonCluster.h"
19 #include "AliCDBPath.h"
20 #include "AliCDBEntry.h"
21 #include "AliCDBManager.h"
22
23 // MUON includes
24 #include "AliMUONTrack.h"
25 #include "AliMUONVTrackStore.h"
26 #include "AliMUONTrackParam.h"
27 #include "AliMUONTrackExtrap.h"
28 #include "AliMUONESDInterface.h"
29 #include "AliMUONRecoCheck.h"
30 #include "AliMUONVCluster.h"
31 #include "AliMUONRecoParam.h"
32 #endif
33
34 /// \ingroup macros
35 /// \file DIMUONFakes.C
36 ///
37 /// \author Ph. Pillot, Subatech, March. 2009
38 ///
39 /// Macro to study the effects of fake tracks on the dimuon spectra
40 /// Results are saved in the root file DiFakes.root
41 /// Results are relevent provided that you use the same recoParams as for the reconstruction
42
43 void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut);
44 TTree* GetESDTree(TFile *esdFile);
45 Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut);
46 AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
47                                 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);
48
49 //-----------------------------------------------------------------------
50 void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
51                  const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/")
52 {
53   
54   //Reset ROOT and connect tree file
55   gROOT->Reset();
56   
57   // File for histograms and histogram booking
58   TFile *histoFile = new TFile("DiFakes.root", "RECREATE");
59   
60   TH1F *hMass = new TH1F("hMass", "Dimuon mass distribution (GeV/c^{2})", 100, 0., 12.);
61   TH1F *hMassM = new TH1F("hMassM", "matched track mass distribution (GeV/c^{2})", 100, 0., 12.);
62   TH1F *hMassF = new TH1F("hMassF", "fake track mass distribution (GeV/c^{2})", 100, 0., 12.);
63   TH1F *hP = new TH1F("hP", "Dimuon P distribution (GeV/c)", 100, 0., 200.);
64   TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
65   TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
66   TH1F *hPt = new TH1F("hPt", "Dimuon Pt distribution (GeV/c)", 100, 0., 20.);
67   TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
68   TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
69   TH1F *hY = new TH1F("hY"," Dimuon rapidity distribution",100,-10,0);
70   TH1F *hYM = new TH1F("hYM"," matched track rapidity distribution",100,-10,0);
71   TH1F *hYF = new TH1F("hYF"," fake track rapidity distribution",100,-10,0);
72   TH1F *hEta = new TH1F("hEta"," Dimuon pseudo-rapidity distribution",100,-10,0);
73   TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
74   TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
75   TH1F *hPhi = new TH1F("hPhi"," Dimuon phi distribution",100,-1.,9.);
76   TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
77   TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
78   
79   // prepare for analysis
80   AliMUONRecoParam *recoParam = 0x0;
81   Double_t sigmaCut = -1;
82   Prepare(recoParam, sigmaCut);
83   
84   // link to reconstructed tracks
85   TFile* esdFile = TFile::Open(esdFileName);
86   TTree* esdTree = GetESDTree(esdFile);
87   AliESDEvent* esd = new AliESDEvent();
88   esd->ReadFromTree(esdTree);
89   
90   // link to simulated tracks
91   AliMUONRecoCheck rc(esdFileName, SimDir);
92   
93   TLorentzVector vMu1, vMu2, vDiMu;
94   
95   // Loop over ESD events
96   FirstEvent = TMath::Max(0, FirstEvent);
97   LastEvent = (LastEvent>=0) ? TMath::Min((Int_t)esdTree->GetEntries() - 1, LastEvent) : (Int_t)esdTree->GetEntries() - 1;
98   for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
99     
100     // get the ESD of current event
101     esdTree->GetEvent(iEvent);
102     if (!esd) {
103       Error("CheckESD", "no ESD object found for event %d", iEvent);
104       return;
105     }
106     
107     // convert TrackRef to MUON tracks
108     AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
109     
110     // loop over ESD tracks and flag them
111     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
112     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
113       
114       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
115       
116       // skip ghosts
117       if (!muonTrack->ContainTrackerData()) continue;
118       
119       // try to match the reconstructed track with a simulated one
120       Float_t fractionOfMatchCluster = 0.;
121       AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*muonTrack, *trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut);
122       
123       // take actions according to matching result
124       if (matchedTrackRef) {
125         
126         // flag matched tracks
127         muonTrack->SetLabel(matchedTrackRef->GetUniqueID());
128         
129         // remove already matched trackRefs
130         trackRefStore->Remove(*matchedTrackRef);
131         
132       } else {
133         
134         // flag fake tracks
135         muonTrack->SetLabel(-1);
136         
137       }
138       
139     }
140     
141     // double loop over ESD tracks, build pairs and fill histograms according to their label
142     for (Int_t iTrack1 = 0; iTrack1 <  nTracks;  iTrack1++) {
143       AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
144       
145       // skip ghosts
146       if (!muonTrack1->ContainTrackerData()) continue;
147       
148       // get track info
149       Short_t charge1 = muonTrack1->Charge();
150       Int_t label1 = muonTrack1->GetLabel();
151       muonTrack1->LorentzP(vMu1);
152       
153       for (Int_t iTrack2 = iTrack1+1; iTrack2 <  nTracks;  iTrack2++) {
154         AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
155         
156         // skip ghosts
157         if (!muonTrack2->ContainTrackerData()) continue;
158         
159         // keep only opposite sign pairs
160         Short_t charge2 = muonTrack2->Charge();
161         if (charge1*charge2 > 0) continue;
162         
163         // get track info
164         Int_t label2 = muonTrack2->GetLabel();
165         muonTrack2->LorentzP(vMu2);
166         
167         // compute kinematics of the pair
168         vDiMu = vMu1 + vMu2;
169         Float_t mass = vDiMu.M();
170         Float_t p = vDiMu.P();
171         Float_t pt = vDiMu.Pt();
172         Float_t y = vDiMu.Rapidity();
173         Float_t eta = vDiMu.Eta();
174         Float_t phi = vDiMu.Phi();
175         if (phi < 0) phi += 2.*TMath::Pi();
176         
177         // fill global histograms
178         hMass->Fill(mass);
179         hP->Fill(p);
180         hPt->Fill(pt);
181         hY->Fill(y);
182         hEta->Fill(eta);
183         hPhi->Fill(phi);
184         
185         // fill histograms according to labels
186         if (label1 >= 0 && label2 >= 0) {
187           
188           hMassM->Fill(mass);
189           hPM->Fill(p);
190           hPtM->Fill(pt);
191           hYM->Fill(y);
192           hEtaM->Fill(eta);
193           hPhiM->Fill(phi);
194           
195         } else {
196           
197           hMassF->Fill(mass);
198           hPF->Fill(p);
199           hPtF->Fill(pt);
200           hYF->Fill(y);
201           hEtaF->Fill(eta);
202           hPhiF->Fill(phi);
203           
204         }
205         
206       }
207       
208     }
209       
210   } // end of loop over events
211   
212   // plot results
213   TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600);
214   cDiFakesSummary.Divide(3,2);
215   cDiFakesSummary.cd(1);
216   cDiFakesSummary.GetPad(1)->SetLogy();
217   hMass->Draw();
218   hMass->SetMinimum(0.5);
219   hMassM->Draw("same");
220   hMassM->SetLineColor(4);
221   hMassF->Draw("same");
222   hMassF->SetLineColor(2);
223   hMassF->SetFillColor(2);
224   hMassF->SetFillStyle(3017);
225   cDiFakesSummary.cd(2);
226   cDiFakesSummary.GetPad(3)->SetLogy();
227   hP->Draw();
228   hP->SetMinimum(0.5);
229   hPM->Draw("same");
230   hPM->SetLineColor(4);
231   hPF->Draw("same");
232   hPF->SetLineColor(2);
233   hPF->SetFillColor(2);
234   hPF->SetFillStyle(3017);
235   cDiFakesSummary.cd(3);
236   cDiFakesSummary.GetPad(4)->SetLogy();
237   hPt->Draw();
238   hPt->SetMinimum(0.5);
239   hPtM->Draw("same");
240   hPtM->SetLineColor(4);
241   hPtF->Draw("same");
242   hPtF->SetLineColor(2);
243   hPtF->SetFillColor(2);
244   hPtF->SetFillStyle(3017);
245   cDiFakesSummary.cd(4);
246   cDiFakesSummary.GetPad(2)->SetLogy();
247   hY->Draw();
248   hY->SetMinimum(0.5);
249   hYM->Draw("same");
250   hYM->SetLineColor(4);
251   hYF->Draw("same");
252   hYF->SetLineColor(2);
253   hYF->SetFillColor(2);
254   hYF->SetFillStyle(3017);
255   cDiFakesSummary.cd(5);
256   cDiFakesSummary.GetPad(5)->SetLogy();
257   hEta->Draw();
258   hEta->SetMinimum(0.5);
259   hEtaM->Draw("same");
260   hEtaM->SetLineColor(4);
261   hEtaF->Draw("same");
262   hEtaF->SetLineColor(2);
263   hEtaF->SetFillColor(2);
264   hEtaF->SetFillStyle(3017);
265   cDiFakesSummary.cd(6);
266   cDiFakesSummary.GetPad(6)->SetLogy();
267   hPhi->Draw();
268   hPhi->SetMinimum(0.5);
269   hPhiM->Draw("same");
270   hPhiM->SetLineColor(4);
271   hPhiF->Draw("same");
272   hPhiF->SetLineColor(2);
273   hPhiF->SetFillColor(2);
274   hPhiF->SetFillStyle(3017);
275   
276   // save results
277   histoFile->cd();
278   histoFile->Write();
279   cDiFakesSummary.Write();
280   histoFile->Close();
281   
282   // clear memory
283   delete esd;
284   esdFile->Close();
285   delete recoParam;
286   
287 }
288
289 //-----------------------------------------------------------------------
290 void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut)
291 {
292   /// Set the magnetic field and return recoParam and sigmaCut to associate cluster and trackRef
293   
294   // prepare OCDB access
295   AliCDBManager* man = AliCDBManager::Instance();
296   man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
297   man->SetRun(0);
298   
299   // set  mag field 
300   // waiting for mag field in CDB 
301   if (!TGeoGlobalMagField::Instance()->GetField()) {
302     printf("Loading field map...\n");
303     AliMagF* field = new AliMagF("Maps","Maps",1.,1.,AliMagF::k5kG);
304     TGeoGlobalMagField::Instance()->SetField(field);
305   }
306   // set the magnetic field for track extrapolations
307   AliMUONTrackExtrap::SetField();
308   
309   // Load initial reconstruction parameters from OCDB
310   AliCDBPath path("MUON","Calib","RecoParam");
311   AliCDBEntry *entry=man->Get(path.GetPath());
312   if(entry) {
313     recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
314     entry->SetOwner(0);
315     AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
316   }
317   if (!recoParam) {
318     printf("Couldn't find RecoParam object in OCDB: create default one");
319     recoParam = AliMUONRecoParam::GetLowFluxParam();
320   }
321   
322   Info("MUONFakes", "\n recontruction parameters:");
323   recoParam->Print("FULL");
324   AliMUONESDInterface::ResetTracker(recoParam);
325   
326   // sigma cut to associate clusters with TrackRefs in case the label are not used
327   sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking(); 
328   
329 }
330
331 //-----------------------------------------------------------------------
332 TTree* GetESDTree(TFile *esdFile)
333 {
334   /// Check that the file is properly open
335   /// Return pointer to the ESD Tree
336   
337   if (!esdFile || !esdFile->IsOpen()) {
338     Error("GetESDTree", "opening ESD file failed");
339     exit(-1);
340   }
341   
342   TTree* tree = (TTree*) esdFile->Get("esdTree");
343   if (!tree) {
344     Error("GetESDTree", "no ESD tree found");
345     exit(-1);
346   }
347   
348   return tree;
349   
350 }
351
352 //-----------------------------------------------------------------------
353 Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut)
354 {
355   /// Try to match 2 tracks
356   
357   Bool_t compTrack[10];
358   Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack);
359   fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters());
360   
361   if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2
362       (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5
363       fractionOfMatchCluster > 0.5) return kTRUE;                       // more than 50% of clusters matched
364   else return kFALSE;
365   
366 }
367
368 //-----------------------------------------------------------------------
369 AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
370                                 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut)
371 {
372   /// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters
373   
374   AliMUONTrack *matchedTrackRef = 0x0;
375   fractionOfMatchCluster = 0.;
376   
377   if (useLabel) { // by using the MC label
378     
379     // get the corresponding simulated track if any
380     Int_t label = muonTrack.GetLabel();
381     matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label);
382     
383     // get the fraction of matched clusters
384     if (matchedTrackRef) {
385       Int_t nMatchClusters = 0;
386       if (muonTrack.ClustersStored()) {
387         AliESDMuonCluster* cluster = (AliESDMuonCluster*) muonTrack.GetClusters().First();
388         while (cluster) {
389           if (cluster->GetLabel() == label) nMatchClusters++;
390           cluster = (AliESDMuonCluster*) muonTrack.GetClusters().After(cluster);
391         }
392       }
393       fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)muonTrack.GetNClusters());
394     }
395     
396   } else { // by comparing cluster/TrackRef positions
397     
398     // convert ESD track to MUON track
399     AliMUONTrack track;
400     AliMUONESDInterface::ESDToMUON(muonTrack,track);
401     
402     // look for the corresponding simulated track if any
403     TIter next(trackRefStore.CreateIterator());
404     AliMUONTrack* trackRef;
405     while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
406       
407       // check compatibility
408       Float_t f = 0.;
409       if (TrackMatched(track, *trackRef, f, sigmaCut)) {
410         matchedTrackRef = trackRef;
411         fractionOfMatchCluster = f;
412         break;
413       }
414       
415     }
416     
417   }
418   
419   return matchedTrackRef;
420   
421 }
422