Coding conventions (Laurent)
[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
23 #include "AliMUONReconstructor.h"
24
25 #include "AliESD.h"
26 #include "AliESDMuonTrack.h"
27 #include "AliLog.h"
28 #include "AliMUON.h"
29 #include "AliMUONCalibrationData.h"
30 #include "AliMUONClusterFinderAZ.h"
31 #include "AliMUONClusterFinderVS.h"
32 #include "AliMUONClusterReconstructor.h"
33 #include "AliMUONData.h"
34 #include "AliMUONDigitCalibrator.h"
35 #include "AliMUONEventRecoCombi.h" 
36 #include "AliMUONRawReader.h"
37 #include "AliMUONTrack.h"
38 #include "AliMUONTrackParam.h"
39 #include "AliMUONTrackReconstructor.h"
40 #include "AliMUONTriggerTrack.h"
41 #include "AliRawReader.h"
42 #include "AliRun.h"
43 #include "AliRunLoader.h"
44 #include "TTask.h"
45 #include "TStopwatch.h"
46
47 ClassImp(AliMUONReconstructor)
48
49 //_____________________________________________________________________________
50 AliMUONReconstructor::AliMUONReconstructor()
51   : AliReconstructor(), fCalibrationData(0x0)
52 {
53 /// Default constructor
54
55     AliDebug(1,"");
56 }
57
58 //______________________________________________________________________________
59 AliMUONReconstructor::AliMUONReconstructor(const AliMUONReconstructor& right) 
60   : AliReconstructor(right) 
61 {  
62 /// Protected copy constructor (not implemented)
63
64   AliFatal("Copy constructor not provided.");
65 }
66
67 //_____________________________________________________________________________
68 AliMUONReconstructor::~AliMUONReconstructor()
69 {
70 /// Destructor
71
72   AliDebug(1,"");
73   delete fCalibrationData;
74 }
75
76 //______________________________________________________________________________
77 AliMUONReconstructor& 
78 AliMUONReconstructor::operator=(const AliMUONReconstructor& right)
79 {
80 /// Protected assignement operator (not implemented)
81
82   // check assignement to self
83   if (this == &right) return *this;
84
85   AliFatal("Assignement operator not provided.");
86     
87   return *this;  
88 }    
89
90 //_____________________________________________________________________________
91 TTask* 
92 AliMUONReconstructor::GetCalibrationTask(AliMUONData* data) const
93 {
94 /// Create the calibration task(s). 
95   
96   const AliRun* run = fRunLoader->GetAliRun();
97   
98   // Not really clean, but for the moment we must check whether the
99   // simulation has decalibrated the data or not...
100   const AliMUON* muon = static_cast<AliMUON*>(run->GetModule("MUON"));
101   if ( muon->DigitizerType().Contains("NewDigitizer") )
102   {
103     AliInfo("Calibration will occur.");
104     Int_t runNumber = run->GetRunNumber();     
105     fCalibrationData = new AliMUONCalibrationData(runNumber);
106     if ( !fCalibrationData->IsValid() )
107     {
108       AliError("Could not retrieve calibrations !");
109       delete fCalibrationData;
110       fCalibrationData = 0x0;
111       return 0x0;
112     }    
113     TTask* calibration = new TTask("MUONCalibrator","MUON Digit calibrator");
114     calibration->Add(new AliMUONDigitCalibrator(data,fCalibrationData));
115     //FIXME: calibration->Add(something about dead channels should go here).
116     return calibration;
117   }
118   else
119   {
120     AliInfo("Detected the usage of old digitizer (w/o decalibration). "
121             "Will not calibrate then.");
122     return 0x0;
123   }
124 }
125
126 //_____________________________________________________________________________
127 void
128 AliMUONReconstructor::Init(AliRunLoader* runLoader)
129 {
130 /// Initialize
131
132   fRunLoader = runLoader;
133 }
134
135 //_____________________________________________________________________________
136 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
137 {
138 /// Reconstruct
139 /// \todo add more
140
141   AliLoader* loader = runLoader->GetLoader("MUONLoader");
142   Int_t nEvents = runLoader->GetNumberOfEvents();
143
144   AliMUONData* data = new AliMUONData(loader,"MUON","MUON");
145
146 // passing loader as argument.
147   AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader, data);
148
149   if (strstr(GetOption(),"Original")) 
150     recoEvent->SetTrackMethod(1); // Original tracking
151   else if (strstr(GetOption(),"Combi")) 
152     recoEvent->SetTrackMethod(3); // Combined cluster / track
153   else
154     recoEvent->SetTrackMethod(2); // Kalman
155
156   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(data);
157   
158   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
159
160   if (!strstr(GetOption(),"VS")) {
161     recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
162     recoCluster->SetRecoModel(recModel);
163   }
164   recModel->SetGhostChi2Cut(10);
165
166   loader->LoadDigits("READ");
167   loader->LoadRecPoints("RECREATE");
168   loader->LoadTracks("RECREATE");
169   
170   TTask* calibration = GetCalibrationTask(data);
171   
172   Int_t chBeg = recoEvent->GetTrackMethod() == 3 ? 6 : 0; 
173   //   Loop over events              
174   for(Int_t ievent = 0; ievent < nEvents; ievent++) {
175
176     AliDebug(1,Form("Event %d",ievent));
177     
178     runLoader->GetEvent(ievent);
179
180     //----------------------- digit2cluster & Trigger2Trigger -------------------
181     if (!loader->TreeR()) loader->MakeRecPointsContainer();
182      
183     // tracking branch
184     if (recoEvent->GetTrackMethod() != 3) {
185       data->MakeBranch("RC");
186       data->SetTreeAddress("D,RC");
187     } else {
188       data->SetTreeAddress("D");
189       data->SetTreeAddress("RCC");
190     }
191     // Important for avoiding a memory leak when reading digits ( to be investigated more in detail)
192     // In any case the reading of GLT is needed for the Trigger2Tigger method below
193     data->SetTreeAddress("GLT");
194
195     data->GetDigits();
196     
197     if ( calibration ) 
198     {
199       calibration->ExecuteTask();
200     }
201     
202     recoCluster->Digits2Clusters(chBeg); 
203     
204     if (recoEvent->GetTrackMethod() == 3) {
205       // Combined cluster / track finder
206       AliMUONEventRecoCombi::Instance()->FillEvent(data, (AliMUONClusterFinderAZ*)recModel);
207       ((AliMUONClusterFinderAZ*) recModel)->SetReco(2); 
208     }
209     else data->Fill("RC"); 
210
211     // trigger branch
212     data->MakeBranch("TC");
213     data->SetTreeAddress("TC");
214     recoCluster->Trigger2Trigger(); 
215     data->Fill("TC");
216
217     //AZ loader->WriteRecPoints("OVERWRITE");
218
219     //---------------------------- Track & TriggerTrack ---------------------
220     if (!loader->TreeT()) loader->MakeTracksContainer();
221
222     // trigger branch
223     data->MakeBranch("RL"); //trigger track
224     data->SetTreeAddress("RL");
225     recoEvent->EventReconstructTrigger();
226     data->Fill("RL");
227
228     // tracking branch
229     data->MakeBranch("RT"); //track
230     data->SetTreeAddress("RT");
231     recoEvent->EventReconstruct();
232     data->Fill("RT");
233
234     loader->WriteTracks("OVERWRITE"); 
235   
236     if (recoEvent->GetTrackMethod() == 3) { 
237       // Combined cluster / track
238       ((AliMUONClusterFinderAZ*) recModel)->SetReco(1);
239       data->MakeBranch("RC");
240       data->SetTreeAddress("RC");
241       AliMUONEventRecoCombi::Instance()->FillRecP(data, recoEvent); 
242       data->Fill("RC"); 
243     }
244     loader->WriteRecPoints("OVERWRITE"); 
245
246     //--------------------------- Resetting branches -----------------------
247     data->ResetDigits();
248     data->ResetRawClusters();
249     data->ResetTrigger();
250
251     data->ResetRawClusters();
252     data->ResetTrigger();
253     data->ResetRecTracks();  
254     data->ResetRecTriggerTracks();
255
256   }
257   loader->UnloadDigits();
258   loader->UnloadRecPoints();
259   loader->UnloadTracks();
260
261   delete recoCluster;
262   delete recoEvent;
263   delete data;
264   delete calibration;
265 }
266
267 //_____________________________________________________________________________
268 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawReader) const
269 {
270 /// Recontruct
271 /// \todo add more
272
273   //  AliLoader
274   AliLoader* loader = runLoader->GetLoader("MUONLoader");
275   AliMUONData data(loader,"MUON","MUON");
276
277   // passing loader as argument.
278   AliMUONTrackReconstructor recoEvent(loader, &data);
279
280   AliMUONRawReader rawData(&data);
281
282   AliMUONClusterReconstructor recoCluster(&data);
283
284   if (strstr(GetOption(),"Original")) 
285   {
286     recoEvent.SetTrackMethod(1); // Original tracking
287   }
288   else
289   {
290     recoEvent.SetTrackMethod(2); // Kalman
291   }
292   
293   AliMUONClusterFinderVS *recModel = recoCluster.GetRecoModel();
294   if (!strstr(GetOption(),"VS")) 
295   {
296     recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
297     recoCluster.SetRecoModel(recModel);
298   }
299   recModel->SetGhostChi2Cut(10);
300
301   TTask* calibration = GetCalibrationTask(&data);
302   
303   loader->LoadRecPoints("RECREATE");
304   loader->LoadTracks("RECREATE");
305   loader->LoadDigits("READ");
306   
307   //   Loop over events  
308   Int_t iEvent = 0;
309            
310   TStopwatch totalTimer;
311   TStopwatch rawTimer;
312   TStopwatch calibTimer;
313   TStopwatch clusterTimer;
314   TStopwatch trackingTimer;
315   
316   rawTimer.Start(kTRUE); rawTimer.Stop();
317   calibTimer.Start(kTRUE); calibTimer.Stop();
318   clusterTimer.Start(kTRUE); clusterTimer.Stop();
319   trackingTimer.Start(kTRUE); trackingTimer.Stop();
320   
321   totalTimer.Start(kTRUE);
322   
323   while (rawReader->NextEvent()) 
324   {
325     AliDebug(1,Form("Event %d",iEvent));
326     
327     runLoader->GetEvent(iEvent++);
328
329     //----------------------- raw2digits & raw2trigger-------------------
330     if (!loader->TreeD()) 
331     {
332       AliDebug(1,Form("Making Digit Container for event %d",iEvent));
333       loader->MakeDigitsContainer();
334     }
335     
336     data.SetTreeAddress("D,GLT");
337     rawTimer.Start(kFALSE);
338     rawData.Raw2Digits(rawReader);
339     rawTimer.Stop();
340     
341     if ( calibration )
342     {
343       calibTimer.Start(kFALSE);
344       calibration->ExecuteTask();
345       calibTimer.Stop();
346     }
347   
348     //----------------------- digit2cluster & Trigger2Trigger -------------------
349     clusterTimer.Start(kFALSE);
350
351     if (!loader->TreeR()) loader->MakeRecPointsContainer();
352      
353     // tracking branch
354     data.MakeBranch("RC");
355     data.SetTreeAddress("RC");
356     recoCluster.Digits2Clusters(); 
357     data.Fill("RC"); 
358
359     // trigger branch
360     data.MakeBranch("TC");
361     data.SetTreeAddress("TC");
362 //    recoCluster.Trigger2Trigger(); 
363     data.Fill("TC");
364     
365     loader->WriteRecPoints("OVERWRITE");
366
367     clusterTimer.Stop();
368
369     //---------------------------- Track & TriggerTrack ---------------------
370     trackingTimer.Start(kFALSE);
371     if (!loader->TreeT()) loader->MakeTracksContainer();
372
373     // trigger branch
374     data.MakeBranch("RL"); //trigger track
375     data.SetTreeAddress("RL");
376     recoEvent.EventReconstructTrigger();
377     data.Fill("RL");
378
379     // tracking branch
380     data.MakeBranch("RT"); //track
381     data.SetTreeAddress("RT");
382     recoEvent.EventReconstruct();
383     data.Fill("RT");
384
385     loader->WriteTracks("OVERWRITE");  
386     trackingTimer.Stop();
387     
388     //--------------------------- Resetting branches -----------------------
389     data.ResetDigits();
390     data.ResetRawClusters();
391     data.ResetTrigger();
392
393     data.ResetRawClusters();
394     data.ResetTrigger();
395     data.ResetRecTracks();
396     data.ResetRecTriggerTracks();
397   
398   }
399   
400   totalTimer.Stop();
401   
402   loader->UnloadRecPoints();
403   loader->UnloadTracks();
404   loader->UnloadDigits();
405
406   AliInfo(Form("Execution time for converting RAW data to digits in MUON : R:%.2fs C:%.2fs",
407                rawTimer.RealTime(),rawTimer.CpuTime()));
408   AliInfo(Form("Execution time for calibrating MUON : R:%.2fs C:%.2fs",
409                calibTimer.RealTime(),calibTimer.CpuTime()));
410   AliInfo(Form("Execution time for clusterizing MUON : R:%.2fs C:%.2fs",
411                clusterTimer.RealTime(),clusterTimer.CpuTime()));
412   AliInfo(Form("Execution time for tracking MUON : R:%.2fs C:%.2fs",
413                trackingTimer.RealTime(),trackingTimer.CpuTime()));
414   AliInfo(Form("Total Execution time for Reconstruct(from raw) MUON : R:%.2fs C:%.2fs",
415                totalTimer.RealTime(),totalTimer.CpuTime()));
416 }
417
418 //_____________________________________________________________________________
419 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
420 {
421 /// Fill ESD
422 /// \todo add more
423
424   TClonesArray* recTracksArray = 0;
425   TClonesArray* recTrigTracksArray = 0;
426   
427   AliLoader* loader = runLoader->GetLoader("MUONLoader");
428   loader->LoadTracks("READ");
429   AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");
430
431    // declaration  
432   Int_t iEvent;// nPart;
433   Int_t nTrackHits;// nPrimary;
434   Double_t fitFmin;
435
436   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
437   Double_t xRec, yRec, zRec, chi2MatchTrigger;
438   Bool_t matchTrigger;
439
440   // setting pointer for tracks, triggertracks & trackparam at vertex
441   AliMUONTrack* recTrack = 0;
442   AliMUONTrackParam* trackParam = 0;
443   AliMUONTriggerTrack* recTriggerTrack = 0;
444 //   TParticle* particle = new TParticle();
445 //   AliGenEventHeader* header = 0;
446   iEvent = runLoader->GetEventNumber(); 
447   runLoader->GetEvent(iEvent);
448
449   // Get vertex 
450   Double_t vertex[3] = {0};
451   const AliESDVertex *esdVert = esd->GetVertex(); 
452   if (esdVert) esdVert->GetXYZ(vertex);
453   
454   //  nPrimary = 0;
455 //   if ( (header = runLoader->GetHeader()->GenEventHeader()) ) {
456 //     header->PrimaryVertex(vertex);
457 //   } else {
458 //     runLoader->LoadKinematics("READ");
459 //     runLoader->TreeK()->GetBranch("Particles")->SetAddress(&particle);
460 //     nPart = (Int_t)runLoader->TreeK()->GetEntries();
461 //     for(Int_t iPart = 0; iPart < nPart; iPart++) {
462 //       runLoader->TreeK()->GetEvent(iPart);
463 //       if (particle->GetFirstMother() == -1) {
464 //      vertex[0] += particle->Vx();
465 //      vertex[1] += particle->Vy();
466 //      vertex[2] += particle->Vz();
467 //      nPrimary++;
468 //       }
469 //       if (nPrimary) {
470 //      vertex[0] /= (double)nPrimary;
471 //      vertex[1] /= (double)nPrimary;
472 //      vertex[2] /= (double)nPrimary;
473 //       }
474 //     }
475 //   }
476   // setting ESD MUON class
477   AliESDMuonTrack* theESDTrack = new  AliESDMuonTrack() ;
478
479   //-------------------- trigger tracks-------------
480   Long_t trigPat = 0;
481   muonData->SetTreeAddress("RL");
482   muonData->GetRecTriggerTracks();
483   recTrigTracksArray = muonData->RecTriggerTracks();
484
485   // ready global trigger pattern from first track
486   if (recTrigTracksArray) 
487     recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
488   if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
489
490   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
491  
492   // -------------------- tracks-------------
493   muonData->SetTreeAddress("RT");
494   muonData->GetRecTracks();
495   recTracksArray = muonData->RecTracks();
496         
497   Int_t nRecTracks = 0;
498   if (recTracksArray)
499     nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
500   
501   // loop over tracks
502   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
503
504     // reading info from tracks
505     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
506
507     trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
508     trackParam->ExtrapToVertex(vertex[0],vertex[1],vertex[2]);
509
510     bendingSlope            = trackParam->GetBendingSlope();
511     nonBendingSlope         = trackParam->GetNonBendingSlope();
512     inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
513     xRec  = trackParam->GetNonBendingCoor();
514     yRec  = trackParam->GetBendingCoor();
515     zRec  = trackParam->GetZ();
516
517     nTrackHits       = recTrack->GetNTrackHits();
518     fitFmin          = recTrack->GetFitFMin();
519     matchTrigger     = recTrack->GetMatchTrigger();
520     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
521
522     // setting data member of ESD MUON
523     theESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
524     theESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
525     theESDTrack->SetThetaY(TMath::ATan(bendingSlope));
526     theESDTrack->SetZ(zRec);
527     theESDTrack->SetBendingCoor(yRec); // calculate vertex at ESD or Tracking level ?
528     theESDTrack->SetNonBendingCoor(xRec);
529     theESDTrack->SetChi2(fitFmin);
530     theESDTrack->SetNHit(nTrackHits);
531     theESDTrack->SetMatchTrigger(matchTrigger);
532     theESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
533
534     // storing ESD MUON Track into ESD Event 
535     if (nRecTracks != 0)  
536       esd->AddMuonTrack(theESDTrack);
537   } // end loop tracks
538
539   // add global trigger pattern
540   if (nRecTracks != 0)  
541     esd->SetTrigger(trigPat);
542
543   // reset muondata
544   muonData->ResetRecTracks();
545   muonData->ResetRecTriggerTracks();
546
547   //} // end loop on event  
548   loader->UnloadTracks(); 
549  //  if (!header)
550 //     runLoader->UnloadKinematics();
551   delete theESDTrack;
552   delete muonData;
553   // delete particle;
554 }//_____________________________________________________________________________
555 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliRawReader* /*rawReader*/, AliESD* esd) const
556 {
557 /// Fill ESD
558 /// \todo add more
559
560   // don't need rawReader ???
561   FillESD(runLoader, esd);
562 }