o add Reset function to CalPad and CalROC o Add functionality to AliTPCdataQA - Reset...
[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 Double_t sigmaCut = -1.;
37
38 //-----------------------------------------------------------------------
39 void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
40                  const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
41                  const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
42 {
43   
44   //Reset ROOT and connect tree file
45   gROOT->Reset();
46   
47   // File for histograms and histogram booking
48   TFile *histoFile = new TFile("DiFakes.root", "RECREATE");
49   
50   TH1F *hMass = new TH1F("hMass", "Dimuon mass distribution (GeV/c^{2})", 100, 0., 12.);
51   TH1F *hMassM = new TH1F("hMassM", "matched track mass distribution (GeV/c^{2})", 100, 0., 12.);
52   TH1F *hMassF = new TH1F("hMassF", "fake track mass distribution (GeV/c^{2})", 100, 0., 12.);
53   TH1F *hP = new TH1F("hP", "Dimuon P distribution (GeV/c)", 100, 0., 200.);
54   TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
55   TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
56   TH1F *hPt = new TH1F("hPt", "Dimuon Pt distribution (GeV/c)", 100, 0., 20.);
57   TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
58   TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
59   TH1F *hY = new TH1F("hY"," Dimuon rapidity distribution",100,-10,0);
60   TH1F *hYM = new TH1F("hYM"," matched track rapidity distribution",100,-10,0);
61   TH1F *hYF = new TH1F("hYF"," fake track rapidity distribution",100,-10,0);
62   TH1F *hEta = new TH1F("hEta"," Dimuon pseudo-rapidity distribution",100,-10,0);
63   TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
64   TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
65   TH1F *hPhi = new TH1F("hPhi"," Dimuon phi distribution",100,-1.,9.);
66   TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
67   TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
68   
69   // link to reconstructed and simulated tracks
70   AliMUONRecoCheck rc(esdFileName, SimDir);
71   
72   // load necessary data from OCDB
73   AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
74   AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
75   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
76   if (!recoParam) return;
77   
78   // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
79   sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
80   
81   TLorentzVector vMu1, vMu2, vDiMu;
82   
83   // Loop over ESD events
84   FirstEvent = TMath::Max(0, FirstEvent);
85   LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
86   for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
87     
88     // get reconstructed and simulated tracks
89     AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
90     AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
91     if (!muonTrackStore || !trackRefStore) continue;
92     
93     // loop over ESD tracks and flag them
94     const AliESDEvent* esd = rc.GetESDEvent();
95     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
96     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
97       
98       AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
99       
100       // skip ghosts
101       if (!esdTrack->ContainTrackerData()) continue;
102       
103       // find the corresponding MUON track
104       AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
105       
106       // try to match the reconstructed track with a simulated one
107       Int_t nMatchClusters = 0;
108       AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut);
109       
110       // take actions according to matching result
111       if (matchedTrackRef) {
112         
113         // flag matched tracks
114         esdTrack->SetLabel(matchedTrackRef->GetUniqueID());
115         
116         // remove already matched trackRefs
117         trackRefStore->Remove(*matchedTrackRef);
118         
119       } else {
120         
121         // flag fake tracks
122         esdTrack->SetLabel(-1);
123         
124       }
125       
126     }
127     
128     // double loop over ESD tracks, build pairs and fill histograms according to their label
129     for (Int_t iTrack1 = 0; iTrack1 <  nTracks;  iTrack1++) {
130       AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
131       
132       // skip ghosts
133       if (!muonTrack1->ContainTrackerData()) continue;
134       
135       // get track info
136       Short_t charge1 = muonTrack1->Charge();
137       Int_t label1 = muonTrack1->GetLabel();
138       muonTrack1->LorentzP(vMu1);
139       
140       for (Int_t iTrack2 = iTrack1+1; iTrack2 <  nTracks;  iTrack2++) {
141         AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
142         
143         // skip ghosts
144         if (!muonTrack2->ContainTrackerData()) continue;
145         
146         // keep only opposite sign pairs
147         Short_t charge2 = muonTrack2->Charge();
148         if (charge1*charge2 > 0) continue;
149         
150         // get track info
151         Int_t label2 = muonTrack2->GetLabel();
152         muonTrack2->LorentzP(vMu2);
153         
154         // compute kinematics of the pair
155         vDiMu = vMu1 + vMu2;
156         Float_t mass = vDiMu.M();
157         Float_t p = vDiMu.P();
158         Float_t pt = vDiMu.Pt();
159         Float_t y = vDiMu.Rapidity();
160         Float_t eta = vDiMu.Eta();
161         Float_t phi = vDiMu.Phi();
162         if (phi < 0) phi += 2.*TMath::Pi();
163         
164         // fill global histograms
165         hMass->Fill(mass);
166         hP->Fill(p);
167         hPt->Fill(pt);
168         hY->Fill(y);
169         hEta->Fill(eta);
170         hPhi->Fill(phi);
171         
172         // fill histograms according to labels
173         if (label1 >= 0 && label2 >= 0) {
174           
175           hMassM->Fill(mass);
176           hPM->Fill(p);
177           hPtM->Fill(pt);
178           hYM->Fill(y);
179           hEtaM->Fill(eta);
180           hPhiM->Fill(phi);
181           
182         } else {
183           
184           hMassF->Fill(mass);
185           hPF->Fill(p);
186           hPtF->Fill(pt);
187           hYF->Fill(y);
188           hEtaF->Fill(eta);
189           hPhiF->Fill(phi);
190           
191         }
192         
193       }
194       
195     }
196       
197   } // end of loop over events
198   
199   // plot results
200   TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600);
201   cDiFakesSummary.Divide(3,2);
202   cDiFakesSummary.cd(1);
203   cDiFakesSummary.GetPad(1)->SetLogy();
204   hMass->Draw();
205   hMass->SetMinimum(0.5);
206   hMassM->Draw("same");
207   hMassM->SetLineColor(4);
208   hMassF->Draw("same");
209   hMassF->SetLineColor(2);
210   hMassF->SetFillColor(2);
211   hMassF->SetFillStyle(3017);
212   cDiFakesSummary.cd(2);
213   cDiFakesSummary.GetPad(3)->SetLogy();
214   hP->Draw();
215   hP->SetMinimum(0.5);
216   hPM->Draw("same");
217   hPM->SetLineColor(4);
218   hPF->Draw("same");
219   hPF->SetLineColor(2);
220   hPF->SetFillColor(2);
221   hPF->SetFillStyle(3017);
222   cDiFakesSummary.cd(3);
223   cDiFakesSummary.GetPad(4)->SetLogy();
224   hPt->Draw();
225   hPt->SetMinimum(0.5);
226   hPtM->Draw("same");
227   hPtM->SetLineColor(4);
228   hPtF->Draw("same");
229   hPtF->SetLineColor(2);
230   hPtF->SetFillColor(2);
231   hPtF->SetFillStyle(3017);
232   cDiFakesSummary.cd(4);
233   cDiFakesSummary.GetPad(2)->SetLogy();
234   hY->Draw();
235   hY->SetMinimum(0.5);
236   hYM->Draw("same");
237   hYM->SetLineColor(4);
238   hYF->Draw("same");
239   hYF->SetLineColor(2);
240   hYF->SetFillColor(2);
241   hYF->SetFillStyle(3017);
242   cDiFakesSummary.cd(5);
243   cDiFakesSummary.GetPad(5)->SetLogy();
244   hEta->Draw();
245   hEta->SetMinimum(0.5);
246   hEtaM->Draw("same");
247   hEtaM->SetLineColor(4);
248   hEtaF->Draw("same");
249   hEtaF->SetLineColor(2);
250   hEtaF->SetFillColor(2);
251   hEtaF->SetFillStyle(3017);
252   cDiFakesSummary.cd(6);
253   cDiFakesSummary.GetPad(6)->SetLogy();
254   hPhi->Draw();
255   hPhi->SetMinimum(0.5);
256   hPhiM->Draw("same");
257   hPhiM->SetLineColor(4);
258   hPhiF->Draw("same");
259   hPhiF->SetLineColor(2);
260   hPhiF->SetFillColor(2);
261   hPhiF->SetFillStyle(3017);
262   
263   // save results
264   histoFile->cd();
265   histoFile->Write();
266   cDiFakesSummary.Write();
267   histoFile->Close();
268   
269 }
270