- AliMUONRecoParam.cxx:
[u/mrichter/AliRoot.git] / MUON / AliMUONRecoCheck.cxx
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
16 /* $Id$ */
17
18 //-----------------------------------------------------------------------------
19 /// \class AliMUONRecoCheck
20 /// Utility class to check reconstruction
21 /// Reconstructed tracks are compared to reference tracks. 
22 /// The reference tracks are built from AliTrackReference for the
23 /// hit in chamber (0..9) and from kinematics for the vertex parameters.     
24 //-----------------------------------------------------------------------------
25
26 // 13 Nov 2007:
27 // Added a method to create a list of reconstructed AliMUONTrack objects from
28 // ESD data. This is necessary since the track objects that are actually created
29 // during offline reconstruction are no longer stored to disk.
30 //  - Artur Szostak <artursz@iafrica.com>
31 // 25 Jan 2008:
32 // Use the new ESDInterface to create MUON objects from ESD data
33 // - Philippe Pillot
34
35 #include "AliMUONRecoCheck.h"
36 #include "AliMUONTrack.h"
37 #include "AliMUONVTrackStore.h"
38 #include "AliMUONVCluster.h"
39 #include "AliMUONVClusterStore.h"
40 #include "AliMUONConstants.h"
41 #include "AliMUONESDInterface.h"
42 #include "AliMCEventHandler.h"
43 #include "AliMCEvent.h"
44 #include "AliStack.h"
45 #include "AliTrackReference.h"
46 #include "AliLog.h"
47 #include "AliESDEvent.h"
48 #include "AliESDMuonTrack.h"
49
50 #include <TFile.h>
51 #include <TTree.h>
52 #include <TParticle.h>
53 #include <TParticlePDG.h>
54 #include <Riostream.h>
55
56 /// \cond CLASSIMP
57 ClassImp(AliMUONRecoCheck)
58 /// \endcond
59
60 //_____________________________________________________________________________
61 AliMUONRecoCheck::AliMUONRecoCheck(Char_t *esdFileName, Char_t *pathSim)
62 : TObject(),
63 fMCEventHandler(new AliMCEventHandler()),
64 fESDEvent(new AliESDEvent()),
65 fESDTree (0x0),
66 fESDFile (0x0),
67 fCurrentEvent(0),
68 fTrackRefStore(0x0),
69 fRecoTrackRefStore(0x0),
70 fRecoTrackStore(0x0)
71 {
72   /// Normal ctor
73   
74   // TrackRefs and Particules
75   fMCEventHandler->SetInputPath(pathSim);
76   fMCEventHandler->InitIO("");
77   
78   // ESD MUON Tracks
79   fESDFile = TFile::Open(esdFileName); // open the file
80   if (!fESDFile || !fESDFile->IsOpen()) {
81     AliError(Form("opening ESD file %s failed", esdFileName));
82     fESDFile = 0x0;
83     return;
84   }
85   fESDTree = (TTree*) fESDFile->Get("esdTree"); // get the tree
86   if (!fESDTree) {
87     AliError("no ESD tree found");
88     fESDFile->Close();
89     fESDFile = 0x0;
90     return;
91   }
92   fESDEvent->ReadFromTree(fESDTree); // link fESDEvent to the tree
93 }
94
95 //_____________________________________________________________________________
96 AliMUONRecoCheck::~AliMUONRecoCheck()
97 {
98   /// Destructor
99   delete fMCEventHandler;
100   delete fESDEvent;
101   if (fESDFile) fESDFile->Close();
102   ResetStores();
103 }
104
105 //_____________________________________________________________________________
106 void AliMUONRecoCheck::ResetStores()
107 {
108   /// Deletes all the store objects that have been created and resets the pointers to 0x0
109   delete fTrackRefStore;      fTrackRefStore = 0x0;
110   delete fRecoTrackRefStore;  fRecoTrackRefStore = 0x0;
111   delete fRecoTrackStore;     fRecoTrackStore = 0x0;
112 }
113
114 //_____________________________________________________________________________
115 Int_t AliMUONRecoCheck::NumberOfEvents() const
116 {
117   /// Return the number of events
118   if (fESDTree) return fESDTree->GetEntries();
119   return 0;
120 }
121
122 //_____________________________________________________________________________
123 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event)
124 {
125   /// Return a track store containing the reconstructed tracks (converted into 
126   /// MUONTrack objects) for a given event
127   if (event != fCurrentEvent) {
128     ResetStores();
129     fCurrentEvent = event;
130   }
131   
132   if (fRecoTrackStore != 0x0) return fRecoTrackStore;
133   else {
134     if (!fESDTree) return 0x0;
135     if (fESDTree->GetEvent(event) <= 0) {
136       AliError(Form("fails to read ESD object for event %d", event));
137       return 0x0;
138     }
139     MakeReconstructedTracks();
140     return fRecoTrackStore;
141   }
142 }
143
144 //_____________________________________________________________________________
145 AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
146 {
147   /// Return a track store containing the track references (converted into 
148   /// MUONTrack objects) for a given event
149   if (event != fCurrentEvent) {
150     ResetStores();
151     fCurrentEvent = event;
152   }
153   
154   if (fTrackRefStore != 0x0) return fTrackRefStore;
155   else {
156     if (!fMCEventHandler->GetEvent(event)) {
157       AliError(Form("fails to read MC objects for event %d", event));
158       return 0x0;
159     }
160     MakeTrackRefs();
161     return fTrackRefStore;
162   }
163 }
164
165 //_____________________________________________________________________________
166 AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
167 {
168   /// Return a track store containing the reconstructible tracks for a given event
169   if (event != fCurrentEvent) {
170     ResetStores();
171     fCurrentEvent = event;
172   }
173   
174   if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
175   else {
176     if (TrackRefs(event) == 0x0) return 0x0;
177     MakeReconstructibleTracks();
178     return fRecoTrackRefStore;
179   }
180 }
181
182 //_____________________________________________________________________________
183 void AliMUONRecoCheck::MakeReconstructedTracks()
184 {
185   /// Make reconstructed tracks
186   if (!(fRecoTrackStore = AliMUONESDInterface::NewTrackStore())) return;
187   
188   // loop over all reconstructed tracks and add them to the store (skip ghosts)
189   Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
190   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
191     AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
192     if (esdTrack->ContainTrackerData()) AliMUONESDInterface::Add(*esdTrack, *fRecoTrackStore);
193   }
194   
195 }
196
197 //_____________________________________________________________________________
198 void AliMUONRecoCheck::MakeTrackRefs()
199 {
200   /// Make reconstructible tracks
201   AliMUONVTrackStore *tmpTrackRefStore = AliMUONESDInterface::NewTrackStore();
202   if (!tmpTrackRefStore) return;
203   
204   Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
205   TParticle* particle;
206   TClonesArray* trackRefs;
207   Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
208   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
209   if (!cStore) return;
210   AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
211   
212   // loop over simulated tracks
213   for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) {
214     Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
215     
216     // skip empty trackRefs
217     if (nHits < 1) continue;
218     
219     // get the particle charge for further calculation
220     TParticlePDG* ppdg = particle->GetPDG();
221     Int_t charge = ppdg != NULL ? (Int_t)(ppdg->Charge()/3.0) : 0;
222     
223     AliMUONTrack track;
224     
225     // loop over simulated track hits
226     for (Int_t iHit = 0; iHit < nHits; ++iHit) {        
227       AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
228       
229       // skip trackRefs not in MUON
230       if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
231     
232       // Get track parameters of current hit
233       x = trackReference->X();
234       y = trackReference->Y();
235       z = trackReference->Z();
236       pX = trackReference->Px();
237       pY = trackReference->Py();
238       pZ = trackReference->Pz();
239       
240       // check chamberId of current trackReference
241       Int_t chamberId = AliMUONConstants::ChamberNumber(z);
242       if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
243       
244       // set hit parameters
245       hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(chamberId, 0, 0));
246       hit->SetXYZ(x,y,z);
247       hit->SetErrXY(0.,0.);
248       
249       // compute track parameters at hit
250       bendingSlope = 0;
251       nonBendingSlope = 0;
252       inverseBendingMomentum = 0;
253       if (TMath::Abs(pZ) > 0) {
254         bendingSlope = pY/pZ;
255         nonBendingSlope = pX/pZ;
256       }
257       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
258       if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
259       inverseBendingMomentum *= charge;
260       
261       // set track parameters at hit
262       AliMUONTrackParam trackParam;
263       trackParam.SetNonBendingCoor(x);
264       trackParam.SetBendingCoor(y);
265       trackParam.SetZ(z);
266       trackParam.SetBendingSlope(bendingSlope);
267       trackParam.SetNonBendingSlope(nonBendingSlope);
268       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
269       
270       // add track parameters at current hit to the track
271       track.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
272     }
273     
274     // if none of the track hits was in MUON, goto the next track
275     if (track.GetNClusters() < 1) continue;
276     
277     // get track parameters at particle's vertex
278     x = particle->Vx();
279     y = particle->Vy();
280     z = particle->Vz();
281     pX = particle->Px();
282     pY = particle->Py();
283     pZ = particle->Pz();
284     
285     // compute rest of track parameters at particle's vertex
286     bendingSlope = 0;
287     nonBendingSlope = 0;
288     inverseBendingMomentum = 0;
289     if (TMath::Abs(pZ) > 0) {
290       bendingSlope = pY/pZ;
291       nonBendingSlope = pX/pZ;
292     }
293     Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
294     if (pYZ >0) inverseBendingMomentum = 1/pYZ;
295     inverseBendingMomentum *= charge;
296     
297     // set track parameters at particle's vertex
298     AliMUONTrackParam trackParamAtVertex;
299     trackParamAtVertex.SetNonBendingCoor(x);
300     trackParamAtVertex.SetBendingCoor(y);
301     trackParamAtVertex.SetZ(z);
302     trackParamAtVertex.SetBendingSlope(bendingSlope);
303     trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
304     trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
305     
306     // add track parameters at vertex
307     track.SetTrackParamAtVertex(&trackParamAtVertex);
308     
309     // store the track
310     track.SetTrackID(iTrackRef);
311     tmpTrackRefStore->Add(track);
312   }
313   
314   CleanMuonTrackRef(tmpTrackRefStore);
315   
316   delete hit;
317   delete cStore;
318   delete tmpTrackRefStore;
319 }
320
321 //_____________________________________________________________________________
322 void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
323 {
324   /// Re-calculate hits parameters because two AliTrackReferences are recorded for
325   /// each chamber (one when particle is entering + one when particle is leaving 
326   /// the sensitive volume) 
327   if (!(fTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
328   
329   Double_t maxGasGap = 1.; // cm 
330   Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
331   Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
332   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
333   if (!cStore) return;
334   AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
335   
336   // create iterator
337   TIter next(tmpTrackRefStore->CreateIterator());
338   
339   // loop over tmpTrackRef
340   AliMUONTrack* track;
341   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
342     
343     AliMUONTrack newTrack;
344     
345     // loop over tmpTrackRef's hits
346     Int_t iHit1 = 0;
347     Int_t nTrackHits = track->GetNClusters();
348     while (iHit1 < nTrackHits) {
349       AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
350       
351       // get track parameters at hit1
352       x1  = trackParam1->GetNonBendingCoor();
353       y1  = trackParam1->GetBendingCoor();
354       z1  = trackParam1->GetZ();
355       pX1 = trackParam1->Px();
356       pY1 = trackParam1->Py();
357       pZ1 = trackParam1->Pz();
358       
359       // prepare new track parameters
360       x  = x1;
361       y  = y1;
362       z  = z1;
363       pX = pX1;
364       pY = pY1;
365       pZ = pZ1;
366       
367       // loop over next tmpTrackRef's hits
368       Int_t nCombinedHits = 1;
369       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
370         AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
371         
372         // get z position of hit2
373         z2 = trackParam2->GetZ();
374         
375         // complete new track parameters if hit2 is on the same detection element
376         if ( TMath::Abs(z2-z1) < maxGasGap ) {
377           x  += trackParam2->GetNonBendingCoor();
378           y  += trackParam2->GetBendingCoor();
379           z  += z2;
380           pX += trackParam2->Px();
381           pY += trackParam2->Py();
382           pZ += trackParam2->Pz();
383           nCombinedHits++;
384           iHit1 = iHit2;
385         }
386         
387       }
388       
389       // finalize new track parameters
390       x  /= (Double_t)nCombinedHits;
391       y  /= (Double_t)nCombinedHits;
392       z  /= (Double_t)nCombinedHits;
393       pX /= (Double_t)nCombinedHits;
394       pY /= (Double_t)nCombinedHits;
395       pZ /= (Double_t)nCombinedHits;
396       bendingSlope = 0;
397       nonBendingSlope = 0;
398       inverseBendingMomentum = 0;
399       if (TMath::Abs(pZ) > 0) {
400         bendingSlope = pY/pZ;
401         nonBendingSlope = pX/pZ;
402       }
403       Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
404       if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
405       inverseBendingMomentum *= trackParam1->GetCharge();
406       
407       // set hit parameters
408       hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(trackParam1->GetClusterPtr()->GetChamberId(), 0, 0));
409       hit->SetXYZ(x,y,z);
410       hit->SetErrXY(0.,0.);
411       
412       // set new track parameters at new hit
413       AliMUONTrackParam trackParam;
414       trackParam.SetNonBendingCoor(x);
415       trackParam.SetBendingCoor(y);
416       trackParam.SetZ(z);
417       trackParam.SetBendingSlope(bendingSlope);
418       trackParam.SetNonBendingSlope(nonBendingSlope);
419       trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
420       
421       // add track parameters at current hit to the track
422       newTrack.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
423       
424       iHit1++;
425     }
426     
427     newTrack.SetTrackID(track->GetTrackID());
428     newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
429     fTrackRefStore->Add(newTrack);
430     
431   }
432   
433   delete hit;
434   delete cStore;
435 }
436
437 //_____________________________________________________________________________
438 void AliMUONRecoCheck::MakeReconstructibleTracks()
439 {
440   /// Isolate the reconstructible tracks
441   if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
442   
443   // create iterator on trackRef
444   TIter next(fTrackRefStore->CreateIterator());
445   
446   // loop over trackRef
447   AliMUONTrack* track;
448   while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
449     
450     Bool_t* chamberInTrack = new Bool_t(AliMUONConstants::NTrackingCh());
451     for (Int_t iCh = 0; iCh < AliMUONConstants::NTrackingCh(); iCh++) chamberInTrack[iCh] = kFALSE;
452     
453     // loop over trackRef's hits to get hit chambers
454     Int_t nTrackHits = track->GetNClusters();
455     for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
456       AliMUONVCluster* hit = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit))->GetClusterPtr(); 
457       chamberInTrack[hit->GetChamberId()] = kTRUE;
458     } 
459     
460     // track is reconstructible if the particle is depositing a hit
461     // in the following chamber combinations:
462     Bool_t trackOK = kTRUE;
463     if (!chamberInTrack[0] && !chamberInTrack[1]) trackOK = kFALSE;
464     if (!chamberInTrack[2] && !chamberInTrack[3]) trackOK = kFALSE;
465     if (!chamberInTrack[4] && !chamberInTrack[5]) trackOK = kFALSE;
466     Int_t nHitsInLastStations = 0;
467     for (Int_t iCh = 6; iCh < AliMUONConstants::NTrackingCh(); iCh++)
468       if (chamberInTrack[iCh]) nHitsInLastStations++; 
469     if(nHitsInLastStations < 3) trackOK = kFALSE;
470     
471     // Add reconstructible tracks to fRecoTrackRefStore
472     if (trackOK) fRecoTrackRefStore->Add(*track);
473     
474     delete [] chamberInTrack;
475   }
476
477 }
478