Coding conventions
[u/mrichter/AliRoot.git] / MUON / AliMUONRecoCheck.cxx
CommitLineData
b8dc484b 1/**************************************************************************
48217459 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**************************************************************************/
b8dc484b 15
d19b6003 16/* $Id$ */
17
3d1463c8 18//-----------------------------------------------------------------------------
5398f946 19/// \class AliMUONRecoCheck
78649106 20/// Utility class to check reconstruction
5398f946 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.
3d1463c8 24//-----------------------------------------------------------------------------
b8dc484b 25
cc056de1 26#include "AliMUONRecoCheck.h"
96ebe67e 27#include "AliMUONRawClusterV2.h"
cc056de1 28#include "AliMUONTrack.h"
cc056de1 29#include "AliMUONConstants.h"
48217459 30#include "AliMUONDataInterface.h"
96ebe67e 31#include "AliMUONTrackStoreV1.h"
32
33#include "AliMCEventHandler.h"
34#include "AliMCEvent.h"
b8dc484b 35#include "AliStack.h"
cc056de1 36#include "AliTrackReference.h"
37#include "AliLog.h"
96ebe67e 38
cc056de1 39#include <TParticle.h>
5b1afea0 40#include <TParticlePDG.h>
b8dc484b 41
96ebe67e 42#include <Riostream.h>
43
5398f946 44/// \cond CLASSIMP
b8dc484b 45ClassImp(AliMUONRecoCheck)
5398f946 46/// \endcond
b8dc484b 47
48//_____________________________________________________________________________
96ebe67e 49AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader, Char_t *pathSim)
48217459 50: TObject(),
96ebe67e 51fMCEventHandler(new AliMCEventHandler()),
52fDataInterface(new AliMUONDataInterface(chLoader)),
53fCurrentEvent(0),
54fTrackRefStore(0x0),
55fRecoTrackRefStore(0x0),
56fRecoTrackStore(0x0)
068e4c36 57{
48217459 58 /// Normal ctor
96ebe67e 59 fMCEventHandler->SetInputPath(pathSim);
60 fMCEventHandler->InitIO("");
068e4c36 61}
62
63//_____________________________________________________________________________
48217459 64AliMUONRecoCheck::~AliMUONRecoCheck()
b8dc484b 65{
48217459 66 /// Destructor
96ebe67e 67 delete fMCEventHandler;
48217459 68 delete fDataInterface;
96ebe67e 69 ResetStores();
70}
71
72//_____________________________________________________________________________
73void AliMUONRecoCheck::ResetStores()
74{
75 /// Deletes all the store objects that have been created and resets the pointers to 0x0
76 delete fTrackRefStore; fTrackRefStore = 0x0;
77 delete fRecoTrackRefStore; fRecoTrackRefStore = 0x0;
78 delete fRecoTrackStore; fRecoTrackStore = 0x0;
79}
80
81//_____________________________________________________________________________
82Int_t AliMUONRecoCheck::NumberOfEvents() const
83{
84 /// Return the number of events
85 if (fDataInterface->IsValid()) return fDataInterface->NumberOfEvents();
86 return 0;
b8dc484b 87}
88
48217459 89//_____________________________________________________________________________
96ebe67e 90AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event)
48217459 91{
92 /// Return the reconstructed track store for a given event
96ebe67e 93 if (!fDataInterface->IsValid()) return new AliMUONTrackStoreV1();
94
95 if (event != fCurrentEvent) {
96 ResetStores();
97 fCurrentEvent = event;
98 }
99
100 if (fRecoTrackStore != 0x0) return fRecoTrackStore;
101
102 fRecoTrackStore = new AliMUONTrackStoreV1();
103
104 return fRecoTrackStore;
48217459 105}
45a05f07 106
b8dc484b 107//_____________________________________________________________________________
96ebe67e 108AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
b8dc484b 109{
48217459 110 /// Return a track store containing the track references (converted into
111 /// MUONTrack objects) for a given event
96ebe67e 112 if (event != fCurrentEvent) {
113 ResetStores();
114 fCurrentEvent = event;
115 }
116
117 if (fTrackRefStore != 0x0) return fTrackRefStore;
118 else {
119 if (!fMCEventHandler->GetEvent(event)) return 0x0;
120 MakeTrackRefs();
121 return fTrackRefStore;
122 }
b8dc484b 123}
124
125//_____________________________________________________________________________
96ebe67e 126AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
b8dc484b 127{
48217459 128 /// Return a track store containing the reconstructible tracks for a given event
96ebe67e 129 if (event != fCurrentEvent) {
130 ResetStores();
131 fCurrentEvent = event;
132 }
133
134 if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
135 else {
136 if (TrackRefs(event) == 0x0) return 0x0;
137 MakeReconstructibleTracks();
138 return fRecoTrackRefStore;
139 }
48217459 140}
b8dc484b 141
48217459 142//_____________________________________________________________________________
96ebe67e 143void AliMUONRecoCheck::MakeTrackRefs()
48217459 144{
145 /// Make reconstructible tracks
96ebe67e 146 AliMUONVTrackStore *tmpTrackRefStore = new AliMUONTrackStoreV1();
48217459 147
96ebe67e 148 Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
149 TParticle* particle;
150 TClonesArray* trackRefs;
151 Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
b8dc484b 152
96ebe67e 153 // loop over simulated tracks
154 for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; ++iTrackRef) {
155 Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
48217459 156
96ebe67e 157 // skip empty trackRefs
158 if (nHits < 1) continue;
159
160 // skip trackRefs not in MUON
161 AliTrackReference* trackReference0 = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(0));
162 if (trackReference0->DetectorId() != AliTrackReference::kMUON) continue;
163
164 AliMUONTrack track;
b8dc484b 165
96ebe67e 166 // get track parameters at particle's vertex
167 x = particle->Vx();
168 y = particle->Vy();
169 z = particle->Vz();
170 pX = particle->Px();
171 pY = particle->Py();
172 pZ = particle->Pz();
dba3c809 173
96ebe67e 174 // compute rest of track parameters at particle's vertex
175 bendingSlope = 0;
176 nonBendingSlope = 0;
177 inverseBendingMomentum = 0;
178 if (TMath::Abs(pZ) > 0) {
179 bendingSlope = pY/pZ;
180 nonBendingSlope = pX/pZ;
181 }
182 Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
183 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
184 TParticlePDG* ppdg = particle->GetPDG(1);
185 Int_t charge = (Int_t)(ppdg->Charge()/3.0);
186 inverseBendingMomentum *= charge;
187
188 // set track parameters at particle's vertex
189 AliMUONTrackParam trackParamAtVertex;
190 trackParamAtVertex.SetNonBendingCoor(x);
191 trackParamAtVertex.SetBendingCoor(y);
192 trackParamAtVertex.SetZ(z);
193 trackParamAtVertex.SetBendingSlope(bendingSlope);
194 trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
195 trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
196
197 // add track parameters at vertex
198 track.SetTrackParamAtVertex(&trackParamAtVertex);
48217459 199
96ebe67e 200 // loop over simulated track hits
201 for (Int_t iHit = 0; iHit < nHits; ++iHit) {
202 AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
068e4c36 203
96ebe67e 204 // Get track parameters of current hit
205 x = trackReference->X();
206 y = trackReference->Y();
207 z = trackReference->Z();
208 pX = trackReference->Px();
209 pY = trackReference->Py();
210 pZ = trackReference->Pz();
b8dc484b 211
96ebe67e 212 // check chamberId of current trackReference
213 Int_t chamberId = AliMUONConstants::ChamberNumber(z);
214 if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
48217459 215
96ebe67e 216 // set hit parameters
217 AliMUONRawClusterV2 hit(chamberId, 0, 0);
218 hit.SetXYZ(x,y,z);
219 hit.SetErrXY(0.,0.);
b8dc484b 220
96ebe67e 221 // compute track parameters at hit
222 Double_t bendingSlope = 0;
223 Double_t nonBendingSlope = 0;
224 Double_t inverseBendingMomentum = 0;
225 if (TMath::Abs(pZ) > 0) {
226 bendingSlope = pY/pZ;
227 nonBendingSlope = pX/pZ;
b8dc484b 228 }
96ebe67e 229 Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
230 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
231 inverseBendingMomentum *= charge;
232
233 // set track parameters at hit
234 AliMUONTrackParam trackParam;
235 trackParam.SetNonBendingCoor(x);
236 trackParam.SetBendingCoor(y);
237 trackParam.SetZ(z);
238 trackParam.SetBendingSlope(bendingSlope);
239 trackParam.SetNonBendingSlope(nonBendingSlope);
240 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
241
242 // add track parameters at current hit to the track
243 track.AddTrackParamAtCluster(trackParam,hit,kTRUE);
244 }
245
246 track.GetTrackParamAtCluster()->Sort();
247 track.SetTrackID(iTrackRef);
248 tmpTrackRefStore->Add(track);
b8dc484b 249 }
b8dc484b 250
96ebe67e 251 CleanMuonTrackRef(tmpTrackRefStore);
48217459 252
96ebe67e 253 delete tmpTrackRefStore;
b8dc484b 254}
255
256//_____________________________________________________________________________
96ebe67e 257void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
b8dc484b 258{
48217459 259 /// Re-calculate hits parameters because two AliTrackReferences are recorded for
260 /// each chamber (one when particle is entering + one when particle is leaving
261 /// the sensitive volume)
96ebe67e 262 fTrackRefStore = new AliMUONTrackStoreV1();
48217459 263
96ebe67e 264 Double_t maxGasGap = 1.; // cm
265 Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
266 Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
48217459 267
96ebe67e 268 // create iterator
269 TIter next(tmpTrackRefStore->CreateIterator());
48217459 270
96ebe67e 271 // loop over tmpTrackRef
48217459 272 AliMUONTrack* track;
96ebe67e 273 while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
274
275 AliMUONTrack newTrack;
276
277 // loop over tmpTrackRef's hits
278 Int_t iHit1 = 0;
279 Int_t nTrackHits = track->GetNClusters();
280 while (iHit1 < nTrackHits) {
281 AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
282
283 // get track parameters at hit1
284 x1 = trackParam1->GetNonBendingCoor();
285 y1 = trackParam1->GetBendingCoor();
286 z1 = trackParam1->GetZ();
287 pX1 = trackParam1->Px();
288 pY1 = trackParam1->Py();
289 pZ1 = trackParam1->Pz();
290
291 // prepare new track parameters
292 x = x1;
293 y = y1;
294 z = z1;
295 pX = pX1;
296 pY = pY1;
297 pZ = pZ1;
298
299 // loop over next tmpTrackRef's hits
300 Int_t nCombinedHits = 1;
301 for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
302 AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
48217459 303
96ebe67e 304 // get z position of hit2
305 z2 = trackParam2->GetZ();
306
307 // complete new track parameters if hit2 is on the same detection element
308 if ( TMath::Abs(z2-z1) < maxGasGap ) {
309 x += trackParam2->GetNonBendingCoor();
310 y += trackParam2->GetBendingCoor();
311 z += z2;
312 pX += trackParam2->Px();
313 pY += trackParam2->Py();
314 pZ += trackParam2->Pz();
315 nCombinedHits++;
48217459 316 iHit1 = iHit2;
317 }
318
96ebe67e 319 }
320
321 // finalize new track parameters
322 x /= (Double_t)nCombinedHits;
323 y /= (Double_t)nCombinedHits;
324 z /= (Double_t)nCombinedHits;
325 pX /= (Double_t)nCombinedHits;
326 pY /= (Double_t)nCombinedHits;
327 pZ /= (Double_t)nCombinedHits;
328 bendingSlope = 0;
329 nonBendingSlope = 0;
330 inverseBendingMomentum = 0;
331 if (TMath::Abs(pZ) > 0) {
332 bendingSlope = pY/pZ;
333 nonBendingSlope = pX/pZ;
334 }
335 Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
336 if (pYZ >0) inverseBendingMomentum = 1/pYZ;
337 inverseBendingMomentum *= trackParam1->GetCharge();
b8dc484b 338
96ebe67e 339 // set hit parameters
340 AliMUONRawClusterV2 hit(trackParam1->GetClusterPtr()->GetChamberId(), 0, 0);
341 hit.SetXYZ(x,y,z);
342 hit.SetErrXY(0.,0.);
48217459 343
96ebe67e 344 // set new track parameters at new hit
345 AliMUONTrackParam trackParam;
346 trackParam.SetNonBendingCoor(x);
347 trackParam.SetBendingCoor(y);
348 trackParam.SetZ(z);
349 trackParam.SetBendingSlope(bendingSlope);
350 trackParam.SetNonBendingSlope(nonBendingSlope);
351 trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
352
353 // add track parameters at current hit to the track
354 newTrack.AddTrackParamAtCluster(trackParam,hit,kTRUE);
b8dc484b 355
356 iHit1++;
96ebe67e 357 }
48217459 358
96ebe67e 359 newTrack.GetTrackParamAtCluster()->Sort();
360 newTrack.SetTrackID(track->GetTrackID());
361 newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
362 fTrackRefStore->Add(newTrack);
b8dc484b 363
96ebe67e 364 }
b8dc484b 365
b8dc484b 366}
367
368//_____________________________________________________________________________
96ebe67e 369void AliMUONRecoCheck::MakeReconstructibleTracks()
b8dc484b 370{
48217459 371 /// Calculate the number of reconstructible tracks
96ebe67e 372 fRecoTrackRefStore = new AliMUONTrackStoreV1();
48217459 373
96ebe67e 374 // create iterator on trackRef
375 TIter next(fTrackRefStore->CreateIterator());
b8dc484b 376
96ebe67e 377 // loop over trackRef
48217459 378 AliMUONTrack* track;
96ebe67e 379 while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
380
381 Bool_t* chamberInTrack = new Bool_t(AliMUONConstants::NTrackingCh());
382 for (Int_t iCh = 0; iCh < AliMUONConstants::NTrackingCh(); iCh++) chamberInTrack[iCh] = kFALSE;
48217459 383
96ebe67e 384 // loop over trackRef's hits to get hit chambers
385 Int_t nTrackHits = track->GetNClusters();
386 for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
387 AliMUONVCluster* hit = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit))->GetClusterPtr();
388 chamberInTrack[hit->GetChamberId()] = kTRUE;
cc056de1 389 }
96ebe67e 390
cc056de1 391 // track is reconstructible if the particle is depositing a hit
392 // in the following chamber combinations:
96ebe67e 393 Bool_t trackOK = kTRUE;
394 if (!chamberInTrack[0] && !chamberInTrack[1]) trackOK = kFALSE;
395 if (!chamberInTrack[2] && !chamberInTrack[3]) trackOK = kFALSE;
396 if (!chamberInTrack[4] && !chamberInTrack[5]) trackOK = kFALSE;
397 Int_t nHitsInLastStations = 0;
398 for (Int_t iCh = 6; iCh < AliMUONConstants::NTrackingCh(); iCh++)
399 if (chamberInTrack[iCh]) nHitsInLastStations++;
400 if(nHitsInLastStations < 3) trackOK = kFALSE;
48217459 401
96ebe67e 402 // Add reconstructible tracks to fRecoTrackRefStore
403 if (trackOK) fRecoTrackRefStore->Add(*track);
48217459 404
96ebe67e 405 delete [] chamberInTrack;
b8dc484b 406 }
48217459 407
b8dc484b 408}
409