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