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