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