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