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