Undefined method removed.
[u/mrichter/AliRoot.git] / MUON / MUONRecoCheck.C
CommitLineData
b8dc484b 1/**************************************************************************
48217459 2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use, copy, modify and distribute this software and its *
8* documentation strictly for non-commercial purposes is hereby granted *
9* without fee, provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
b8dc484b 15
e54bf126 16// $Id$
17
18/// \ingroup macros
19/// \file MUONRecoCheck.C
20/// \brief Utility macro to check the muon reconstruction.
21///
22/// Reconstructed tracks are compared to reference tracks. The reference tracks
23/// are built from AliTrackReference for the hit in chamber (0..9) and from
24/// kinematics (TreeK) for the vertex parameters.
25///
26/// \author Jean-Pierre Cussonneau, Subatech
27
b8dc484b 28// ROOT includes
f202486b 29#include "TMath.h"
b8dc484b 30#include "TClonesArray.h"
31#include "TH1.h"
f202486b 32#include "TH2.h"
33#include "TGraphErrors.h"
34#include "TF1.h"
b8dc484b 35#include "TFile.h"
36
37// STEER includes
a99c3449 38#include "AliCDBManager.h"
b8dc484b 39
40// MUON includes
a99c3449 41#include "AliMUONCDB.h"
b8dc484b 42#include "AliMUONConstants.h"
43#include "AliMUONTrack.h"
44#include "AliMUONRecoCheck.h"
45#include "AliMUONTrackParam.h"
a99c3449 46#include "AliMUONRecoParam.h"
48217459 47#include "AliMUONVTrackStore.h"
f202486b 48#include "AliMUONVCluster.h"
b8dc484b 49
a99c3449 50void MUONRecoCheck (Int_t nEvent = -1, char * pathSim="./generated/", char * esdFileName="AliESDs.root",
f202486b 51 const char* ocdbPath = "local://$ALICE_ROOT/OCDB")
e54bf126 52{
b8dc484b 53
b8dc484b 54 // File for histograms and histogram booking
55 TFile *histoFile = new TFile("MUONRecoCheck.root", "RECREATE");
48217459 56
b8dc484b 57 TH1F *hReconstructible = new TH1F("hReconstructible"," Nb of reconstructible tracks ",15,-0.5,14.5);
58 TH1F *hReco = new TH1F("hReco"," Nb of reconstructed tracks / evt",15,-0.5,14.5);
96ebe67e 59 TH1F *hNClusterComp = new TH1F("hNClusterComp"," Nb of compatible clusters / track ",15,-0.5,14.5);
b8dc484b 60 TH1F *hTrackRefID = new TH1F("hTrackRefID"," track reference ID ",100,-0.5,99.5);
61
f202486b 62 // momentum resolution
96ebe67e 63 TH1F *hResMomVertex = new TH1F("hResMomVertex"," delta P vertex (GeV/c)",100,-10.,10);
64 TH1F *hResMomFirstCluster = new TH1F("hResMomFirstCluster"," delta P first cluster (GeV/c)",100,-10.,10);
f202486b 65 TH2D *hResMomVertexVsMom = new TH2D("hResMomVertexVsMom","#Delta_{p} at vertex versus p (GeV/c)",30,0.,300.,800,-20.,20.);
66 TH2D *hResMomFirstClusterVsMom = new TH2D("hResMomFirstClusterVsMom","#Delta_{p} at first cluster versus p (GeV/c)",30,0.,300.,800,-20.,20.);
67 TGraphErrors* gResMomVertexVsMom = new TGraphErrors(30);
68 gResMomVertexVsMom->SetName("gResMomVertexVsMom");
69 gResMomVertexVsMom->SetTitle("#Delta_{p}/p at vertex versus p;p (GeV/c);#sigma_{p}/p (%)");
70 TGraphErrors* gResMomFirstClusterVsMom = new TGraphErrors(30);
71 gResMomFirstClusterVsMom->SetName("gResMomFirstClusterVsMom");
72 gResMomFirstClusterVsMom->SetTitle("#Delta_{p}/p at first cluster versus p;p (GeV/c);#sigma_{p}/p (%)");
73 TF1* f = new TF1("f","gaus");
74
75 // residuals at clusters
76 histoFile->mkdir("residuals","residuals");
77 histoFile->cd("residuals");
78 TH1F* hResidualXInCh[AliMUONConstants::NTrackingCh()];
79 TH1F* hResidualYInCh[AliMUONConstants::NTrackingCh()];
80 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
81 hResidualXInCh[i] = new TH1F(Form("hResidualXInCh%d",i+1), Form("cluster-track residual-X distribution in chamber %d;#Delta_{X} (cm)",i+1), 1000, -1., 1.);
82 hResidualYInCh[i] = new TH1F(Form("hResidualYInCh%d",i+1), Form("cluster-track residual-Y distribution in chamber %d;#Delta_{Y} (cm)",i+1), 1000, -0.5, 0.5);
83 }
84 TGraphErrors* gResidualXPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh());
85 gResidualXPerChMean->SetName("gResidualXPerChMean");
86 gResidualXPerChMean->SetTitle("cluster-track residual-X per Ch: mean;chamber ID;<#Delta_{Y}> (cm)");
87 gResidualXPerChMean->SetMarkerStyle(kFullDotLarge);
88 TGraphErrors* gResidualYPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh());
89 gResidualYPerChMean->SetName("gResidualYPerChMean");
90 gResidualYPerChMean->SetTitle("cluster-track residual-Y per Ch: mean;chamber ID;<#Delta_{Y}> (cm)");
91 gResidualYPerChMean->SetMarkerStyle(kFullDotLarge);
92 TGraphErrors* gResidualXPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh());
93 gResidualXPerChSigma->SetName("gResidualXPerChSigma");
94 gResidualXPerChSigma->SetTitle("cluster-track residual-X per Ch: sigma;chamber ID;#sigma_{X} (cm)");
95 gResidualXPerChSigma->SetMarkerStyle(kFullDotLarge);
96 TGraphErrors* gResidualYPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh());
97 gResidualYPerChSigma->SetName("gResidualYPerChSigma");
98 gResidualYPerChSigma->SetTitle("cluster-track residual-Y per Ch: sigma;chamber ID;#sigma_{X} (cm)");
99 gResidualYPerChSigma->SetMarkerStyle(kFullDotLarge);
48217459 100
a99c3449 101 AliMUONRecoCheck rc(esdFileName, pathSim);
102
103 // load necessary data from OCDB
104 AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
105 AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
a99c3449 106 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
107 if (!recoParam) return;
a0bfad0e 108
f202486b 109 // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
a99c3449 110 Double_t sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
f202486b 111 // compute the mask of requested stations from recoParam
112 UInt_t requestedStationMask = 0;
113 for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
114 // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
115 Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
48217459 116
117 Int_t nevents = rc.NumberOfEvents();
b8dc484b 118
edcc5b56 119 if (nevents < nEvent || nEvent < 0) nEvent = nevents;
b8dc484b 120
121 Int_t ievent;
122 Int_t nReconstructibleTracks = 0;
123 Int_t nReconstructedTracks = 0;
124 Int_t nReconstructibleTracksCheck = 0;
f202486b 125 AliMUONTrackParam *trackParam;
126 Double_t x1,y1,z1,pX1,pY1,pZ1,p1,pT1;
127 Double_t x2,y2,z2,pX2,pY2,pZ2,p2,pT2;
48217459 128
129 for (ievent=0; ievent<nEvent; ievent++)
130 {
b8dc484b 131 if (!(ievent%10)) printf(" **** event # %d \n",ievent);
b8dc484b 132
a99c3449 133 AliMUONVTrackStore* trackStore = rc.ReconstructedTracks(ievent, kFALSE);
f202486b 134 AliMUONVTrackStore* trackRefStore = rc.ReconstructibleTracks(ievent, requestedStationMask, request2ChInSameSt45);
b8dc484b 135
48217459 136 hReconstructible->Fill(trackRefStore->GetSize());
137 hReco->Fill(trackStore->GetSize());
138
139 nReconstructibleTracks += trackRefStore->GetSize();
140 nReconstructedTracks += trackStore->GetSize();
141
f202486b 142 // loop over trackRef
48217459 143 TIter next(trackRefStore->CreateIterator());
144 AliMUONTrack* trackRef;
48217459 145 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) )
146 {
48217459 147
f202486b 148 hTrackRefID->Fill(trackRef->GetMCLabel());
149
150 AliMUONTrack* trackMatched = 0x0;
151 Int_t nMatchClusters = 0;
152
153 // loop over trackReco and look for compatible track
48217459 154 TIter next2(trackStore->CreateIterator());
155 AliMUONTrack* trackReco;
48217459 156 while ( ( trackReco = static_cast<AliMUONTrack*>(next2()) ) )
157 {
f202486b 158
159 // check if trackReco is compatible with trackRef
160 Int_t n = 0;
161 if (trackReco->Match(*trackRef, sigmaCut, nMatchClusters)) {
162 trackMatched = trackReco;
163 nMatchClusters = n;
164 break;
61fed964 165 }
f202486b 166
48217459 167 }
168
f202486b 169 if (trackMatched) { // tracking requirements verified, track is found
48217459 170 nReconstructibleTracksCheck++;
f202486b 171 hNClusterComp->Fill(nMatchClusters);
48217459 172 trackParam = trackRef->GetTrackParamAtVertex();
173 x1 = trackParam->GetNonBendingCoor();
174 y1 = trackParam->GetBendingCoor();
175 z1 = trackParam->GetZ();
176 pX1 = trackParam->Px();
177 pY1 = trackParam->Py();
178 pZ1 = trackParam->Pz();
179 p1 = trackParam->P();
f202486b 180 pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
48217459 181
182 // printf(" Ref. track at vertex: x,y,z: %f %f %f px,py,pz,p: %f %f %f %f \n",x1,y1,z1,pX1,pY1,pZ1,p1);
f202486b 183 trackParam = trackMatched->GetTrackParamAtVertex();
48217459 184 x2 = trackParam->GetNonBendingCoor();
185 y2 = trackParam->GetBendingCoor();
186 z2 = trackParam->GetZ();
187 pX2 = trackParam->Px();
188 pY2 = trackParam->Py();
189 pZ2 = trackParam->Pz();
190 p2 = trackParam->P();
f202486b 191 pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
48217459 192 // printf(" Reconst. track at vertex: x,y,z: %f %f %f px,py,pz: %f %f %f %f \n",x2,y2,z2,pX2,pY2,pZ2,p2);
193
194 hResMomVertex->Fill(p2-p1);
f202486b 195 hResMomVertexVsMom->Fill(p1,p2-p1);
196
96ebe67e 197 trackParam = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
48217459 198 x1 = trackParam->GetNonBendingCoor();
199 y1 = trackParam->GetBendingCoor();
200 z1 = trackParam->GetZ();
201 pX1 = trackParam->Px();
202 pY1 = trackParam->Py();
203 pZ1 = trackParam->Pz();
204 p1 = trackParam->P();
f202486b 205 pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
206
96ebe67e 207 // printf(" Ref. track at 1st cluster: x,y,z: %f %f %f px,py,pz: %f %f %f \n",x1,y1,z1,pX1,pY1,pZ1);
f202486b 208 trackParam = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
48217459 209 x2 = trackParam->GetNonBendingCoor();
210 y2 = trackParam->GetBendingCoor();
211 z2 = trackParam->GetZ();
212 pX2 = trackParam->Px();
213 pY2 = trackParam->Py();
214 pZ2 = trackParam->Pz();
215 p2 = trackParam->P();
f202486b 216 pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
96ebe67e 217 // printf(" Reconst. track at 1st cluster: x,y,z: %f %f %f px,py,pz: %f %f %f \n",x2,y2,z2,pX2,pY2,pZ2);
48217459 218
96ebe67e 219 hResMomFirstCluster->Fill(p2-p1);
f202486b 220 hResMomFirstClusterVsMom->Fill(p1,p2-p1);
221
222 // Fill residuals
223 // Loop over clusters of first track
224 AliMUONTrackParam* trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
225 while (trackParamAtCluster1) {
226 AliMUONVCluster* cluster1 = trackParamAtCluster1->GetClusterPtr();
227 AliMUONTrackParam* trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
228 while (trackParamAtCluster2) {
229 AliMUONVCluster* cluster2 = trackParamAtCluster2->GetClusterPtr();
230 if (cluster1->GetDetElemId() == cluster2->GetDetElemId()) {
231 hResidualXInCh[cluster1->GetChamberId()]->Fill(cluster1->GetX() - cluster2->GetX());
232 hResidualYInCh[cluster1->GetChamberId()]->Fill(cluster1->GetY() - cluster2->GetY());
233 break;
234 }
235 trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->After(trackParamAtCluster2);
236 }
237 trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->After(trackParamAtCluster1);
238 }
239
b8dc484b 240 }
f202486b 241
b8dc484b 242 } // end loop track ref.
b8dc484b 243
48217459 244 } // end loop on event
245
f202486b 246 // compute momentum resolution versus p
247 for (Int_t i = 1; i <= hResMomVertexVsMom->GetNbinsX(); i++) {
248 TH1D *tmp = hResMomVertexVsMom->ProjectionY("tmp",i,i,"e");
249 Double_t p = hResMomVertexVsMom->GetBinCenter(i);
250 gResMomVertexVsMom->SetPoint(i-1,p,100.*tmp->GetRMS()/p);
251 gResMomVertexVsMom->SetPointError(i-1,hResMomVertexVsMom->GetBinWidth(i)/2.,100.*tmp->GetRMSError()/p);
252 delete tmp;
253 }
254 for (Int_t i = 1; i <= hResMomFirstClusterVsMom->GetNbinsX(); i++) {
255 TH1D *tmp = hResMomFirstClusterVsMom->ProjectionY("tmp",i,i,"e");
256 Double_t p = hResMomFirstClusterVsMom->GetBinCenter(i);
257 f->SetParameters(1.,0.,1.);
258 f->SetParLimits(0,0.,1.e3);
259 tmp->Fit("f","WWN");
260 gResMomFirstClusterVsMom->SetPoint(i-1,p,100.*f->GetParameter(2)/p);
261 gResMomFirstClusterVsMom->SetPointError(i-1,hResMomFirstClusterVsMom->GetBinWidth(i)/2.,100.*f->GetParError(2)/p);
262 delete tmp;
263 }
264
265 // compute residual mean and dispersion
266 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
267 hResidualXInCh[i]->GetXaxis()->SetRangeUser(-3.*hResidualXInCh[i]->GetRMS(), 3.*hResidualXInCh[i]->GetRMS());
268 gResidualXPerChMean->SetPoint(i, i+1, hResidualXInCh[i]->GetMean());
269 gResidualXPerChMean->SetPointError(i, 0., hResidualXInCh[i]->GetMeanError());
270 gResidualXPerChSigma->SetPoint(i, i+1, hResidualXInCh[i]->GetRMS());
271 gResidualXPerChSigma->SetPointError(i, 0., hResidualXInCh[i]->GetRMSError());
272 hResidualXInCh[i]->GetXaxis()->SetRange(0,0);
273 hResidualYInCh[i]->GetXaxis()->SetRangeUser(-3.*hResidualYInCh[i]->GetRMS(), 3.*hResidualYInCh[i]->GetRMS());
274 gResidualYPerChMean->SetPoint(i, i+1, hResidualYInCh[i]->GetMean());
275 gResidualYPerChMean->SetPointError(i, 0., hResidualYInCh[i]->GetMeanError());
276 gResidualYPerChSigma->SetPoint(i, i+1, hResidualYInCh[i]->GetRMS());
277 gResidualYPerChSigma->SetPointError(i, 0., hResidualYInCh[i]->GetRMSError());
278 hResidualYInCh[i]->GetXaxis()->SetRange(0,0);
279 }
280
b8dc484b 281 printf(" nb of reconstructible tracks: %d \n", nReconstructibleTracks);
282 printf(" nb of reconstructed tracks: %d \n", nReconstructedTracks);
283 printf(" nb of reconstructible tracks which are reconstructed: %d \n", nReconstructibleTracksCheck);
284
285 histoFile->Write();
f202486b 286 histoFile->cd();
287 gResMomVertexVsMom->Write();
288 gResMomFirstClusterVsMom->Write();
289 gResidualXPerChMean->Write();
290 gResidualXPerChSigma->Write();
291 gResidualYPerChMean->Write();
292 gResidualYPerChSigma->Write();
b8dc484b 293 histoFile->Close();
294}
295