Geometry builder classes moved from base to sim.
[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
cc056de1 25#include "AliMUON.h"
26#include "AliMUONRecoCheck.h"
27#include "AliMUONTrack.h"
28#include "AliMUONData.h"
29#include "AliMUONConstants.h"
b8dc484b 30
31#include "AliRun.h" // for gAlice
b8dc484b 32#include "AliLoader.h"
33#include "AliRunLoader.h"
b8dc484b 34#include "AliHeader.h"
b8dc484b 35#include "AliStack.h"
cc056de1 36#include "AliTrackReference.h"
37#include "AliLog.h"
38
39#include <TParticle.h>
b8dc484b 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;
30d3ea8a 136 Int_t iHitMin, iChamber, detElemId;
b8dc484b 137
138 trackParam = new AliMUONTrackParam();
139 hitForRec = new AliMUONHitForRec();
140 muonTrack = new AliMUONTrack();
141
cc056de1 142 Int_t max = fRunLoader->GetHeader()->Stack()->GetNtrack();
b8dc484b 143 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; iTrackRef++) {
cc056de1 144
de876dd1 145 branch->GetEntry(iTrackRef);
b8dc484b 146
147 iHitMin = 0;
148 isNewTrack = kTRUE;
dba3c809 149
29fc2c86 150 if (!trackRefs->GetEntries()) continue;
b8dc484b 151
29fc2c86 152 while (isNewTrack) {
b8dc484b 153 for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
154
155 trackReference = (AliTrackReference*)trackRefs->At(iHit);
156 x = trackReference->X();
157 y = trackReference->Y();
158 z = trackReference->Z();
159 pX = trackReference->Px();
160 pY = trackReference->Py();
161 pZ = trackReference->Pz();
162
163 track = trackReference->GetTrack();
164
cc056de1 165 if(track >= max){
166 AliWarningStream()
167 << "Track ID " << track
168 << " larger than max number of particles " << max << endl;
169 isNewTrack = kFALSE;
170 break;
171 }
b8dc484b 172 if (track != trackSave && iHit != 0) {
173 iHitMin = iHit;
174 trackSave = track;
175 break;
176 }
177
178 // track parameters at hit
179 trackParam->SetBendingCoor(y);
180 trackParam->SetNonBendingCoor(x);
181 trackParam->SetZ(z);
cc056de1 182
b8dc484b 183 if (TMath::Abs(pZ) > 0) {
184 bendingSlope = pY/pZ;
185 nonBendingSlope = pX/pZ;
186 }
187 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
188 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
189
190 trackParam->SetBendingSlope(bendingSlope);
191 trackParam->SetNonBendingSlope(nonBendingSlope);
192 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
193
194 hitForRec->SetBendingCoor(y);
195 hitForRec->SetNonBendingCoor(x);
196 hitForRec->SetZ(z);
197 hitForRec->SetBendingReso2(0.0);
30d3ea8a 198 hitForRec->SetNonBendingReso2(0.0);
199 detElemId = hitForRec->GetDetElemId();
200 if (detElemId) iChamber = detElemId / 100 - 1;
201 else iChamber = AliMUONConstants::ChamberNumber(z);
b8dc484b 202 hitForRec->SetChamberNumber(iChamber);
203
204 muonTrack->AddTrackParamAtHit(trackParam);
205 muonTrack->AddHitForRecAtHit(hitForRec);
206 muonTrack->SetTrackID(track);
cc056de1 207
b8dc484b 208 trackSave = track;
209 if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
210 }
cc056de1 211
b8dc484b 212 // track parameters at vertex
213 particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
cc056de1 214
b8dc484b 215 if (particle) {
cc056de1 216
b8dc484b 217 x = particle->Vx();
218 y = particle->Vy();
219 z = particle->Vz();
220 pX = particle->Px();
221 pY = particle->Py();
222 pZ = particle->Pz();
223
224 trackParam->SetBendingCoor(y);
225 trackParam->SetNonBendingCoor(x);
226 trackParam->SetZ(z);
227 if (TMath::Abs(pZ) > 0) {
228 bendingSlope = pY/pZ;
229 nonBendingSlope = pX/pZ;
230 }
231 pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
232 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
233 trackParam->SetBendingSlope(bendingSlope);
234 trackParam->SetNonBendingSlope(nonBendingSlope);
235 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
236
237 muonTrack->SetTrackParamAtVertex(trackParam);
238 }
cc056de1 239
b8dc484b 240 AddMuonTrackReference(muonTrack);
241 muonTrack->ResetTrackParamAtHit();
242 muonTrack->ResetHitForRecAtHit();
243
244 } // end while isNewTrack
b8dc484b 245 }
246
247 CleanMuonTrackRef();
248
249 ReconstructibleTracks();
250
251 delete muonTrack;
252 delete trackParam;
253 delete hitForRec;
b8dc484b 254
255}
256
257//____________________________________________________________________________
258TClonesArray* AliMUONRecoCheck::GetTrackReco()
259{
d19b6003 260/// Return TClonesArray of reconstructed tracks
b8dc484b 261
29fc2c86 262 GetMUONData()->ResetRecTracks();
b8dc484b 263 GetMUONData()->SetTreeAddress("RT");
264 fTrackReco = GetMUONData()->RecTracks();
265 GetMUONData()->GetRecTracks();
266 fRecoTracks = fTrackReco->GetEntriesFast();
267 return fTrackReco;
268}
269//_____________________________________________________________________________
270void AliMUONRecoCheck::PrintEvent() const
271{
d19b6003 272/// Debug facility
b8dc484b 273
274 AliMUONTrack *track;
275 AliMUONHitForRec *hitForRec;
276 TClonesArray * hitForRecAtHit = 0;
277 Float_t xRec,yRec,zRec;
278
279 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
280
281 printf(" ******************************* \n");
282 printf(" nb of tracks %d \n",nTrackRef);
283
284 for (Int_t index = 0; index < nTrackRef; index++) {
285 track = (AliMUONTrack*)fMuonTrackRef->At(index);
286 hitForRecAtHit = track->GetHitForRecAtHit();
287 Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
288 printf(" track number %d \n",index);
289 for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
290 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
291 xRec = hitForRec->GetNonBendingCoor();
292 yRec = hitForRec->GetBendingCoor();
293 zRec = hitForRec->GetZ();
294 printf(" x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
295 }
296 }
297}
298
299//_____________________________________________________________________________
300void AliMUONRecoCheck::ResetTracks() const
301{
d19b6003 302/// Reset tracks
303
b8dc484b 304 if (fMuonTrackRef) fMuonTrackRef->Clear();
305}
306//_____________________________________________________________________________
307void AliMUONRecoCheck::CleanMuonTrackRef()
308{
d19b6003 309/// Re-calculate hits parameters because two AliTrackReferences are recorded for
310/// each chamber (one when particle is entering + one when particle is leaving
311/// the sensitive volume)
b8dc484b 312
313 Float_t maxGasGap = 1.; // cm
314 AliMUONTrack *track, *trackNew;
315 AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
316 AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
317 TClonesArray * hitForRecAtHit = 0;
318 TClonesArray * trackParamAtHit = 0;
319 Float_t xRec,yRec,zRec;
320 Float_t xRec1,yRec1,zRec1;
321 Float_t xRec2,yRec2,zRec2;
322 Float_t bendingSlope,nonBendingSlope,bendingMomentum;
323 Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
324 Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
325 TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
326 Int_t iHit1;
30d3ea8a 327 Int_t iChamber = 0, detElemId = 0;
b8dc484b 328 Int_t nRec = 0;
329 Int_t nTrackHits = 0;
330
331 hitForRec = new AliMUONHitForRec();
332 trackParam = new AliMUONTrackParam();
333 trackNew = new AliMUONTrack();
334
335 Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
b8dc484b 336 for (Int_t index = 0; index < nTrackRef; index++) {
337 track = (AliMUONTrack*)fMuonTrackRef->At(index);
338 hitForRecAtHit = track->GetHitForRecAtHit();
339 trackParamAtHit = track->GetTrackParamAtHit();
340 trackParamAtVertex = track->GetTrackParamAtVertex();
341 nTrackHits = hitForRecAtHit->GetEntriesFast();
342 iHit1 = 0;
343 while (iHit1 < nTrackHits) {
344 hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
345 trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
346 xRec1 = hitForRec1->GetNonBendingCoor();
347 yRec1 = hitForRec1->GetBendingCoor();
348 zRec1 = hitForRec1->GetZ();
349 xRec = xRec1;
350 yRec = yRec1;
351 zRec = zRec1;
352 bendingSlope1 = trackParam1->GetBendingSlope();
353 nonBendingSlope1 = trackParam1->GetNonBendingSlope();
de876dd1 354 bendingMomentum1 = 0;
355 if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
356 bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
b8dc484b 357 bendingSlope = bendingSlope1;
358 nonBendingSlope = nonBendingSlope1;
359 bendingMomentum = bendingMomentum1;
360 nRec = 1;
361 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
362 hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
363 trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
364 xRec2 = hitForRec2->GetNonBendingCoor();
365 yRec2 = hitForRec2->GetBendingCoor();
366 zRec2 = hitForRec2->GetZ();
367 bendingSlope2 = trackParam2->GetBendingSlope();
368 nonBendingSlope2 = trackParam2->GetNonBendingSlope();
de876dd1 369 bendingMomentum2 = 0;
370 if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
371 bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
b8dc484b 372
373 if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
cc056de1 374
b8dc484b 375 nRec++;
376 xRec += xRec2;
377 yRec += yRec2;
378 zRec += zRec2;
379 bendingSlope += bendingSlope2;
380 nonBendingSlope += nonBendingSlope2;
381 bendingMomentum += bendingMomentum2;
382 iHit1 = iHit2;
383 }
384
385 } // end iHit2
386 xRec /= (Float_t)nRec;
387 yRec /= (Float_t)nRec;
388 zRec /= (Float_t)nRec;
389 bendingSlope /= (Float_t)nRec;
390 nonBendingSlope /= (Float_t)nRec;
391 bendingMomentum /= (Float_t)nRec;
392
393 hitForRec->SetNonBendingCoor(xRec);
394 hitForRec->SetBendingCoor(yRec);
395 hitForRec->SetZ(zRec);
30d3ea8a 396 detElemId = hitForRec->GetDetElemId();
397 if (detElemId) iChamber = detElemId / 100 - 1;
398 else iChamber = AliMUONConstants::ChamberNumber(zRec);
b8dc484b 399 hitForRec->SetChamberNumber(iChamber);
400 hitForRec->SetBendingReso2(0.0);
401 hitForRec->SetNonBendingReso2(0.0);
402 trackParam->SetNonBendingCoor(xRec);
403 trackParam->SetBendingCoor(yRec);
404 trackParam->SetZ(zRec);
405 trackParam->SetNonBendingSlope(nonBendingSlope);
406 trackParam->SetBendingSlope(bendingSlope);
de876dd1 407 if (TMath::Abs(bendingMomentum) > 0)
408 trackParam->SetInverseBendingMomentum(1./bendingMomentum);
b8dc484b 409
410 trackNew->AddHitForRecAtHit(hitForRec);
411 trackNew->AddTrackParamAtHit(trackParam);
412
413 iHit1++;
414 } // end iHit1
415
416 trackNew->SetTrackID(track->GetTrackID());
417 trackNew->SetTrackParamAtVertex(trackParamAtVertex);
418 {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}
419 trackNew->ResetHitForRecAtHit();
420 trackNew->ResetTrackParamAtHit();
421
422 } // end trackRef
423
424 fMuonTrackRef->Clear();
425 nTrackRef = newMuonTrackRef->GetEntriesFast();
426 for (Int_t index = 0; index < nTrackRef; index++) {
427 track = (AliMUONTrack*)newMuonTrackRef->At(index);
428 AddMuonTrackReference(track);
429 }
430
431
432 delete trackNew;
433 delete hitForRec;
434 delete trackParam;
29fc2c86 435 newMuonTrackRef->Delete();
b8dc484b 436 delete newMuonTrackRef;
437
438}
439
440//_____________________________________________________________________________
441void AliMUONRecoCheck::ReconstructibleTracks()
442{
d19b6003 443/// Calculate the number of reconstructible tracks
b8dc484b 444
445 AliMUONTrack* track;
446 TClonesArray* hitForRecAtHit = NULL;
447 AliMUONHitForRec* hitForRec;
448 Float_t zRec;
449 Int_t nTrackRef, nTrackHits;
450 Int_t isChamberInTrack[10];
451 Int_t iChamber = 0;
452 Bool_t isTrackOK = kTRUE;
453
454 fReconstructibleTracks = 0;
455
456 nTrackRef = fMuonTrackRef->GetEntriesFast();
457 for (Int_t index = 0; index < nTrackRef; index++) {
458 track = (AliMUONTrack*)fMuonTrackRef->At(index);
459 hitForRecAtHit = track->GetHitForRecAtHit();
460 nTrackHits = hitForRecAtHit->GetEntriesFast();
461 for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
462
463 for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
464 hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
465 zRec = hitForRec->GetZ();
466 iChamber = hitForRec->GetChamberNumber();
467 if (iChamber < 0 || iChamber > 10) continue;
468 isChamberInTrack[iChamber] = 1;
cc056de1 469 }
470 // track is reconstructible if the particle is depositing a hit
471 // in the following chamber combinations:
472
b8dc484b 473 isTrackOK = kTRUE;
cc056de1 474 if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
475 if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
476 if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
477 Int_t nHitsInLastStations=0;
478 for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
479 if (isChamberInTrack[ch]) nHitsInLastStations++;
480 if(nHitsInLastStations < 3) isTrackOK = kFALSE;
481
482 if (isTrackOK) fReconstructibleTracks++;
b8dc484b 483 if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
484 }
485 fMuonTrackRef->Compress();
486}
487