Converting PWG/TRD to native cmake
[u/mrichter/AliRoot.git] / MUON / DIMUONFakes.C
CommitLineData
30b3fe0f 1#if !defined(__CINT__) || defined(__MAKECINT__)
2// ROOT includes
30b3fe0f 3#include <TFile.h>
4#include <TH1.h>
5#include <TCanvas.h>
6#include <Riostream.h>
7#include <TROOT.h>
8#include <TClonesArray.h>
30b3fe0f 9#include <TLorentzVector.h>
10
11// STEER includes
12#include "AliLog.h"
30b3fe0f 13#include "AliESDEvent.h"
14#include "AliESDMuonTrack.h"
30b3fe0f 15#include "AliCDBManager.h"
16
17// MUON includes
a99c3449 18#include "AliMUONCDB.h"
30b3fe0f 19#include "AliMUONTrack.h"
20#include "AliMUONVTrackStore.h"
21#include "AliMUONTrackParam.h"
30b3fe0f 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
f202486b 36Double_t sigmaCut = -1.;
30b3fe0f 37
38//-----------------------------------------------------------------------
39void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
a99c3449 40 const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
41 const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
30b3fe0f 42{
43
44 //Reset ROOT and connect tree file
45 gROOT->Reset();
46
590a335d 47 AliLog::SetClassDebugLevel("AliMCEvent",-1);
48
30b3fe0f 49 // File for histograms and histogram booking
50 TFile *histoFile = new TFile("DiFakes.root", "RECREATE");
51
23035434 52 TH1F *hMass = new TH1F("hMass", "Dimuon mass distribution (GeV/c^{2})", 100, 0., 12.);
30b3fe0f 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.);
23035434 55 TH1F *hP = new TH1F("hP", "Dimuon P distribution (GeV/c)", 100, 0., 200.);
30b3fe0f 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.);
23035434 58 TH1F *hPt = new TH1F("hPt", "Dimuon Pt distribution (GeV/c)", 100, 0., 20.);
30b3fe0f 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.);
23035434 61 TH1F *hY = new TH1F("hY"," Dimuon rapidity distribution",100,-10,0);
30b3fe0f 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);
23035434 64 TH1F *hEta = new TH1F("hEta"," Dimuon pseudo-rapidity distribution",100,-10,0);
30b3fe0f 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);
23035434 67 TH1F *hPhi = new TH1F("hPhi"," Dimuon phi distribution",100,-1.,9.);
30b3fe0f 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
a99c3449 71 // link to reconstructed and simulated tracks
72 AliMUONRecoCheck rc(esdFileName, SimDir);
30b3fe0f 73
a99c3449 74 // load necessary data from OCDB
75 AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
590a335d 76 AliCDBManager::Instance()->SetSpecificStorage("GRP/GRP/Data","local://.");
a99c3449 77 AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
78 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
79 if (!recoParam) return;
30b3fe0f 80
a99c3449 81 // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
f202486b 82 sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
590a335d 83 // compute the mask of requested stations from recoParam
84 UInt_t requestedStationMask = 0;
85 for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
86 // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
87 Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
30b3fe0f 88
89 TLorentzVector vMu1, vMu2, vDiMu;
590a335d 90 Int_t nReconstructiblePairs = 0;
91 Int_t nReconstructedPairs = 0;
92 Int_t nMatchedPairs = 0;
30b3fe0f 93
94 // Loop over ESD events
95 FirstEvent = TMath::Max(0, FirstEvent);
a99c3449 96 LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
30b3fe0f 97 for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
98
a99c3449 99 // get reconstructed and simulated tracks
100 AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
30b3fe0f 101 AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
a99c3449 102 if (!muonTrackStore || !trackRefStore) continue;
30b3fe0f 103
590a335d 104 // count the number of reconstructible pairs
105 Int_t nMuPlus = 0, nMuMinus = 0;
106 TIter next(trackRefStore->CreateIterator());
107 AliMUONTrack* trackRef;
108 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
109 if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
110 if (trackRef->GetTrackParamAtVertex()->GetCharge() > 0) nMuPlus++;
111 else nMuMinus++;
112 }
113 }
114 nReconstructiblePairs += nMuPlus*nMuMinus;
115
30b3fe0f 116 // loop over ESD tracks and flag them
a99c3449 117 const AliESDEvent* esd = rc.GetESDEvent();
30b3fe0f 118 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
119 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
120
a99c3449 121 AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
30b3fe0f 122
123 // skip ghosts
a99c3449 124 if (!esdTrack->ContainTrackerData()) continue;
125
126 // find the corresponding MUON track
127 AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
30b3fe0f 128
129 // try to match the reconstructed track with a simulated one
f202486b 130 Int_t nMatchClusters = 0;
131 AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut);
30b3fe0f 132
133 // take actions according to matching result
134 if (matchedTrackRef) {
135
136 // flag matched tracks
a99c3449 137 esdTrack->SetLabel(matchedTrackRef->GetUniqueID());
30b3fe0f 138
139 // remove already matched trackRefs
140 trackRefStore->Remove(*matchedTrackRef);
141
142 } else {
143
144 // flag fake tracks
a99c3449 145 esdTrack->SetLabel(-1);
30b3fe0f 146
147 }
148
149 }
150
151 // double loop over ESD tracks, build pairs and fill histograms according to their label
152 for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
153 AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
154
155 // skip ghosts
156 if (!muonTrack1->ContainTrackerData()) continue;
157
158 // get track info
159 Short_t charge1 = muonTrack1->Charge();
160 Int_t label1 = muonTrack1->GetLabel();
161 muonTrack1->LorentzP(vMu1);
162
163 for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
164 AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
165
166 // skip ghosts
167 if (!muonTrack2->ContainTrackerData()) continue;
168
169 // keep only opposite sign pairs
170 Short_t charge2 = muonTrack2->Charge();
171 if (charge1*charge2 > 0) continue;
172
590a335d 173 nReconstructedPairs++;
174
30b3fe0f 175 // get track info
176 Int_t label2 = muonTrack2->GetLabel();
177 muonTrack2->LorentzP(vMu2);
178
179 // compute kinematics of the pair
180 vDiMu = vMu1 + vMu2;
181 Float_t mass = vDiMu.M();
182 Float_t p = vDiMu.P();
183 Float_t pt = vDiMu.Pt();
184 Float_t y = vDiMu.Rapidity();
185 Float_t eta = vDiMu.Eta();
23035434 186 Float_t phi = vDiMu.Phi();
30b3fe0f 187 if (phi < 0) phi += 2.*TMath::Pi();
188
189 // fill global histograms
190 hMass->Fill(mass);
191 hP->Fill(p);
192 hPt->Fill(pt);
193 hY->Fill(y);
194 hEta->Fill(eta);
195 hPhi->Fill(phi);
196
197 // fill histograms according to labels
198 if (label1 >= 0 && label2 >= 0) {
199
590a335d 200 nMatchedPairs++;
201
30b3fe0f 202 hMassM->Fill(mass);
203 hPM->Fill(p);
204 hPtM->Fill(pt);
205 hYM->Fill(y);
206 hEtaM->Fill(eta);
207 hPhiM->Fill(phi);
208
209 } else {
210
211 hMassF->Fill(mass);
212 hPF->Fill(p);
213 hPtF->Fill(pt);
214 hYF->Fill(y);
215 hEtaF->Fill(eta);
216 hPhiF->Fill(phi);
217
218 }
219
220 }
221
222 }
223
224 } // end of loop over events
225
226 // plot results
227 TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600);
228 cDiFakesSummary.Divide(3,2);
229 cDiFakesSummary.cd(1);
230 cDiFakesSummary.GetPad(1)->SetLogy();
231 hMass->Draw();
232 hMass->SetMinimum(0.5);
233 hMassM->Draw("same");
234 hMassM->SetLineColor(4);
235 hMassF->Draw("same");
236 hMassF->SetLineColor(2);
237 hMassF->SetFillColor(2);
238 hMassF->SetFillStyle(3017);
239 cDiFakesSummary.cd(2);
240 cDiFakesSummary.GetPad(3)->SetLogy();
241 hP->Draw();
242 hP->SetMinimum(0.5);
243 hPM->Draw("same");
244 hPM->SetLineColor(4);
245 hPF->Draw("same");
246 hPF->SetLineColor(2);
247 hPF->SetFillColor(2);
248 hPF->SetFillStyle(3017);
249 cDiFakesSummary.cd(3);
250 cDiFakesSummary.GetPad(4)->SetLogy();
251 hPt->Draw();
252 hPt->SetMinimum(0.5);
253 hPtM->Draw("same");
254 hPtM->SetLineColor(4);
255 hPtF->Draw("same");
256 hPtF->SetLineColor(2);
257 hPtF->SetFillColor(2);
258 hPtF->SetFillStyle(3017);
259 cDiFakesSummary.cd(4);
260 cDiFakesSummary.GetPad(2)->SetLogy();
261 hY->Draw();
262 hY->SetMinimum(0.5);
263 hYM->Draw("same");
264 hYM->SetLineColor(4);
265 hYF->Draw("same");
266 hYF->SetLineColor(2);
267 hYF->SetFillColor(2);
268 hYF->SetFillStyle(3017);
269 cDiFakesSummary.cd(5);
270 cDiFakesSummary.GetPad(5)->SetLogy();
271 hEta->Draw();
272 hEta->SetMinimum(0.5);
273 hEtaM->Draw("same");
274 hEtaM->SetLineColor(4);
275 hEtaF->Draw("same");
276 hEtaF->SetLineColor(2);
277 hEtaF->SetFillColor(2);
278 hEtaF->SetFillStyle(3017);
279 cDiFakesSummary.cd(6);
280 cDiFakesSummary.GetPad(6)->SetLogy();
281 hPhi->Draw();
282 hPhi->SetMinimum(0.5);
283 hPhiM->Draw("same");
284 hPhiM->SetLineColor(4);
285 hPhiF->Draw("same");
286 hPhiF->SetLineColor(2);
287 hPhiF->SetFillColor(2);
288 hPhiF->SetFillStyle(3017);
289
290 // save results
291 histoFile->cd();
292 histoFile->Write();
293 cDiFakesSummary.Write();
294 histoFile->Close();
295
590a335d 296 // print results
297 printf("\n");
298 printf("nb of reconstructible OS pairs: %d \n", nReconstructiblePairs);
299 printf("nb of reconstructed OS pairs: %d \n", nReconstructedPairs);
300 printf("nb of reconstructed OS pairs matched with trackRefs: %d \n", nMatchedPairs);
301 printf("\n");
302
30b3fe0f 303}
304