Removing warnings ...
[u/mrichter/AliRoot.git] / MUON / AliMUONRecoEvent.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 //Authors: Mihaela Gheata, Andrei Gheata 09/10/00
19 ////////////////////////////////////////////////////////////////////
20 //                                                                //
21 // AliMUONRecoEvent, AliMUONRecoTrack (and AliMUONRecoDisplay)    //
22 //                                                                //
23 // Theses classes are used to store and retrieve                  //
24 // MUON reconstructed events.                                     //
25 // The corresponding tree is constructed and filled               //
26 // during the FillEvent() method of AliMUONEventReconstructor,    //
27 // when all reconstruction information is available.              //
28 //                                                                //
29 ////////////////////////////////////////////////////////////////////
30 ////////////////////////////////////////////////////////////////////
31 //                                                                //
32 // AliMUONRecoEvent                                               //
33 //                                                                //
34 // This class handles an array of reconstructed tracks.           //
35 // It provides :                                                  //
36 //      - filling the tracks array according to the information   //
37 //        stored in AliMUONEventReconstructor class ;             //
38 //      - printing event and track informations : event numer,    //
39 //        number of tracks, hits positions, reconstr. mometum.    //
40 //                                                                //
41 ////////////////////////////////////////////////////////////////////
42 ////////////////////////////////////////////////////////////////////
43 //                                                                //
44 // AliMUONRecoTrack                                               //
45 //                                                                //
46 // This class represents a reconstructed muon track               //
47 // in the tree of reconstructed events.                           //
48 //                                                                //
49 ////////////////////////////////////////////////////////////////////
50
51 #include <Riostream.h>
52 #include <AliRun.h>
53 #include <TClonesArray.h>
54 #include <TClass.h>
55
56 #include <TFile.h>
57 #include <TMatrixD.h>
58 #include <TParticle.h>
59
60 #include "AliMUONRecoEvent.h"
61 #include "AliMUONEventReconstructor.h"
62 #include "AliMUONTrack.h"
63 #include "AliMUONTrackK.h"
64 #include "AliMUONTrackParam.h"
65 #include "AliMUONHitForRec.h"
66 #include "AliMUONTrackHit.h"
67 #include "AliHeader.h"
68
69 ClassImp(AliMUONRecoTrack)
70 ClassImp(AliMUONRecoEvent)
71
72 //-------------------------------------------------------------------
73 AliMUONRecoEvent::AliMUONRecoEvent(Int_t eventNo) 
74 {
75 // Reconstructed event constructor
76    fTracks      = new TClonesArray("AliMUONRecoTrack",200);
77    fNevr        = eventNo;
78    fNtracks = 0;
79 }
80
81 //-------------------------------------------------------------------
82 AliMUONRecoEvent::~AliMUONRecoEvent() 
83 {
84 // Destructor of AliMUONRecoEvent
85    fTracks->Delete();
86    delete fTracks;
87    fTracks = 0;
88 }
89
90 //-------------------------------------------------------------------
91 AliMUONRecoTrack* AliMUONRecoEvent::AddEmptyTrack()
92 {
93 // Add a empty AliMUONRecoTrackObject to the track list
94    TClonesArray &dumptracks = *fTracks;
95    return (new(dumptracks[fNtracks++])AliMUONRecoTrack(kTRUE));
96 }
97
98 //-------------------------------------------------------------------
99 void AliMUONRecoEvent::Clear(Option_t * /*option*/)
100 {
101 // Clears all track pointers from the list
102 //   fTracks->Clear(option);
103    fTracks->Delete();
104    fNtracks=0;
105 }
106
107 //-------------------------------------------------------------------
108 void AliMUONRecoEvent::EventInfo()
109 {
110 // Prints reconstructed event information
111    cout << "*********************Reco Dumper**************************" << endl;
112    cout << "Event number : " << fNevr << endl;
113    cout << "   Number of tracks : " << fNtracks << endl;
114    AliMUONRecoTrack *currentTrack =0;
115    Int_t trackIndex = 0;
116    for(trackIndex=0; trackIndex<fNtracks; trackIndex++) {
117       currentTrack = (AliMUONRecoTrack*)fTracks->UncheckedAt(trackIndex);
118       cout << "Track : " << trackIndex << endl;
119       cout << "   Sign : " << currentTrack->GetSign() << endl;
120       cout << "   Vertex position    : " << currentTrack->GetVertexPos() << endl;
121       Double_t momreco[3];
122       for (Int_t axis=0; axis<3; axis++) {
123          momreco[axis] = currentTrack->GetMomReconstr(axis);
124       }
125       cout << "   Reconstructed mom. : " << "Px=" << momreco[0] << "Py=" << momreco[1] << "Pz=" << momreco[2] << endl;
126       cout << "   Chi squared        : " << currentTrack->GetChi2r() << endl;
127       cout << "   Hits positions     : " << endl;
128       Double_t xhit, yhit, zhit;
129       for (Int_t chamber=0; chamber<10; chamber++) {
130          xhit = currentTrack->GetPosX(chamber);
131          yhit = currentTrack->GetPosY(chamber);
132          zhit = currentTrack->GetPosZ(chamber);
133 //         cout <<"      chamber" << chamber << " X=" << xhit << " Y=" << yhit << " Z=" << zhit << endl;
134       }  
135    }
136    cout << "**********************************************************" << endl;
137 }
138
139 //-------------------------------------------------------------------
140 Bool_t AliMUONRecoEvent::MakeDumpTracks(Int_t muons, TClonesArray *tracksPtr, 
141   AliMUONEventReconstructor *EventReco)
142 {
143 // This method takes the pointer of the list of reconstructed tracks from
144 // AliMUONEventReconstructor and fill the reconstructed AliMUONRecoEvent
145 // fields.
146
147         cout << "Enter MakeDumpTracks..." << endl;
148    Int_t nTracks = tracksPtr->GetEntriesFast();
149    cout << "nTracks = "<< nTracks << endl;
150    if (nTracks == 0) {
151       cout << "AliMUONRecoEvent::MakeDumpTracks: Number of tracks is zero !" << endl;
152       //AZ return kFALSE;
153    }
154    cout << tracksPtr << endl;
155    if (!tracksPtr) {
156       cout << "AliMUONRecoEvent::MakeDumpTracks() : You should call SetRecoTracksPtr() first..." << endl;
157       return kFALSE;
158    }
159         // Get event number
160    Int_t noEvent = gAlice->GetHeader()->GetEvent();
161    cout << "noEvent = "<< nTracks << endl;
162    tracksPtr->Compress();  // simple loop
163    AliMUONRecoTrack *currentTrack;
164    Int_t trackIndex, nTrackHits = 0;
165    Double_t z, pYZ, bendingSlope, nonBendingSlope;
166    Double_t pX, pY, pZ;                 // reconstructed momentum components
167    Int_t isign, flag=0;         // charge sign, flag of reconstructed track
168    Double_t alpha, beta;
169    TObjArray *hitsOnTrack = 0;
170    AliMUONTrackHit *trackHit = 0;
171    AliMUONTrack *track = 0;
172    AliMUONTrackK *trackK = 0;
173    TMatrixD *trackParamK; //AZnon
174    AliMUONTrackParam *trackParam = 0;
175    // Fill event number and number of tracks
176    fNevr = noEvent;
177    fMuons = muons; //AZ - number of muons within acceptance
178    // Loop over reconstructed tracks
179    for (trackIndex=0; trackIndex<nTracks; trackIndex++) {
180       cout << " trackIndex = " << trackIndex << endl;
181       currentTrack = AddEmptyTrack();
182       cout << " currentTrack = " << currentTrack << endl;
183
184       if (EventReco->GetTrackMethod() == 2) { // Kalman
185
186         trackK = (AliMUONTrackK*) ((*tracksPtr)[trackIndex]);
187         nTrackHits = trackK->GetNTrackHits();
188         trackParamK = trackK->GetTrackParameters();
189         isign = Int_t(TMath::Sign(1., (*trackParamK)(4,0)));
190         z = trackK->GetZ();
191         alpha = (*trackParamK)(2,0);
192         beta = (*trackParamK)(3,0);
193         pYZ = TMath::Cos(beta)/TMath::Abs((*trackParamK)(4,0));
194         pZ = pYZ/TMath::Cos(alpha);
195         pX = TMath::Sin(beta)/TMath::Abs((*trackParamK)(4,0));
196         pY = pYZ*TMath::Sin(alpha);
197
198         currentTrack->SetVertexPos(z);
199         currentTrack->SetMomReconstr(pX,pY,pZ);
200         currentTrack->SetSign(isign);
201         currentTrack->SetChi2r(trackK->GetTrackQuality());
202
203         // Check hits on the track
204         hitsOnTrack = trackK->GetHitOnTrack();
205         Float_t signal = 0;
206         Float_t tht = 0;
207         for (int ihit = 0; ihit < nTrackHits; ihit++) {
208           signal += ((AliMUONHitForRec*)((*hitsOnTrack)[ihit]))->GetGeantSignal();
209           tht += TMath::Min (1,((AliMUONHitForRec*)((*hitsOnTrack)[ihit]))->GetTHTrack());
210         }
211         signal /= nTrackHits;
212         tht /= nTrackHits;
213         flag = 0;
214         if (TMath::Nint(signal) > 0) { // signal muon
215           for (int ihit = 0; ihit < nTrackHits ; ihit++) {
216             if (((AliMUONHitForRec*)((*hitsOnTrack)[ihit]))->GetTHTrack() != TMath::Nint(tht)) flag++;
217           }
218         } else flag = -9; // background track
219         //cout << TMath::Nint(signal) << " " << TMath::Nint(tht) << " " << recTrackNt->fFlag << endl;
220         currentTrack->SetFlag(flag);
221       } else { // default tracking
222
223         track = (AliMUONTrack*) ((*tracksPtr)[trackIndex]);
224         nTrackHits = track->GetNTrackHits();
225         // track parameters at Vertex
226         trackParam = track->GetTrackParamAtVertex();
227         bendingSlope = trackParam->GetBendingSlope();
228         nonBendingSlope = trackParam->GetNonBendingSlope();
229
230         z = trackParam->GetZ();
231         pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
232         pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
233         pX = pZ * nonBendingSlope;
234         pY = pZ * bendingSlope;
235         
236         if (trackParam->GetInverseBendingMomentum()<0) isign=-1; else isign=1;
237         currentTrack->SetVertexPos(z);
238         currentTrack->SetMomReconstr(pX,pY,pZ);
239         currentTrack->SetSign(isign);
240         //         currentTrack->SetChi2r(trackParam->GetChi2());
241         currentTrack->SetChi2r(0);
242
243         // Check hits on the track
244         hitsOnTrack = track->GetTrackHitsPtr();
245         Float_t signal = 0;
246         Float_t tht = 0;
247         AliMUONHitForRec *hitForRec = 0;
248         for (int ihit = 0; ihit < nTrackHits; ihit++) {
249           hitForRec = ((AliMUONTrackHit*)(*hitsOnTrack)[ihit])->GetHitForRecPtr();
250           signal += hitForRec->GetGeantSignal();
251           tht += TMath::Min (1,hitForRec->GetTHTrack());
252         }
253         signal /= nTrackHits;
254         tht /= nTrackHits;
255         flag = 0;
256         if (TMath::Nint(signal) > 0) { // signal muon
257           for (int ihit = 0; ihit < nTrackHits ; ihit++) {
258             hitForRec = ((AliMUONTrackHit*)(*hitsOnTrack)[ihit])->GetHitForRecPtr();
259             if (hitForRec->GetTHTrack() != TMath::Nint(tht)) flag++;
260           }
261         } else flag = -9; // background track
262         //cout << TMath::Nint(signal) << " " << TMath::Nint(tht) << " " << recTrackNt->fFlag << endl;
263         currentTrack->SetFlag(flag);
264       }
265      
266       Double_t xhit,yhit,zhit;
267       // Loop over track hits
268       for (Int_t trackHitIndex = 0; trackHitIndex < nTrackHits; trackHitIndex++) {
269         if (EventReco->GetTrackMethod() == 2) { // Kalman
270           xhit = ((AliMUONHitForRec*)((*hitsOnTrack)[trackHitIndex]))->GetNonBendingCoor();
271           yhit = ((AliMUONHitForRec*)((*hitsOnTrack)[trackHitIndex]))->GetBendingCoor();
272           zhit = ((AliMUONHitForRec*)((*hitsOnTrack)[trackHitIndex]))->GetZ();    
273         } else {
274           trackHit = (AliMUONTrackHit*) (*(track->GetTrackHitsPtr()))[trackHitIndex];
275           xhit = trackHit->GetHitForRecPtr()->GetNonBendingCoor();
276           yhit = trackHit->GetHitForRecPtr()->GetBendingCoor();
277           zhit = trackHit->GetHitForRecPtr()->GetZ();
278         }
279         if (trackHitIndex >= 0 && trackHitIndex < 10) {
280           currentTrack->SetHitPosition(trackHitIndex,xhit,yhit,zhit);
281         } else { cout << "track " << trackIndex << " hit out of range" << endl;} 
282       }
283    
284    }
285    cout << "Leave MakeDumpTracks..." << endl;
286    return kTRUE;
287 }
288
289 //-------------------------------------------------------------------
290 void AliMUONRecoEvent::Streamer(TBuffer &R__b)
291 {
292 // Streams an object of class AliMUONRecoEvent
293    if (R__b.IsReading()) {
294       fTracks->Clear();
295       AliMUONRecoEvent::Class()->ReadBuffer(R__b, this);
296    } else {
297       cout << "...writing event to file...\n";
298       AliMUONRecoEvent::Class()->WriteBuffer(R__b, this);
299    }
300 }
301
302 //-------------------------------------------------------------------
303 AliMUONRecoTrack::AliMUONRecoTrack(Bool_t active)
304 {
305 //Constructor of AliMUONRecoTrack
306    fSign  = 0;
307    fZvr   = 0.0;
308    fChi2r = 0.0;
309    if (active) {
310         for (Int_t axis=0; axis<3; axis++) {
311         fPr[axis] = 0.0;
312         }
313         for (Int_t chamber=0; chamber<10; chamber++) {
314         fPosX[chamber] = 0.0;
315         fPosY[chamber] = 0.0;
316         fPosZ[chamber] = 0.0;
317         }
318    }
319 }
320
321 //-------------------------------------------------------------------
322 const Double_t AliMUONRecoTrack::Phi()
323 {
324 // Return trach phi angle
325         return TMath::ATan2(fPr[2], fPr[1]);
326 }
327
328 //-------------------------------------------------------------------
329 const Double_t AliMUONRecoTrack::Theta()
330 {
331 // Return trach theta angle
332    return TMath::ACos(fPr[2] / P());
333 }
334
335 //-------------------------------------------------------------------
336 void AliMUONRecoTrack::SetMomReconstr(Double_t px, Double_t py, Double_t pz)
337 {
338 // Set the track reconstructed momentum 
339    fPr[0] = px;
340    fPr[1] = py;
341    fPr[2] = pz;            
342
343    
344 //-------------------------------------------------------------------           
345 void AliMUONRecoTrack::SetHitPosition(Int_t chamber, Double_t x, Double_t y, Double_t z)
346 {
347 // Set hit coordinates in chamber[0..9]
348    fPosX[chamber] = x;
349    fPosY[chamber] = y;
350    fPosZ[chamber] = z;
351 }
352 //-------------------------------------------------------------------           
353 void AliMUONRecoTrack::TrackInfo()
354 {
355 // Prints momentum info for this track
356    cout << "Px=" << GetMomReconstr(0) << " Py=" << GetMomReconstr(1) <<
357           " Pz=" << GetMomReconstr(2) << " P=" << P() << endl;
358 }