Class description correction/addition/cleanup (Ivana)
[u/mrichter/AliRoot.git] / MUON / AliMUONRecoCheck.cxx
CommitLineData
b8dc484b 1/**************************************************************************
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 **************************************************************************/
15
16//////////////////////////////////////////////////////////////////////////////////////
17//
18// Utility class to check the muon reconstruction. Reconstructed tracks are compared
19// to reference tracks. The reference tracks are built from AliTrackReference for the
20// hit in chamber (0..9) and from kinematics for the vertex parameters.
21//
22//////////////////////////////////////////////////////////////////////////////////////
23
24
25#include <TParticle.h>
26
27#include "AliRun.h" // for gAlice
28#include "AliLog.h"
29#include "AliLoader.h"
30#include "AliRunLoader.h"
31#include "AliTrackReference.h"
32#include "AliHeader.h"
33#include "AliMC.h"
34#include "AliStack.h"
35#include "AliMUON.h"
36#include "AliMUONRecoCheck.h"
37#include "AliMUONTrack.h"
38#include "AliMUONData.h"
b8dc484b 39#include "AliMUONConstants.h"
40
41ClassImp(AliMUONRecoCheck)
42
43//_____________________________________________________________________________
44AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader)
45{
46// Constructor
47 fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
48
49 // open the run loader
50 fRunLoader = AliRunLoader::Open(chLoader);
51 if (!fRunLoader) {
52 AliError(Form("no run loader found in file %s","galice.root" ));
53 return;
54 }
55 // initialize loader's
56 AliLoader *loader = fRunLoader->GetLoader("MUONLoader");
57
58 // initialize container
59 fMUONData = new AliMUONData(loader,"MUON","MUON");
60
61 // Loading AliRun master
62 if (fRunLoader->GetAliRun() == 0x0) fRunLoader->LoadgAlice();
63
64 fRunLoader->LoadKinematics("READ");
65 fRunLoader->LoadTrackRefs("READ");
66 loader->LoadTracks("READ");
67
68 fReconstructibleTracks = 0;
69 fRecoTracks = 0;
70}
71
f29ba3e1 72//____________________________________________________________________
73AliMUONRecoCheck::AliMUONRecoCheck(const AliMUONRecoCheck& rhs)
74 : TObject(rhs)
75{
76// Protected copy constructor
77
78 AliFatal("Not implemented.");
79}
80
b8dc484b 81//_____________________________________________________________________________
82AliMUONRecoCheck::~AliMUONRecoCheck()
83{
84 fRunLoader->UnloadKinematics();
85 fRunLoader->UnloadTrackRefs();
86 fRunLoader->UnloadTracks();
29fc2c86 87 fMuonTrackRef->Delete();
b8dc484b 88 delete fMuonTrackRef;
89 delete fMUONData;
90}
91
f29ba3e1 92//________________________________________________________________________
93AliMUONRecoCheck& AliMUONRecoCheck::operator = (const AliMUONRecoCheck& rhs)
94{
95// Protected assignement operator
96
97 if (this == &rhs) return *this;
98
99 AliFatal("Not implemented.");
100
101 return *this;
102}
103
b8dc484b 104//_____________________________________________________________________________
105void AliMUONRecoCheck::MakeTrackRef()
106{
107 // Make reconstructible tracks
108
109 AliTrackReference *trackReference;
110 AliMUONTrack *muonTrack;
111 AliMUONTrackParam *trackParam;
112 AliMUONHitForRec *hitForRec;
113 Float_t x, y, z, pX, pY, pZ, pYZ;
114 Int_t track, trackSave;
115 TParticle *particle;
116 Float_t bendingSlope = 0;
117 Float_t nonBendingSlope = 0;
118 Float_t inverseBendingMomentum = 0;
119
120 TTree* treeTR = fRunLoader->TreeTR();
121 if (treeTR == NULL) return;
122
123 TBranch* branch = treeTR->GetBranch("MUON");
124 if (branch == NULL) return;
125
de876dd1 126 TClonesArray* trackRefs = 0;
b8dc484b 127 branch->SetAddress(&trackRefs);
de876dd1 128 branch->SetAutoDelete(kTRUE);
b8dc484b 129
de876dd1 130 Int_t nTrackRef = (Int_t)branch->GetEntries();
b8dc484b 131
132 track = trackSave = -999;
133 Bool_t isNewTrack;
134 Int_t iHitMin, iChamber;
135
136 trackParam = new AliMUONTrackParam();
137 hitForRec = new AliMUONHitForRec();
138 muonTrack = new AliMUONTrack();
139
140 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; iTrackRef++) {
de876dd1 141 branch->GetEntry(iTrackRef);
b8dc484b 142
143 iHitMin = 0;
144 isNewTrack = kTRUE;
dba3c809 145
29fc2c86 146 if (!trackRefs->GetEntries()) continue;
b8dc484b 147
29fc2c86 148 while (isNewTrack) {
149
b8dc484b 150 for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
151
152 trackReference = (AliTrackReference*)trackRefs->At(iHit);
153 x = trackReference->X();
154 y = trackReference->Y();
155 z = trackReference->Z();
156 pX = trackReference->Px();
157 pY = trackReference->Py();
158 pZ = trackReference->Pz();
159
160 track = trackReference->GetTrack();
161
162 if (track != trackSave && iHit != 0) {
163 iHitMin = iHit;
164 trackSave = track;
165 break;
166 }
167
168 // track parameters at hit
169 trackParam->SetBendingCoor(y);
170 trackParam->SetNonBendingCoor(x);
171 trackParam->SetZ(z);
172
173 if (TMath::Abs(pZ) > 0) {
174 bendingSlope = pY/pZ;
175 nonBendingSlope = pX/pZ;
176 }
177 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
178 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
179
180 trackParam->SetBendingSlope(bendingSlope);
181 trackParam->SetNonBendingSlope(nonBendingSlope);
182 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
183
184 hitForRec->SetBendingCoor(y);
185 hitForRec->SetNonBendingCoor(x);
186 hitForRec->SetZ(z);
187 hitForRec->SetBendingReso2(0.0);
188 hitForRec->SetNonBendingReso2(0.0);
29fc2c86 189 iChamber = AliMUONConstants::ChamberNumber(z);
b8dc484b 190 hitForRec->SetChamberNumber(iChamber);
191
192 muonTrack->AddTrackParamAtHit(trackParam);
193 muonTrack->AddHitForRecAtHit(hitForRec);
194 muonTrack->SetTrackID(track);
195
196 trackSave = track;
197 if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
198 }
199
200 // track parameters at vertex
201 particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
202 if (particle) {
203 x = particle->Vx();
204 y = particle->Vy();
205 z = particle->Vz();
206 pX = particle->Px();
207 pY = particle->Py();
208 pZ = particle->Pz();
209
210 trackParam->SetBendingCoor(y);
211 trackParam->SetNonBendingCoor(x);
212 trackParam->SetZ(z);
213 if (TMath::Abs(pZ) > 0) {
214 bendingSlope = pY/pZ;
215 nonBendingSlope = pX/pZ;
216 }
217 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
218 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
219 trackParam->SetBendingSlope(bendingSlope);
220 trackParam->SetNonBendingSlope(nonBendingSlope);
221 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
222
223 muonTrack->SetTrackParamAtVertex(trackParam);
224 }
225
226 AddMuonTrackReference(muonTrack);
227 muonTrack->ResetTrackParamAtHit();
228 muonTrack->ResetHitForRecAtHit();
229
230 } // end while isNewTrack
231
232 }
233
234 CleanMuonTrackRef();
235
236 ReconstructibleTracks();
237
238 delete muonTrack;
239 delete trackParam;
240 delete hitForRec;
b8dc484b 241
242}
243
244//____________________________________________________________________________
245TClonesArray* AliMUONRecoCheck::GetTrackReco()
246{
247 // Return TClonesArray of reconstructed tracks
248
29fc2c86 249 GetMUONData()->ResetRecTracks();
b8dc484b 250 GetMUONData()->SetTreeAddress("RT");
251 fTrackReco = GetMUONData()->RecTracks();
252 GetMUONData()->GetRecTracks();
253 fRecoTracks = fTrackReco->GetEntriesFast();
254 return fTrackReco;
255}
256//_____________________________________________________________________________
257void AliMUONRecoCheck::PrintEvent() const
258{
259 // debug facility
260
261 AliMUONTrack *track;
262 AliMUONHitForRec *hitForRec;
263 TClonesArray * hitForRecAtHit = 0;
264 Float_t xRec,yRec,zRec;
265
266 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
267
268 printf(" ******************************* \n");
269 printf(" nb of tracks %d \n",nTrackRef);
270
271 for (Int_t index = 0; index < nTrackRef; index++) {
272 track = (AliMUONTrack*)fMuonTrackRef->At(index);
273 hitForRecAtHit = track->GetHitForRecAtHit();
274 Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
275 printf(" track number %d \n",index);
276 for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
277 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
278 xRec = hitForRec->GetNonBendingCoor();
279 yRec = hitForRec->GetBendingCoor();
280 zRec = hitForRec->GetZ();
281 printf(" x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
282 }
283 }
284}
285
286//_____________________________________________________________________________
287void AliMUONRecoCheck::ResetTracks() const
288{
289 if (fMuonTrackRef) fMuonTrackRef->Clear();
290}
291//_____________________________________________________________________________
292void AliMUONRecoCheck::CleanMuonTrackRef()
293{
294 // Re-calculate hits parameters because two AliTrackReferences are recorded for
295 // each chamber (one when particle is entering + one when particle is leaving
296 // the sensitive volume)
297
298 Float_t maxGasGap = 1.; // cm
299 AliMUONTrack *track, *trackNew;
300 AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
301 AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
302 TClonesArray * hitForRecAtHit = 0;
303 TClonesArray * trackParamAtHit = 0;
304 Float_t xRec,yRec,zRec;
305 Float_t xRec1,yRec1,zRec1;
306 Float_t xRec2,yRec2,zRec2;
307 Float_t bendingSlope,nonBendingSlope,bendingMomentum;
308 Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
309 Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
310 TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
311 Int_t iHit1;
312 Int_t iChamber = 0;
313 Int_t nRec = 0;
314 Int_t nTrackHits = 0;
315
316 hitForRec = new AliMUONHitForRec();
317 trackParam = new AliMUONTrackParam();
318 trackNew = new AliMUONTrack();
319
320 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
321
322 for (Int_t index = 0; index < nTrackRef; index++) {
323 track = (AliMUONTrack*)fMuonTrackRef->At(index);
324 hitForRecAtHit = track->GetHitForRecAtHit();
325 trackParamAtHit = track->GetTrackParamAtHit();
326 trackParamAtVertex = track->GetTrackParamAtVertex();
327 nTrackHits = hitForRecAtHit->GetEntriesFast();
328 iHit1 = 0;
329 while (iHit1 < nTrackHits) {
330 hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
331 trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
332 xRec1 = hitForRec1->GetNonBendingCoor();
333 yRec1 = hitForRec1->GetBendingCoor();
334 zRec1 = hitForRec1->GetZ();
335 xRec = xRec1;
336 yRec = yRec1;
337 zRec = zRec1;
338 bendingSlope1 = trackParam1->GetBendingSlope();
339 nonBendingSlope1 = trackParam1->GetNonBendingSlope();
de876dd1 340 bendingMomentum1 = 0;
341 if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
342 bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
b8dc484b 343 bendingSlope = bendingSlope1;
344 nonBendingSlope = nonBendingSlope1;
345 bendingMomentum = bendingMomentum1;
346 nRec = 1;
347 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
348 hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
349 trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
350 xRec2 = hitForRec2->GetNonBendingCoor();
351 yRec2 = hitForRec2->GetBendingCoor();
352 zRec2 = hitForRec2->GetZ();
353 bendingSlope2 = trackParam2->GetBendingSlope();
354 nonBendingSlope2 = trackParam2->GetNonBendingSlope();
de876dd1 355 bendingMomentum2 = 0;
356 if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
357 bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
b8dc484b 358
359 if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
360 nRec++;
361 xRec += xRec2;
362 yRec += yRec2;
363 zRec += zRec2;
364 bendingSlope += bendingSlope2;
365 nonBendingSlope += nonBendingSlope2;
366 bendingMomentum += bendingMomentum2;
367 iHit1 = iHit2;
368 }
369
370 } // end iHit2
371 xRec /= (Float_t)nRec;
372 yRec /= (Float_t)nRec;
373 zRec /= (Float_t)nRec;
374 bendingSlope /= (Float_t)nRec;
375 nonBendingSlope /= (Float_t)nRec;
376 bendingMomentum /= (Float_t)nRec;
377
378 hitForRec->SetNonBendingCoor(xRec);
379 hitForRec->SetBendingCoor(yRec);
380 hitForRec->SetZ(zRec);
29fc2c86 381 iChamber = AliMUONConstants::ChamberNumber(zRec);
b8dc484b 382 hitForRec->SetChamberNumber(iChamber);
383 hitForRec->SetBendingReso2(0.0);
384 hitForRec->SetNonBendingReso2(0.0);
385 trackParam->SetNonBendingCoor(xRec);
386 trackParam->SetBendingCoor(yRec);
387 trackParam->SetZ(zRec);
388 trackParam->SetNonBendingSlope(nonBendingSlope);
389 trackParam->SetBendingSlope(bendingSlope);
de876dd1 390 if (TMath::Abs(bendingMomentum) > 0)
391 trackParam->SetInverseBendingMomentum(1./bendingMomentum);
b8dc484b 392
393 trackNew->AddHitForRecAtHit(hitForRec);
394 trackNew->AddTrackParamAtHit(trackParam);
395
396 iHit1++;
397 } // end iHit1
398
399 trackNew->SetTrackID(track->GetTrackID());
400 trackNew->SetTrackParamAtVertex(trackParamAtVertex);
401 {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}
402 trackNew->ResetHitForRecAtHit();
403 trackNew->ResetTrackParamAtHit();
404
405 } // end trackRef
406
407 fMuonTrackRef->Clear();
408 nTrackRef = newMuonTrackRef->GetEntriesFast();
409 for (Int_t index = 0; index < nTrackRef; index++) {
410 track = (AliMUONTrack*)newMuonTrackRef->At(index);
411 AddMuonTrackReference(track);
412 }
413
414
415 delete trackNew;
416 delete hitForRec;
417 delete trackParam;
29fc2c86 418 newMuonTrackRef->Delete();
b8dc484b 419 delete newMuonTrackRef;
420
421}
422
423//_____________________________________________________________________________
424void AliMUONRecoCheck::ReconstructibleTracks()
425{
426 // calculate the number of reconstructible tracks
427
428 AliMUONTrack* track;
429 TClonesArray* hitForRecAtHit = NULL;
430 AliMUONHitForRec* hitForRec;
431 Float_t zRec;
432 Int_t nTrackRef, nTrackHits;
433 Int_t isChamberInTrack[10];
434 Int_t iChamber = 0;
435 Bool_t isTrackOK = kTRUE;
436
437 fReconstructibleTracks = 0;
438
439 nTrackRef = fMuonTrackRef->GetEntriesFast();
440 for (Int_t index = 0; index < nTrackRef; index++) {
441 track = (AliMUONTrack*)fMuonTrackRef->At(index);
442 hitForRecAtHit = track->GetHitForRecAtHit();
443 nTrackHits = hitForRecAtHit->GetEntriesFast();
444 for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
445
446 for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
447 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
448 zRec = hitForRec->GetZ();
449 iChamber = hitForRec->GetChamberNumber();
450 if (iChamber < 0 || iChamber > 10) continue;
451 isChamberInTrack[iChamber] = 1;
452
453 }
454 // track is reconstructible if the particle is crossing every tracking chambers
455 isTrackOK = kTRUE;
456 for (Int_t ch = 0; ch < 10; ch++) {
457 if (!isChamberInTrack[ch]) isTrackOK = kFALSE;
29fc2c86 458 } if (isTrackOK) fReconstructibleTracks++;
b8dc484b 459 if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
460 }
461 fMuonTrackRef->Compress();
462}
463
464