Correcting MUONGenerateGeometryData.C macro description
[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   loader->LoadDigits("RECREATE");
154
155
156   //   Loop over events  
157   Int_t iEvent = 0;
158             
159   while (rawReader->NextEvent()) {
160     printf("Event %d\n",iEvent);
161     runLoader->GetEvent(iEvent++);
162
163     //----------------------- raw2digits & raw2trigger-------------------
164     if (!loader->TreeD()) loader->MakeDigitsContainer();
165
166     // tracking branch
167     dataCluster->MakeBranch("D");
168     dataCluster->SetTreeAddress("D");
169     rawData->ReadTrackerDDL(rawReader);
170     dataCluster->Fill("D"); 
171
172     // trigger branch
173     dataCluster->MakeBranch("GLT");
174     dataCluster->SetTreeAddress("GLT");
175     rawData->ReadTriggerDDL(rawReader);
176     dataCluster->Fill("GLT"); 
177
178     loader->WriteDigits("OVERWRITE");
179
180     //----------------------- digit2cluster & Trigger2Trigger -------------------
181     if (!loader->TreeR()) loader->MakeRecPointsContainer();
182      
183     // tracking branch
184     dataCluster->MakeBranch("RC");
185     dataCluster->SetTreeAddress("RC");
186     recoCluster->Digits2Clusters(); 
187     dataCluster->Fill("RC"); 
188
189     // trigger branch
190     dataCluster->MakeBranch("TC");
191     dataCluster->SetTreeAddress("TC");
192     recoCluster->Trigger2Trigger(); 
193     dataCluster->Fill("TC");
194
195     loader->WriteRecPoints("OVERWRITE");
196
197     //---------------------------- Track & TriggerTrack ---------------------
198     if (!loader->TreeT()) loader->MakeTracksContainer();
199
200     // trigger branch
201     dataEvent->MakeBranch("RL"); //trigger track
202     dataEvent->SetTreeAddress("RL");
203     recoEvent->EventReconstructTrigger();
204     dataEvent->Fill("RL");
205
206     // tracking branch
207     dataEvent->MakeBranch("RT"); //track
208     dataEvent->SetTreeAddress("RT");
209     recoEvent->EventReconstruct();
210     dataEvent->Fill("RT");
211
212     loader->WriteTracks("OVERWRITE");  
213   
214     //--------------------------- Resetting branches -----------------------
215     dataCluster->ResetDigits();
216     dataCluster->ResetRawClusters();
217     dataCluster->ResetTrigger();
218
219     dataEvent->ResetRawClusters();
220     dataEvent->ResetTrigger();
221     dataEvent->ResetRecTracks();
222     dataEvent->ResetRecTriggerTracks();
223   
224   }
225   loader->UnloadRecPoints();
226   loader->UnloadTracks();
227   loader->UnloadDigits();
228
229   delete recoCluster;
230   delete recoEvent;
231 }
232
233 //_____________________________________________________________________________
234 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
235 {
236   TClonesArray* recTracksArray = 0;
237   TClonesArray* recTrigTracksArray = 0;
238   
239   AliLoader* loader = runLoader->GetLoader("MUONLoader");
240   loader->LoadTracks("READ");
241   AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");
242
243    // declaration  
244   Int_t iEvent;// nPart;
245   Int_t nTrackHits;// nPrimary;
246   Double_t fitFmin;
247   TArrayF vertex(3);
248
249   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
250   Double_t xRec, yRec, zRec, chi2MatchTrigger;
251   Bool_t matchTrigger;
252
253   // setting pointer for tracks, triggertracks & trackparam at vertex
254   AliMUONTrack* recTrack = 0;
255   AliMUONTrackParam* trackParam = 0;
256   AliMUONTriggerTrack* recTriggerTrack = 0;
257 //   TParticle* particle = new TParticle();
258 //   AliGenEventHeader* header = 0;
259   iEvent = runLoader->GetEventNumber(); 
260   runLoader->GetEvent(iEvent);
261
262   // vertex calculation (maybe it exists already somewhere else)
263   vertex[0] = vertex[1] = vertex[2] = 0.;
264  //  nPrimary = 0;
265 //   if ( (header = runLoader->GetHeader()->GenEventHeader()) ) {
266 //     header->PrimaryVertex(vertex);
267 //   } else {
268 //     runLoader->LoadKinematics("READ");
269 //     runLoader->TreeK()->GetBranch("Particles")->SetAddress(&particle);
270 //     nPart = (Int_t)runLoader->TreeK()->GetEntries();
271 //     for(Int_t iPart = 0; iPart < nPart; iPart++) {
272 //       runLoader->TreeK()->GetEvent(iPart);
273 //       if (particle->GetFirstMother() == -1) {
274 //      vertex[0] += particle->Vx();
275 //      vertex[1] += particle->Vy();
276 //      vertex[2] += particle->Vz();
277 //      nPrimary++;
278 //       }
279 //       if (nPrimary) {
280 //      vertex[0] /= (double)nPrimary;
281 //      vertex[1] /= (double)nPrimary;
282 //      vertex[2] /= (double)nPrimary;
283 //       }
284 //     }
285 //   }
286   // setting ESD MUON class
287   AliESDMuonTrack* theESDTrack = new  AliESDMuonTrack() ;
288
289   //-------------------- trigger tracks-------------
290   Long_t trigPat = 0;
291   muonData->SetTreeAddress("RL");
292   muonData->GetRecTriggerTracks();
293   recTrigTracksArray = muonData->RecTriggerTracks();
294
295   // ready global trigger pattern from first track
296   if (recTrigTracksArray) 
297     recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
298   if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
299
300   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
301  
302   // -------------------- tracks-------------
303   muonData->SetTreeAddress("RT");
304   muonData->GetRecTracks();
305   recTracksArray = muonData->RecTracks();
306         
307   Int_t nRecTracks = 0;
308   if (recTracksArray)
309     nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
310   
311   // loop over tracks
312   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
313
314     // reading info from tracks
315     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
316
317     trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
318     trackParam->ExtrapToVertex(vertex[0],vertex[1],vertex[2]);
319
320     bendingSlope            = trackParam->GetBendingSlope();
321     nonBendingSlope         = trackParam->GetNonBendingSlope();
322     inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
323     xRec  = trackParam->GetNonBendingCoor();
324     yRec  = trackParam->GetBendingCoor();
325     zRec  = trackParam->GetZ();
326
327     nTrackHits       = recTrack->GetNTrackHits();
328     fitFmin          = recTrack->GetFitFMin();
329     matchTrigger     = recTrack->GetMatchTrigger();
330     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
331
332     // setting data member of ESD MUON
333     theESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
334     theESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
335     theESDTrack->SetThetaY(TMath::ATan(bendingSlope));
336     theESDTrack->SetZ(zRec);
337     theESDTrack->SetBendingCoor(yRec); // calculate vertex at ESD or Tracking level ?
338     theESDTrack->SetNonBendingCoor(xRec);
339     theESDTrack->SetChi2(fitFmin);
340     theESDTrack->SetNHit(nTrackHits);
341     theESDTrack->SetMatchTrigger(matchTrigger);
342     theESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
343
344     // storing ESD MUON Track into ESD Event 
345     if (nRecTracks != 0)  
346       esd->AddMuonTrack(theESDTrack);
347   } // end loop tracks
348
349   // add global trigger pattern
350   if (nRecTracks != 0)  
351     esd->SetTrigger(trigPat);
352
353   // reset muondata
354   muonData->ResetRecTracks();
355   muonData->ResetRecTriggerTracks();
356
357   //} // end loop on event  
358   loader->UnloadTracks(); 
359  //  if (!header)
360 //     runLoader->UnloadKinematics();
361   delete theESDTrack;
362   delete muonData;
363   // delete particle;
364 }//_____________________________________________________________________________
365 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliRawReader* /*rawReader*/, AliESD* esd) const
366 {
367   // don't need rawReader ???
368   FillESD(runLoader, esd);
369 }