1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
19 /// \file MUONRecoCheck.C
20 /// \brief Utility macro to check the muon reconstruction.
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.
26 /// \author Jean-Pierre Cussonneau, Subatech
30 #include "TClonesArray.h"
33 #include "TGraphErrors.h"
38 #include "AliCDBManager.h"
41 #include "AliMUONCDB.h"
42 #include "AliMUONConstants.h"
43 #include "AliMUONTrack.h"
44 #include "AliMUONRecoCheck.h"
45 #include "AliMUONTrackParam.h"
46 #include "AliMUONRecoParam.h"
47 #include "AliMUONVTrackStore.h"
48 #include "AliMUONVCluster.h"
50 void MUONRecoCheck (Int_t nEvent = -1, char * pathSim="./generated/", char * esdFileName="AliESDs.root",
51 const char* ocdbPath = "local://$ALICE_ROOT/OCDB")
54 // File for histograms and histogram booking
55 TFile *histoFile = new TFile("MUONRecoCheck.root", "RECREATE");
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);
59 TH1F *hNClusterComp = new TH1F("hNClusterComp"," Nb of compatible clusters / track ",15,-0.5,14.5);
60 TH1F *hTrackRefID = new TH1F("hTrackRefID"," track reference ID ",100,-0.5,99.5);
62 // momentum resolution
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);
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");
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);
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);
101 AliMUONRecoCheck rc(esdFileName, pathSim);
103 // load necessary data from OCDB
104 AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
105 AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
106 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
107 if (!recoParam) return;
109 // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
110 Double_t sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
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();
117 Int_t nevents = rc.NumberOfEvents();
119 if (nevents < nEvent || nEvent < 0) nEvent = nevents;
122 Int_t nReconstructibleTracks = 0;
123 Int_t nReconstructedTracks = 0;
124 Int_t nReconstructibleTracksCheck = 0;
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;
129 for (ievent=0; ievent<nEvent; ievent++)
131 if (!(ievent%10)) printf(" **** event # %d \n",ievent);
133 AliMUONVTrackStore* trackStore = rc.ReconstructedTracks(ievent, kFALSE);
134 AliMUONVTrackStore* trackRefStore = rc.ReconstructibleTracks(ievent, requestedStationMask, request2ChInSameSt45);
136 hReconstructible->Fill(trackRefStore->GetSize());
137 hReco->Fill(trackStore->GetSize());
139 nReconstructibleTracks += trackRefStore->GetSize();
140 nReconstructedTracks += trackStore->GetSize();
142 // loop over trackRef
143 TIter next(trackRefStore->CreateIterator());
144 AliMUONTrack* trackRef;
145 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) )
148 hTrackRefID->Fill(trackRef->GetMCLabel());
150 AliMUONTrack* trackMatched = 0x0;
151 Int_t nMatchClusters = 0;
153 // loop over trackReco and look for compatible track
154 TIter next2(trackStore->CreateIterator());
155 AliMUONTrack* trackReco;
156 while ( ( trackReco = static_cast<AliMUONTrack*>(next2()) ) )
159 // check if trackReco is compatible with trackRef
161 if (trackReco->Match(*trackRef, sigmaCut, nMatchClusters)) {
162 trackMatched = trackReco;
169 if (trackMatched) { // tracking requirements verified, track is found
170 nReconstructibleTracksCheck++;
171 hNClusterComp->Fill(nMatchClusters);
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();
180 pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
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);
183 trackParam = trackMatched->GetTrackParamAtVertex();
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();
191 pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
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);
194 hResMomVertex->Fill(p2-p1);
195 hResMomVertexVsMom->Fill(p1,p2-p1);
197 trackParam = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
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();
205 pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
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);
208 trackParam = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
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();
216 pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
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);
219 hResMomFirstCluster->Fill(p2-p1);
220 hResMomFirstClusterVsMom->Fill(p1,p2-p1);
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());
235 trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->After(trackParamAtCluster2);
237 trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->After(trackParamAtCluster1);
242 } // end loop track ref.
244 } // end loop on event
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);
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);
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);
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);
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);
287 gResMomVertexVsMom->Write();
288 gResMomFirstClusterVsMom->Write();
289 gResidualXPerChMean->Write();
290 gResidualXPerChSigma->Write();
291 gResidualYPerChMean->Write();
292 gResidualYPerChSigma->Write();