change class name: AliMUONEventReconstructor -> AliMUONTrackReconstructor
[u/mrichter/AliRoot.git] / MUON / AliMUONReconstructor.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 /* $Id$ */
16
17 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
19 // class for MUON reconstruction                                             //
20 //                                                                           //
21 ///////////////////////////////////////////////////////////////////////////////
22 #include <TParticle.h>
23 #include <TArrayF.h>
24
25 #include "AliRunLoader.h"
26 #include "AliHeader.h"
27 #include "AliGenEventHeader.h"
28 #include "AliESD.h"
29 #include "AliMUONReconstructor.h"
30  
31 #include "AliMUONData.h"
32 #include "AliMUONTrackReconstructor.h"
33 #include "AliMUONClusterReconstructor.h"
34 #include "AliMUONClusterFinderVS.h"
35 #include "AliMUONTrack.h"
36 #include "AliMUONTrackParam.h"
37 #include "AliMUONTriggerTrack.h"
38 #include "AliESDMuonTrack.h"
39 #include "AliMUONRawData.h"
40
41 #include "AliRawReader.h"
42
43
44 ClassImp(AliMUONReconstructor)
45 //_____________________________________________________________________________
46 AliMUONReconstructor::AliMUONReconstructor()
47   : AliReconstructor()
48 {
49 }
50 //_____________________________________________________________________________
51 AliMUONReconstructor::~AliMUONReconstructor()
52 {
53 }
54 //_____________________________________________________________________________
55 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
56 {
57 //  AliLoader
58   AliLoader* loader = runLoader->GetLoader("MUONLoader");
59   Int_t nEvents = runLoader->GetNumberOfEvents();
60
61 // used local container for each method
62 // passing fLoader as argument, could be avoided ???
63   AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader);
64   AliMUONData* dataEvent = recoEvent->GetMUONData();
65
66   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader);
67   AliMUONData* dataCluster = recoCluster->GetMUONData();
68   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
69   recModel->SetGhostChi2Cut(10);
70
71   loader->LoadDigits("READ");
72   loader->LoadRecPoints("RECREATE");
73   loader->LoadTracks("RECREATE");
74   
75   //   Loop over events              
76   for(Int_t ievent = 0; ievent < nEvents; ievent++) {
77     printf("Event %d\n",ievent);
78     runLoader->GetEvent(ievent);
79
80     //----------------------- digit2cluster & Trigger2Trigger -------------------
81     if (!loader->TreeR()) loader->MakeRecPointsContainer();
82      
83     // tracking branch
84     dataCluster->MakeBranch("RC");
85     dataCluster->SetTreeAddress("D,RC");
86     recoCluster->Digits2Clusters(); 
87     dataCluster->Fill("RC"); 
88
89     // trigger branch
90     dataCluster->MakeBranch("TC");
91     dataCluster->SetTreeAddress("TC");
92     recoCluster->Trigger2Trigger(); 
93     dataCluster->Fill("TC");
94
95     loader->WriteRecPoints("OVERWRITE");
96
97     //---------------------------- Track & TriggerTrack ---------------------
98     if (!loader->TreeT()) loader->MakeTracksContainer();
99
100     // trigger branch
101     dataEvent->MakeBranch("RL"); //trigger track
102     dataEvent->SetTreeAddress("RL");
103     recoEvent->EventReconstructTrigger();
104     dataEvent->Fill("RL");
105
106     // tracking branch
107     dataEvent->MakeBranch("RT"); //track
108     dataEvent->SetTreeAddress("RT");
109     recoEvent->EventReconstruct();
110     dataEvent->Fill("RT");
111
112     loader->WriteTracks("OVERWRITE"); 
113   
114     //--------------------------- Resetting branches -----------------------
115     dataCluster->ResetDigits();
116     dataCluster->ResetRawClusters();
117     dataCluster->ResetTrigger();
118
119     dataEvent->ResetRawClusters();
120     dataEvent->ResetTrigger();
121     dataEvent->ResetRecTracks();  
122     dataEvent->ResetRecTriggerTracks();
123
124   }
125   loader->UnloadDigits();
126   loader->UnloadRecPoints();
127   loader->UnloadTracks();
128
129   delete recoCluster;
130   delete recoEvent;
131 }
132
133 //_____________________________________________________________________________
134 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawReader) const
135 {
136 //  AliLoader
137   AliLoader* loader = runLoader->GetLoader("MUONLoader");
138
139 // used local container for each method
140 // passing fLoader as argument, could be avoided ???
141   AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader);
142   AliMUONData* dataEvent = recoEvent->GetMUONData();
143
144   AliMUONRawData* rawData = new AliMUONRawData(loader);
145   AliMUONData* dataCluster = rawData->GetMUONData();
146
147   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader, dataCluster);
148   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
149   recModel->SetGhostChi2Cut(10);
150
151   loader->LoadRecPoints("RECREATE");
152   loader->LoadTracks("RECREATE");
153   
154   //   Loop over events  
155   Int_t iEvent = 0;
156             
157   while (rawReader->NextEvent()) {
158     printf("Event %d\n",iEvent);
159     runLoader->GetEvent(iEvent++);
160     rawData->Raw2Digits(rawReader);
161
162     //----------------------- digit2cluster & Trigger2Trigger -------------------
163     if (!loader->TreeR()) loader->MakeRecPointsContainer();
164      
165     // tracking branch
166     dataCluster->MakeBranch("RC");
167     dataCluster->SetTreeAddress("D,RC");
168     recoCluster->Digits2Clusters(); 
169     dataCluster->Fill("RC"); 
170
171     // trigger branch
172     dataCluster->MakeBranch("TC");
173     dataCluster->SetTreeAddress("TC");
174     recoCluster->Trigger2Trigger(); 
175     dataCluster->Fill("TC");
176
177     loader->WriteRecPoints("OVERWRITE");
178
179     //---------------------------- Track & TriggerTrack ---------------------
180     if (!loader->TreeT()) loader->MakeTracksContainer();
181
182     // trigger branch
183     dataEvent->MakeBranch("RL"); //trigger track
184     dataEvent->SetTreeAddress("RL");
185     recoEvent->EventReconstructTrigger();
186     dataEvent->Fill("RL");
187
188     // tracking branch
189     dataEvent->MakeBranch("RT"); //track
190     dataEvent->SetTreeAddress("RT");
191     recoEvent->EventReconstruct();
192     dataEvent->Fill("RT");
193
194     loader->WriteTracks("OVERWRITE");  
195   
196     //--------------------------- Resetting branches -----------------------
197     dataCluster->ResetDigits();
198     dataCluster->ResetRawClusters();
199     dataCluster->ResetTrigger();
200
201     dataEvent->ResetRawClusters();
202     dataEvent->ResetTrigger();
203     dataEvent->ResetRecTracks();
204     dataEvent->ResetRecTriggerTracks();
205   
206   }
207   loader->UnloadRecPoints();
208   loader->UnloadTracks();
209
210   delete recoCluster;
211   delete recoEvent;
212 }
213
214 //_____________________________________________________________________________
215 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
216 {
217   TClonesArray* recTracksArray = 0;
218   TClonesArray* recTrigTracksArray = 0;
219   
220   AliLoader* loader = runLoader->GetLoader("MUONLoader");
221   loader->LoadTracks("READ");
222   AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");
223
224    // declaration  
225   Int_t iEvent, nPart;
226   Int_t nTrackHits, nPrimary;
227   Double_t fitFmin;
228   TArrayF vertex(3);
229
230   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
231   Double_t xRec, yRec, zRec, chi2MatchTrigger;
232   Bool_t matchTrigger;
233
234   // setting pointer for tracks, triggertracks & trackparam at vertex
235   AliMUONTrack* recTrack = 0;
236   AliMUONTrackParam* trackParam = 0;
237   AliMUONTriggerTrack* recTriggerTrack = 0;
238   TParticle* particle = new TParticle();
239   AliGenEventHeader* header = 0;
240   iEvent = runLoader->GetEventNumber(); 
241   runLoader->GetEvent(iEvent);
242
243   // vertex calculation (maybe it exists already somewhere else)
244   vertex[0] = vertex[1] = vertex[2] = 0.;
245   nPrimary = 0;
246   if ( (header = runLoader->GetHeader()->GenEventHeader()) ) {
247     header->PrimaryVertex(vertex);
248   } else {
249     runLoader->LoadKinematics("READ");
250     runLoader->TreeK()->GetBranch("Particles")->SetAddress(&particle);
251     nPart = (Int_t)runLoader->TreeK()->GetEntries();
252     for(Int_t iPart = 0; iPart < nPart; iPart++) {
253       runLoader->TreeK()->GetEvent(iPart);
254       if (particle->GetFirstMother() == -1) {
255         vertex[0] += particle->Vx();
256         vertex[1] += particle->Vy();
257         vertex[2] += particle->Vz();
258         nPrimary++;
259       }
260       if (nPrimary) {
261         vertex[0] /= (double)nPrimary;
262         vertex[1] /= (double)nPrimary;
263         vertex[2] /= (double)nPrimary;
264       }
265     }
266   }
267   // setting ESD MUON class
268   AliESDMuonTrack* theESDTrack = new  AliESDMuonTrack() ;
269
270   //-------------------- trigger tracks-------------
271   Long_t trigPat = 0;
272   muonData->SetTreeAddress("RL");
273   muonData->GetRecTriggerTracks();
274   recTrigTracksArray = muonData->RecTriggerTracks();
275
276   // ready global trigger pattern from first track
277   if (recTrigTracksArray) 
278     recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
279   if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
280
281   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
282  
283   // -------------------- tracks-------------
284   muonData->SetTreeAddress("RT");
285   muonData->GetRecTracks();
286   recTracksArray = muonData->RecTracks();
287         
288   Int_t nRecTracks = 0;
289   if (recTracksArray)
290     nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
291   
292   // loop over tracks
293   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
294
295     // reading info from tracks
296     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
297
298     trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
299     trackParam->ExtrapToVertex(vertex[0],vertex[1],vertex[2]);
300
301     bendingSlope            = trackParam->GetBendingSlope();
302     nonBendingSlope         = trackParam->GetNonBendingSlope();
303     inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
304     xRec  = trackParam->GetNonBendingCoor();
305     yRec  = trackParam->GetBendingCoor();
306     zRec  = trackParam->GetZ();
307
308     nTrackHits       = recTrack->GetNTrackHits();
309     fitFmin          = recTrack->GetFitFMin();
310     matchTrigger     = recTrack->GetMatchTrigger();
311     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
312
313     // setting data member of ESD MUON
314     theESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
315     theESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
316     theESDTrack->SetThetaY(TMath::ATan(bendingSlope));
317     theESDTrack->SetZ(zRec);
318     theESDTrack->SetBendingCoor(yRec); // calculate vertex at ESD or Tracking level ?
319     theESDTrack->SetNonBendingCoor(xRec);
320     theESDTrack->SetChi2(fitFmin);
321     theESDTrack->SetNHit(nTrackHits);
322     theESDTrack->SetMatchTrigger(matchTrigger);
323     theESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
324
325     // storing ESD MUON Track into ESD Event 
326     if (nRecTracks != 0)  
327       esd->AddMuonTrack(theESDTrack);
328   } // end loop tracks
329
330   // add global trigger pattern
331   if (nRecTracks != 0)  
332     esd->SetTrigger(trigPat);
333
334   // reset muondata
335   muonData->ResetRecTracks();
336   muonData->ResetRecTriggerTracks();
337
338   //} // end loop on event  
339   loader->UnloadTracks(); 
340   if (!header)
341     runLoader->UnloadKinematics();
342   delete theESDTrack;
343   delete muonData;
344   // delete particle;
345 }