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