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