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