95b7ace37e86c522c85fa0ff63fa26511ab8fff7
[u/mrichter/AliRoot.git] / MUON / DIMUONFakes.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 // ROOT includes
3 #include <TFile.h>
4 #include <TH1.h>
5 #include <TCanvas.h>
6 #include <Riostream.h>
7 #include <TROOT.h>
8 #include <TClonesArray.h>
9 #include <TLorentzVector.h>
10
11 // STEER includes
12 #include "AliLog.h"
13 #include "AliESDEvent.h"
14 #include "AliESDMuonTrack.h"
15 #include "AliCDBManager.h"
16
17 // MUON includes
18 #include "AliMUONCDB.h"
19 #include "AliMUONTrack.h"
20 #include "AliMUONVTrackStore.h"
21 #include "AliMUONTrackParam.h"
22 #include "AliMUONRecoCheck.h"
23 #include "AliMUONVCluster.h"
24 #include "AliMUONRecoParam.h"
25 #endif
26
27 /// \ingroup macros
28 /// \file DIMUONFakes.C
29 ///
30 /// \author Ph. Pillot, Subatech, March. 2009
31 ///
32 /// Macro to study the effects of fake tracks on the dimuon spectra
33 /// Results are saved in the root file DiFakes.root
34 /// Results are relevent provided that you use the same recoParams as for the reconstruction
35
36 Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut);
37 AliMUONTrack* MatchWithTrackRef(AliMUONTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
38                                 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);
39
40 //-----------------------------------------------------------------------
41 void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
42                  const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
43                  const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
44 {
45   
46   //Reset ROOT and connect tree file
47   gROOT->Reset();
48   
49   // File for histograms and histogram booking
50   TFile *histoFile = new TFile("DiFakes.root", "RECREATE");
51   
52   TH1F *hMass = new TH1F("hMass", "Dimuon mass distribution (GeV/c^{2})", 100, 0., 12.);
53   TH1F *hMassM = new TH1F("hMassM", "matched track mass distribution (GeV/c^{2})", 100, 0., 12.);
54   TH1F *hMassF = new TH1F("hMassF", "fake track mass distribution (GeV/c^{2})", 100, 0., 12.);
55   TH1F *hP = new TH1F("hP", "Dimuon P distribution (GeV/c)", 100, 0., 200.);
56   TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
57   TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
58   TH1F *hPt = new TH1F("hPt", "Dimuon Pt distribution (GeV/c)", 100, 0., 20.);
59   TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
60   TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
61   TH1F *hY = new TH1F("hY"," Dimuon rapidity distribution",100,-10,0);
62   TH1F *hYM = new TH1F("hYM"," matched track rapidity distribution",100,-10,0);
63   TH1F *hYF = new TH1F("hYF"," fake track rapidity distribution",100,-10,0);
64   TH1F *hEta = new TH1F("hEta"," Dimuon pseudo-rapidity distribution",100,-10,0);
65   TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
66   TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
67   TH1F *hPhi = new TH1F("hPhi"," Dimuon phi distribution",100,-1.,9.);
68   TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
69   TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
70   
71   // link to reconstructed and simulated tracks
72   AliMUONRecoCheck rc(esdFileName, SimDir);
73   
74   // load necessary data from OCDB
75   AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
76   AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
77   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
78   if (!recoParam) return;
79   
80   // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
81   Double_t sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
82   
83   TLorentzVector vMu1, vMu2, vDiMu;
84   
85   // Loop over ESD events
86   FirstEvent = TMath::Max(0, FirstEvent);
87   LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
88   for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
89     
90     // get reconstructed and simulated tracks
91     AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
92     AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
93     if (!muonTrackStore || !trackRefStore) continue;
94     
95     // loop over ESD tracks and flag them
96     const AliESDEvent* esd = rc.GetESDEvent();
97     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
98     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
99       
100       AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
101       
102       // skip ghosts
103       if (!esdTrack->ContainTrackerData()) continue;
104       
105       // find the corresponding MUON track
106       AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
107       
108       // try to match the reconstructed track with a simulated one
109       Float_t fractionOfMatchCluster = 0.;
110       AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*muonTrack, *trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut);
111       
112       // take actions according to matching result
113       if (matchedTrackRef) {
114         
115         // flag matched tracks
116         esdTrack->SetLabel(matchedTrackRef->GetUniqueID());
117         
118         // remove already matched trackRefs
119         trackRefStore->Remove(*matchedTrackRef);
120         
121       } else {
122         
123         // flag fake tracks
124         esdTrack->SetLabel(-1);
125         
126       }
127       
128     }
129     
130     // double loop over ESD tracks, build pairs and fill histograms according to their label
131     for (Int_t iTrack1 = 0; iTrack1 <  nTracks;  iTrack1++) {
132       AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
133       
134       // skip ghosts
135       if (!muonTrack1->ContainTrackerData()) continue;
136       
137       // get track info
138       Short_t charge1 = muonTrack1->Charge();
139       Int_t label1 = muonTrack1->GetLabel();
140       muonTrack1->LorentzP(vMu1);
141       
142       for (Int_t iTrack2 = iTrack1+1; iTrack2 <  nTracks;  iTrack2++) {
143         AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
144         
145         // skip ghosts
146         if (!muonTrack2->ContainTrackerData()) continue;
147         
148         // keep only opposite sign pairs
149         Short_t charge2 = muonTrack2->Charge();
150         if (charge1*charge2 > 0) continue;
151         
152         // get track info
153         Int_t label2 = muonTrack2->GetLabel();
154         muonTrack2->LorentzP(vMu2);
155         
156         // compute kinematics of the pair
157         vDiMu = vMu1 + vMu2;
158         Float_t mass = vDiMu.M();
159         Float_t p = vDiMu.P();
160         Float_t pt = vDiMu.Pt();
161         Float_t y = vDiMu.Rapidity();
162         Float_t eta = vDiMu.Eta();
163         Float_t phi = vDiMu.Phi();
164         if (phi < 0) phi += 2.*TMath::Pi();
165         
166         // fill global histograms
167         hMass->Fill(mass);
168         hP->Fill(p);
169         hPt->Fill(pt);
170         hY->Fill(y);
171         hEta->Fill(eta);
172         hPhi->Fill(phi);
173         
174         // fill histograms according to labels
175         if (label1 >= 0 && label2 >= 0) {
176           
177           hMassM->Fill(mass);
178           hPM->Fill(p);
179           hPtM->Fill(pt);
180           hYM->Fill(y);
181           hEtaM->Fill(eta);
182           hPhiM->Fill(phi);
183           
184         } else {
185           
186           hMassF->Fill(mass);
187           hPF->Fill(p);
188           hPtF->Fill(pt);
189           hYF->Fill(y);
190           hEtaF->Fill(eta);
191           hPhiF->Fill(phi);
192           
193         }
194         
195       }
196       
197     }
198       
199   } // end of loop over events
200   
201   // plot results
202   TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600);
203   cDiFakesSummary.Divide(3,2);
204   cDiFakesSummary.cd(1);
205   cDiFakesSummary.GetPad(1)->SetLogy();
206   hMass->Draw();
207   hMass->SetMinimum(0.5);
208   hMassM->Draw("same");
209   hMassM->SetLineColor(4);
210   hMassF->Draw("same");
211   hMassF->SetLineColor(2);
212   hMassF->SetFillColor(2);
213   hMassF->SetFillStyle(3017);
214   cDiFakesSummary.cd(2);
215   cDiFakesSummary.GetPad(3)->SetLogy();
216   hP->Draw();
217   hP->SetMinimum(0.5);
218   hPM->Draw("same");
219   hPM->SetLineColor(4);
220   hPF->Draw("same");
221   hPF->SetLineColor(2);
222   hPF->SetFillColor(2);
223   hPF->SetFillStyle(3017);
224   cDiFakesSummary.cd(3);
225   cDiFakesSummary.GetPad(4)->SetLogy();
226   hPt->Draw();
227   hPt->SetMinimum(0.5);
228   hPtM->Draw("same");
229   hPtM->SetLineColor(4);
230   hPtF->Draw("same");
231   hPtF->SetLineColor(2);
232   hPtF->SetFillColor(2);
233   hPtF->SetFillStyle(3017);
234   cDiFakesSummary.cd(4);
235   cDiFakesSummary.GetPad(2)->SetLogy();
236   hY->Draw();
237   hY->SetMinimum(0.5);
238   hYM->Draw("same");
239   hYM->SetLineColor(4);
240   hYF->Draw("same");
241   hYF->SetLineColor(2);
242   hYF->SetFillColor(2);
243   hYF->SetFillStyle(3017);
244   cDiFakesSummary.cd(5);
245   cDiFakesSummary.GetPad(5)->SetLogy();
246   hEta->Draw();
247   hEta->SetMinimum(0.5);
248   hEtaM->Draw("same");
249   hEtaM->SetLineColor(4);
250   hEtaF->Draw("same");
251   hEtaF->SetLineColor(2);
252   hEtaF->SetFillColor(2);
253   hEtaF->SetFillStyle(3017);
254   cDiFakesSummary.cd(6);
255   cDiFakesSummary.GetPad(6)->SetLogy();
256   hPhi->Draw();
257   hPhi->SetMinimum(0.5);
258   hPhiM->Draw("same");
259   hPhiM->SetLineColor(4);
260   hPhiF->Draw("same");
261   hPhiF->SetLineColor(2);
262   hPhiF->SetFillColor(2);
263   hPhiF->SetFillStyle(3017);
264   
265   // save results
266   histoFile->cd();
267   histoFile->Write();
268   cDiFakesSummary.Write();
269   histoFile->Close();
270   
271 }
272
273 //-----------------------------------------------------------------------
274 Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut)
275 {
276   /// Try to match 2 tracks
277   
278   Bool_t compTrack[10];
279   Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack);
280   fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters());
281   
282   if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2
283       (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5
284       fractionOfMatchCluster > 0.5) return kTRUE;                       // more than 50% of clusters matched
285   else return kFALSE;
286   
287 }
288
289 //-----------------------------------------------------------------------
290 AliMUONTrack* MatchWithTrackRef(AliMUONTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
291                                 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut)
292 {
293   /// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters
294   
295   AliMUONTrack *matchedTrackRef = 0x0;
296   fractionOfMatchCluster = 0.;
297   
298   if (useLabel) { // by using the MC label
299     
300     // get the corresponding simulated track if any
301     Int_t label = muonTrack.GetMCLabel();
302     matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label);
303     
304     // get the fraction of matched clusters
305     if (matchedTrackRef) {
306       Int_t nMatchClusters = 0;
307       Int_t nClusters = muonTrack.GetNClusters();
308       for (Int_t iCl = 0; iCl < nClusters; iCl++)
309         if (((AliMUONTrackParam*) muonTrack.GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
310           nMatchClusters++;
311       fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)nClusters);
312     }
313     
314   } else { // by comparing cluster/TrackRef positions
315     
316     // look for the corresponding simulated track if any
317     TIter next(trackRefStore.CreateIterator());
318     AliMUONTrack* trackRef;
319     while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
320       
321       // check compatibility
322       Float_t f = 0.;
323       if (TrackMatched(muonTrack, *trackRef, f, sigmaCut)) {
324         matchedTrackRef = trackRef;
325         fractionOfMatchCluster = f;
326         break;
327       }
328       
329     }
330     
331   }
332   
333   return matchedTrackRef;
334   
335 }
336