]>
Commit | Line | Data |
---|---|---|
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 | 50 | void 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 |