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