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