AliRoot Core  edcc906 (edcc906)
DIMUONFakes.C
Go to the documentation of this file.
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 
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  AliLog::SetClassDebugLevel("AliMCEvent",-1);
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
76  AliCDBManager::Instance()->SetSpecificStorage("GRP/GRP/Data","local://.");
79  if (!recoParam) return;
80 
81  // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
82  sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
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();
88 
89  TLorentzVector vMu1, vMu2, vDiMu;
90  Int_t nReconstructiblePairs = 0;
91  Int_t nReconstructedPairs = 0;
92  Int_t nMatchedPairs = 0;
93 
94  // Loop over ESD events
95  FirstEvent = TMath::Max(0, FirstEvent);
96  LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
97  for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
98 
99  // get reconstructed and simulated tracks
100  AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
101  AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
102  if (!muonTrackStore || !trackRefStore) continue;
103 
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 
116  // loop over ESD tracks and flag them
117  const AliESDEvent* esd = rc.GetESDEvent();
118  Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
119  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
120 
121  AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
122 
123  // skip ghosts
124  if (!esdTrack->ContainTrackerData()) continue;
125 
126  // find the corresponding MUON track
127  AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
128 
129  // try to match the reconstructed track with a simulated one
130  Int_t nMatchClusters = 0;
131  AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut);
132 
133  // take actions according to matching result
134  if (matchedTrackRef) {
135 
136  // flag matched tracks
137  esdTrack->SetLabel(matchedTrackRef->GetUniqueID());
138 
139  // remove already matched trackRefs
140  trackRefStore->Remove(*matchedTrackRef);
141 
142  } else {
143 
144  // flag fake tracks
145  esdTrack->SetLabel(-1);
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 
173  nReconstructedPairs++;
174 
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();
186  Float_t phi = vDiMu.Phi();
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 
200  nMatchedPairs++;
201 
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 
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 
303 }
304 
Base class of a track container.
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
void ImproveTracks(Bool_t flag)
switch on/off the track improvement and keep the default cut in sigma to apply on cluster (local chi2...
AliESDMuonTrack * GetMuonTrack(Int_t i)
void MakeMoreTrackCandidates(Bool_t flag)
switch on/off the building of track candidates starting from 1 cluster in each of the stations 4 and ...
Double_t GetSigmaCutForTracking() const
return the cut in sigma to apply on cluster (local chi2) and track (global chi2) during tracking ...
Class to describe the MUON tracks in the Event Summary Data class.
Int_t NumberOfEvents() const
Return the total number of events.
virtual TIterator * CreateIterator() const =0
Create an iterator to loop over tracks.
static void SetClassDebugLevel(const char *className, Int_t level)
Definition: AliLog.cxx:513
Int_t GetRunNumber()
Return the run number of the current ESD event.
virtual AliMUONTrack * Remove(AliMUONTrack &track)=0
Remove a track from the store.
TROOT * gROOT
Float_t p[]
Definition: kNNTest.C:133
Class with MUON reconstruction parameters.
Bool_t request2ChInSameSt45
Definition: MUONFakes.C:38
Bool_t ContainTrackerData() const
Utility class to check reconstruction.
AliMUONVTrackStore * TrackRefs(Int_t event)
Return reference muon tracks.
Double_t sigmaCut
Definition: DIMUONFakes.C:36
void SetLabel(Int_t label)
Set the corresponding MC track number.
void SetSpecificStorage(const char *calibType, const char *dbString, Int_t version=-1, Int_t subVersion=-1)
Int_t GetNumberOfMuonTracks() const
Definition: AliESDEvent.h:543
Short_t Charge() const
void LorentzP(TLorentzVector &vP) const
static AliMUONTrack * FindCompatibleTrack(AliMUONTrack &track, AliMUONVTrackStore &trackStore, Int_t &nMatchClusters, Bool_t useLabel=kFALSE, Double_t sigmaCut=10.)
Return the track from the store matched with the given track (or 0x0) and the fraction of matched clu...
void SetRun(Int_t run)
void SetDefaultStorage(const char *dbString)
Int_t GetLabel() const
Return the corresponding MC track number.
AliMUONVTrackStore * ReconstructedTracks(Int_t event, Bool_t refit=kTRUE)
Return the list of reconstructed tracks.
UInt_t requestedStationMask
Definition: MUONFakes.C:37
void RequestStation(Int_t iSt, Bool_t flag)
request or not at least one cluster in the station to validate the track
Reconstructed track in ALICE dimuon spectrometer.
Definition: AliMUONTrack.h:24
Double_t GetSigmaCutForImprovement() const
return the cut in sigma to apply on cluster (local chi2) during track improvement ...
static AliCDBManager * Instance(TMap *entryCache=NULL, Int_t run=-1)
void DIMUONFakes(Bool_t useLabel=kFALSE, Int_t FirstEvent=0, Int_t LastEvent=-1, const TString esdFileName="AliESDs.root", const TString SimDir="./generated/", const TString ocdbPath="local://$ALICE_ROOT/OCDB")
Definition: DIMUONFakes.C:39
virtual TObject * FindObject(const char *name) const
Find an object by name.
const AliESDEvent * GetESDEvent() const
Return the reconstructed data of current event.
AliMUONRecoParam * LoadRecoParam()
Definition: AliMUONCDB.cxx:535