Added DetElemId() method and changed the Print() (Laurent)
[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
d19b6003 16/* $Id$ */
17
5398f946 18/// \class AliMUONRecoCheck
78649106 19/// Utility class to check reconstruction
5398f946 20/// Reconstructed tracks are compared to reference tracks.
21/// The reference tracks are built from AliTrackReference for the
22/// hit in chamber (0..9) and from kinematics for the vertex parameters.
b8dc484b 23
cc056de1 24#include "AliMUON.h"
25#include "AliMUONRecoCheck.h"
26#include "AliMUONTrack.h"
6b092dfc 27#include "AliMUONRecData.h"
cc056de1 28#include "AliMUONConstants.h"
b8dc484b 29
b8dc484b 30#include "AliLoader.h"
31#include "AliRunLoader.h"
b8dc484b 32#include "AliHeader.h"
b8dc484b 33#include "AliStack.h"
cc056de1 34#include "AliTrackReference.h"
35#include "AliLog.h"
36
37#include <TParticle.h>
5b1afea0 38#include <TParticlePDG.h>
b8dc484b 39
5398f946 40/// \cond CLASSIMP
b8dc484b 41ClassImp(AliMUONRecoCheck)
5398f946 42/// \endcond
b8dc484b 43
44//_____________________________________________________________________________
6b092dfc 45 AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader, Char_t *chLoaderSim)
068e4c36 46 : TObject(),
47 fRunLoader(0x0),
48 fMUONData(0x0),
6b092dfc 49 fRunLoaderSim(0x0),
50 fMUONDataSim(0x0),
068e4c36 51 fMuonTrackRef(0x0),
52 fTrackReco(0x0),
53 fReconstructibleTracks(0),
54 fRecoTracks(0),
55 fIsLoadConstructor(kTRUE)
56{
57/// Constructor using "galice.root",
58/// takes care of loading/unloading Kinematics, TrackRefs and Tracks
59
60 fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
61
62 // run loader
6b092dfc 63 fRunLoader = AliRunLoader::Open(chLoader,"MUONFolder","READ");
068e4c36 64 if (!fRunLoader) {
65 AliError(Form("no run loader found " ));
66 return;
67 }
68
69 // initialize loader
70 AliLoader *loader = fRunLoader->GetLoader("MUONLoader");
71
72 // container
6b092dfc 73 fMUONData = new AliMUONRecData(loader,"MUON","MUON");
068e4c36 74 if (!fMUONData) {
75 AliError(Form("no MUONData found " ));
76 return;
77 }
78
6b092dfc 79 // run loader
80 fRunLoaderSim = AliRunLoader::Open(chLoaderSim,"MUONFolderSim","READ");
81 if (!fRunLoaderSim) {
82 AliError(Form("no run sim loader found " ));
83 return;
84 }
85
86 // initialize loader
87 AliLoader *loaderSim = fRunLoaderSim->GetLoader("MUONLoader");
88
89 // container
90 fMUONDataSim = new AliMUONSimData(loaderSim,"MUON","MUON");
91 if (!fMUONDataSim) {
92 AliError(Form("no MUONDataSim found " ));
93 return;
94 }
95
96 fRunLoaderSim->LoadKinematics("READ");
cad8ba75 97 fRunLoaderSim->LoadTrackRefs("READ");
068e4c36 98 loader->LoadTracks("READ");
99
100}
101
102//_____________________________________________________________________________
6b092dfc 103 AliMUONRecoCheck::AliMUONRecoCheck(AliRunLoader *runloader, AliMUONRecData *muondata,
104 AliRunLoader *runloaderSim, AliMUONSimData *muondataSim)
54d7ba50 105 : TObject(),
106 fRunLoader(0x0),
107 fMUONData(0x0),
6b092dfc 108 fRunLoaderSim(0x0),
109 fMUONDataSim(0x0),
54d7ba50 110 fMuonTrackRef(0x0),
111 fTrackReco(0x0),
112 fReconstructibleTracks(0),
068e4c36 113 fRecoTracks(0),
114 fIsLoadConstructor(kFALSE)
b8dc484b 115{
068e4c36 116/// Constructor from AliRunLoader and AliMUONData,
117/// does not load/unload Kinematics, TrackRefs and Tracks internally
118/// it should be done in the execution program (e.g. see MUONRecoCheck.C)
5398f946 119
b8dc484b 120 fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
121
45a05f07 122 // run loader
123 fRunLoader = runloader;
b8dc484b 124 if (!fRunLoader) {
45a05f07 125 AliError(Form("no run loader found " ));
b8dc484b 126 return;
127 }
b8dc484b 128
45a05f07 129 // container
130 fMUONData = muondata;
131 if (!fMUONData) {
132 AliError(Form("no MUONData found " ));
133 return;
134 }
b8dc484b 135
6b092dfc 136 // run loader
137 fRunLoaderSim = runloaderSim;
138 if (!fRunLoaderSim) {
139 AliError(Form("no run sim loader found " ));
140 return;
141 }
142
143 // container
144 fMUONDataSim = muondataSim;
145 if (!fMUONDataSim) {
146 AliError(Form("no MUONDataSim found " ));
147 return;
148 }
149
b8dc484b 150}
151
45a05f07 152
b8dc484b 153//_____________________________________________________________________________
154AliMUONRecoCheck::~AliMUONRecoCheck()
155{
d19b6003 156/// Destructor
b8dc484b 157 delete fMuonTrackRef;
068e4c36 158 if(fIsLoadConstructor){
6b092dfc 159 fRunLoaderSim->UnloadKinematics();
cad8ba75 160 fRunLoaderSim->UnloadTrackRefs();
068e4c36 161 fRunLoader->UnloadTracks();
162 delete fMUONData;
fac70e25 163 delete fRunLoader;
6b092dfc 164 delete fMUONDataSim;
165 delete fRunLoaderSim;
068e4c36 166 }
b8dc484b 167}
168
169//_____________________________________________________________________________
170void AliMUONRecoCheck::MakeTrackRef()
171{
d19b6003 172/// Make reconstructible tracks
b8dc484b 173
174 AliTrackReference *trackReference;
175 AliMUONTrack *muonTrack;
176 AliMUONTrackParam *trackParam;
177 AliMUONHitForRec *hitForRec;
178 Float_t x, y, z, pX, pY, pZ, pYZ;
179 Int_t track, trackSave;
180 TParticle *particle;
181 Float_t bendingSlope = 0;
182 Float_t nonBendingSlope = 0;
183 Float_t inverseBendingMomentum = 0;
184
cad8ba75 185 TTree* treeTR = fRunLoaderSim->TreeTR();
b8dc484b 186 if (treeTR == NULL) return;
187
188 TBranch* branch = treeTR->GetBranch("MUON");
189 if (branch == NULL) return;
190
de876dd1 191 TClonesArray* trackRefs = 0;
b8dc484b 192 branch->SetAddress(&trackRefs);
de876dd1 193 branch->SetAutoDelete(kTRUE);
b8dc484b 194
de876dd1 195 Int_t nTrackRef = (Int_t)branch->GetEntries();
b8dc484b 196
197 track = trackSave = -999;
198 Bool_t isNewTrack;
30d3ea8a 199 Int_t iHitMin, iChamber, detElemId;
b8dc484b 200
201 trackParam = new AliMUONTrackParam();
202 hitForRec = new AliMUONHitForRec();
b8dc484b 203
5b1afea0 204 Int_t charge;
205 TParticlePDG *ppdg;
206
6b092dfc 207 Int_t max = fRunLoaderSim->GetHeader()->Stack()->GetNtrack();
b8dc484b 208 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; iTrackRef++) {
cc056de1 209
de876dd1 210 branch->GetEntry(iTrackRef);
b8dc484b 211
212 iHitMin = 0;
213 isNewTrack = kTRUE;
dba3c809 214
29fc2c86 215 if (!trackRefs->GetEntries()) continue;
b8dc484b 216
29fc2c86 217 while (isNewTrack) {
068e4c36 218
219 muonTrack = new AliMUONTrack();
220
b8dc484b 221 for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
222
223 trackReference = (AliTrackReference*)trackRefs->At(iHit);
224 x = trackReference->X();
225 y = trackReference->Y();
226 z = trackReference->Z();
227 pX = trackReference->Px();
228 pY = trackReference->Py();
229 pZ = trackReference->Pz();
230
231 track = trackReference->GetTrack();
232
cc056de1 233 if(track >= max){
234 AliWarningStream()
235 << "Track ID " << track
236 << " larger than max number of particles " << max << endl;
237 isNewTrack = kFALSE;
238 break;
239 }
b8dc484b 240 if (track != trackSave && iHit != 0) {
241 iHitMin = iHit;
242 trackSave = track;
243 break;
244 }
245
246 // track parameters at hit
247 trackParam->SetBendingCoor(y);
248 trackParam->SetNonBendingCoor(x);
249 trackParam->SetZ(z);
cc056de1 250
b8dc484b 251 if (TMath::Abs(pZ) > 0) {
252 bendingSlope = pY/pZ;
253 nonBendingSlope = pX/pZ;
254 }
255 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
256 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
257
258 trackParam->SetBendingSlope(bendingSlope);
259 trackParam->SetNonBendingSlope(nonBendingSlope);
260 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
261
262 hitForRec->SetBendingCoor(y);
263 hitForRec->SetNonBendingCoor(x);
264 hitForRec->SetZ(z);
265 hitForRec->SetBendingReso2(0.0);
30d3ea8a 266 hitForRec->SetNonBendingReso2(0.0);
267 detElemId = hitForRec->GetDetElemId();
268 if (detElemId) iChamber = detElemId / 100 - 1;
269 else iChamber = AliMUONConstants::ChamberNumber(z);
b8dc484b 270 hitForRec->SetChamberNumber(iChamber);
271
068e4c36 272 muonTrack->AddTrackParamAtHit(trackParam,0);
b8dc484b 273 muonTrack->AddHitForRecAtHit(hitForRec);
274 muonTrack->SetTrackID(track);
cc056de1 275
b8dc484b 276 trackSave = track;
277 if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
278 }
cc056de1 279
b8dc484b 280 // track parameters at vertex
6b092dfc 281 particle = fRunLoaderSim->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
cc056de1 282
b8dc484b 283 if (particle) {
cc056de1 284
b8dc484b 285 x = particle->Vx();
286 y = particle->Vy();
287 z = particle->Vz();
288 pX = particle->Px();
289 pY = particle->Py();
290 pZ = particle->Pz();
291
292 trackParam->SetBendingCoor(y);
c61c826b 293 trackParam->SetNonBendingCoor(x);
b8dc484b 294 trackParam->SetZ(z);
295 if (TMath::Abs(pZ) > 0) {
296 bendingSlope = pY/pZ;
297 nonBendingSlope = pX/pZ;
298 }
5b1afea0 299
b8dc484b 300 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
301 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
5b1afea0 302
303 ppdg = particle->GetPDG(1);
304 charge = (Int_t)(ppdg->Charge()/3.0);
305 inverseBendingMomentum *= charge;
306
b8dc484b 307 trackParam->SetBendingSlope(bendingSlope);
308 trackParam->SetNonBendingSlope(nonBendingSlope);
309 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
310
311 muonTrack->SetTrackParamAtVertex(trackParam);
312 }
cc056de1 313
b8dc484b 314 AddMuonTrackReference(muonTrack);
068e4c36 315 delete muonTrack;
316 muonTrack = NULL;
b8dc484b 317
318 } // end while isNewTrack
b8dc484b 319 }
2782bbd7 320 delete trackRefs;
b8dc484b 321 CleanMuonTrackRef();
322
323 ReconstructibleTracks();
324
b8dc484b 325 delete trackParam;
326 delete hitForRec;
b8dc484b 327
328}
329
330//____________________________________________________________________________
331TClonesArray* AliMUONRecoCheck::GetTrackReco()
332{
d19b6003 333/// Return TClonesArray of reconstructed tracks
b8dc484b 334
45a05f07 335 fMUONData->ResetRecTracks();
336 fMUONData->SetTreeAddress("RT");
337 fTrackReco = fMUONData->RecTracks();
338 fMUONData->GetRecTracks();
b8dc484b 339 fRecoTracks = fTrackReco->GetEntriesFast();
340 return fTrackReco;
341}
342//_____________________________________________________________________________
343void AliMUONRecoCheck::PrintEvent() const
344{
d19b6003 345/// Debug facility
b8dc484b 346
347 AliMUONTrack *track;
348 AliMUONHitForRec *hitForRec;
349 TClonesArray * hitForRecAtHit = 0;
350 Float_t xRec,yRec,zRec;
351
352 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
353
354 printf(" ******************************* \n");
355 printf(" nb of tracks %d \n",nTrackRef);
356
357 for (Int_t index = 0; index < nTrackRef; index++) {
358 track = (AliMUONTrack*)fMuonTrackRef->At(index);
359 hitForRecAtHit = track->GetHitForRecAtHit();
360 Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
361 printf(" track number %d \n",index);
362 for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
363 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
364 xRec = hitForRec->GetNonBendingCoor();
365 yRec = hitForRec->GetBendingCoor();
366 zRec = hitForRec->GetZ();
367 printf(" x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
368 }
369 }
370}
371
372//_____________________________________________________________________________
373void AliMUONRecoCheck::ResetTracks() const
374{
d19b6003 375/// Reset tracks
376
2782bbd7 377 if (fMuonTrackRef) fMuonTrackRef->Delete();
b8dc484b 378}
379//_____________________________________________________________________________
380void AliMUONRecoCheck::CleanMuonTrackRef()
381{
d19b6003 382/// Re-calculate hits parameters because two AliTrackReferences are recorded for
383/// each chamber (one when particle is entering + one when particle is leaving
384/// the sensitive volume)
b8dc484b 385
386 Float_t maxGasGap = 1.; // cm
387 AliMUONTrack *track, *trackNew;
388 AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
389 AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
390 TClonesArray * hitForRecAtHit = 0;
391 TClonesArray * trackParamAtHit = 0;
392 Float_t xRec,yRec,zRec;
393 Float_t xRec1,yRec1,zRec1;
394 Float_t xRec2,yRec2,zRec2;
395 Float_t bendingSlope,nonBendingSlope,bendingMomentum;
396 Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
397 Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
398 TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
399 Int_t iHit1;
30d3ea8a 400 Int_t iChamber = 0, detElemId = 0;
b8dc484b 401 Int_t nRec = 0;
402 Int_t nTrackHits = 0;
403
404 hitForRec = new AliMUONHitForRec();
405 trackParam = new AliMUONTrackParam();
b8dc484b 406
407 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
b8dc484b 408 for (Int_t index = 0; index < nTrackRef; index++) {
409 track = (AliMUONTrack*)fMuonTrackRef->At(index);
410 hitForRecAtHit = track->GetHitForRecAtHit();
411 trackParamAtHit = track->GetTrackParamAtHit();
412 trackParamAtVertex = track->GetTrackParamAtVertex();
413 nTrackHits = hitForRecAtHit->GetEntriesFast();
068e4c36 414 trackNew = new AliMUONTrack();
b8dc484b 415 iHit1 = 0;
416 while (iHit1 < nTrackHits) {
417 hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
418 trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
419 xRec1 = hitForRec1->GetNonBendingCoor();
420 yRec1 = hitForRec1->GetBendingCoor();
421 zRec1 = hitForRec1->GetZ();
422 xRec = xRec1;
423 yRec = yRec1;
424 zRec = zRec1;
425 bendingSlope1 = trackParam1->GetBendingSlope();
426 nonBendingSlope1 = trackParam1->GetNonBendingSlope();
de876dd1 427 bendingMomentum1 = 0;
428 if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
429 bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
b8dc484b 430 bendingSlope = bendingSlope1;
431 nonBendingSlope = nonBendingSlope1;
432 bendingMomentum = bendingMomentum1;
433 nRec = 1;
434 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
435 hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
436 trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
437 xRec2 = hitForRec2->GetNonBendingCoor();
438 yRec2 = hitForRec2->GetBendingCoor();
439 zRec2 = hitForRec2->GetZ();
440 bendingSlope2 = trackParam2->GetBendingSlope();
441 nonBendingSlope2 = trackParam2->GetNonBendingSlope();
de876dd1 442 bendingMomentum2 = 0;
443 if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
444 bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
b8dc484b 445
446 if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
cc056de1 447
b8dc484b 448 nRec++;
449 xRec += xRec2;
450 yRec += yRec2;
451 zRec += zRec2;
452 bendingSlope += bendingSlope2;
453 nonBendingSlope += nonBendingSlope2;
454 bendingMomentum += bendingMomentum2;
455 iHit1 = iHit2;
456 }
457
458 } // end iHit2
459 xRec /= (Float_t)nRec;
460 yRec /= (Float_t)nRec;
461 zRec /= (Float_t)nRec;
462 bendingSlope /= (Float_t)nRec;
463 nonBendingSlope /= (Float_t)nRec;
464 bendingMomentum /= (Float_t)nRec;
465
466 hitForRec->SetNonBendingCoor(xRec);
467 hitForRec->SetBendingCoor(yRec);
468 hitForRec->SetZ(zRec);
30d3ea8a 469 detElemId = hitForRec->GetDetElemId();
470 if (detElemId) iChamber = detElemId / 100 - 1;
471 else iChamber = AliMUONConstants::ChamberNumber(zRec);
b8dc484b 472 hitForRec->SetChamberNumber(iChamber);
473 hitForRec->SetBendingReso2(0.0);
474 hitForRec->SetNonBendingReso2(0.0);
475 trackParam->SetNonBendingCoor(xRec);
476 trackParam->SetBendingCoor(yRec);
477 trackParam->SetZ(zRec);
478 trackParam->SetNonBendingSlope(nonBendingSlope);
479 trackParam->SetBendingSlope(bendingSlope);
de876dd1 480 if (TMath::Abs(bendingMomentum) > 0)
481 trackParam->SetInverseBendingMomentum(1./bendingMomentum);
b8dc484b 482
483 trackNew->AddHitForRecAtHit(hitForRec);
068e4c36 484 trackNew->AddTrackParamAtHit(trackParam,0);
b8dc484b 485
486 iHit1++;
487 } // end iHit1
488
489 trackNew->SetTrackID(track->GetTrackID());
490 trackNew->SetTrackParamAtVertex(trackParamAtVertex);
068e4c36 491 {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}
492 delete trackNew;
493 trackNew = NULL;
b8dc484b 494
495 } // end trackRef
496
2782bbd7 497 fMuonTrackRef->Delete();
b8dc484b 498 nTrackRef = newMuonTrackRef->GetEntriesFast();
499 for (Int_t index = 0; index < nTrackRef; index++) {
500 track = (AliMUONTrack*)newMuonTrackRef->At(index);
501 AddMuonTrackReference(track);
502 }
503
b8dc484b 504 delete hitForRec;
505 delete trackParam;
29fc2c86 506 newMuonTrackRef->Delete();
b8dc484b 507 delete newMuonTrackRef;
508
509}
510
511//_____________________________________________________________________________
512void AliMUONRecoCheck::ReconstructibleTracks()
513{
d19b6003 514/// Calculate the number of reconstructible tracks
b8dc484b 515
516 AliMUONTrack* track;
517 TClonesArray* hitForRecAtHit = NULL;
518 AliMUONHitForRec* hitForRec;
519 Float_t zRec;
520 Int_t nTrackRef, nTrackHits;
521 Int_t isChamberInTrack[10];
522 Int_t iChamber = 0;
523 Bool_t isTrackOK = kTRUE;
524
525 fReconstructibleTracks = 0;
526
527 nTrackRef = fMuonTrackRef->GetEntriesFast();
528 for (Int_t index = 0; index < nTrackRef; index++) {
529 track = (AliMUONTrack*)fMuonTrackRef->At(index);
530 hitForRecAtHit = track->GetHitForRecAtHit();
531 nTrackHits = hitForRecAtHit->GetEntriesFast();
532 for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
533
534 for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
535 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
536 zRec = hitForRec->GetZ();
537 iChamber = hitForRec->GetChamberNumber();
538 if (iChamber < 0 || iChamber > 10) continue;
539 isChamberInTrack[iChamber] = 1;
cc056de1 540 }
541 // track is reconstructible if the particle is depositing a hit
542 // in the following chamber combinations:
543
b8dc484b 544 isTrackOK = kTRUE;
cc056de1 545 if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
546 if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
547 if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
548 Int_t nHitsInLastStations=0;
549 for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
550 if (isChamberInTrack[ch]) nHitsInLastStations++;
551 if(nHitsInLastStations < 3) isTrackOK = kFALSE;
552
553 if (isTrackOK) fReconstructibleTracks++;
b8dc484b 554 if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
555 }
556 fMuonTrackRef->Compress();
557}
558