Adding an option to disabling status map usage (for debug only)
[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 // Class AliMUONReconstructor  
19 //-----------------------------
20 // Class for the 
21 // MUON track reconstruction
22
23 #include "AliMUONReconstructor.h"
24
25 #include "AliMUONConstants.h"
26 #include "AliMUONCalibrationData.h"
27 #include "AliMUONClusterFinderAZ.h"
28 #include "AliMUONClusterReconstructor.h"
29 #include "AliMUONData.h"
30 #include "AliMUONDigitCalibrator.h"
31 #include "AliMUONEventRecoCombi.h" 
32 #include "AliMUONDigitMaker.h"
33 #include "AliMUONTrack.h"
34 #include "AliMUONTrackParam.h"
35 #include "AliMUONTrackExtrap.h"
36 #include "AliMUONVTrackReconstructor.h"
37 #include "AliMUONTrackReconstructor.h"
38 #include "AliMUONTrackReconstructorK.h"
39 #include "AliMUONTriggerTrack.h"
40 #include "AliMUONTriggerCircuit.h"
41 #include "AliMUONTriggerCrateStore.h"
42 #include "AliMUONSegFactory.h"
43 #include "AliMUONSegmentation.h"
44 #include "AliMUONPreClusterFinder.h"
45 #include "AliMUONClusterFinderCOG.h"
46 #include "AliMUONClusterFinderSimpleFit.h"
47 #include "AliMUONClusterFinderMLEM.h"
48
49 #include "AliESD.h"
50 #include "AliESDMuonTrack.h"
51 #include "AliLog.h"
52 #include "AliRawReader.h"
53 #include "AliRunLoader.h"
54 #include "AliCDBManager.h"
55
56 #include "TTask.h"
57 #include "TStopwatch.h"
58
59 /// \cond CLASSIMP
60 ClassImp(AliMUONReconstructor)
61 /// \endcond
62
63 //_____________________________________________________________________________
64 AliMUONReconstructor::AliMUONReconstructor()
65   : AliReconstructor(), 
66     fDigitMaker(new AliMUONDigitMaker()), 
67     fCalibrationData(0x0),
68     fCrateManager(new AliMUONTriggerCrateStore()),
69     fTriggerCircuit(new TClonesArray("AliMUONTriggerCircuit", 234)),
70     fTransformer(new AliMUONGeometryTransformer(kTRUE)),
71     fSegmentation(0x0)
72
73 {
74 /// Default constructor
75
76     AliDebug(1,"");
77     // Crate manager
78     fCrateManager->ReadFromFile();
79
80     // set to digit maker
81     fDigitMaker->SetCrateManager(fCrateManager);
82
83     // transformater
84     fTransformer->ReadGeometryData("volpath.dat", "geometry.root");
85     
86     // create segmentation and pass it to EventRecoCombi
87     AliMUONSegFactory factory(fTransformer);
88     fSegmentation = factory.CreateSegmentation();
89     AliMUONEventRecoCombi::Instance(fSegmentation); 
90
91     // trigger circuit
92     for (Int_t i = 0; i < AliMUONConstants::NTriggerCircuit(); i++)  {
93       AliMUONTriggerCircuit* c = new AliMUONTriggerCircuit();
94       c->SetTransformer(fTransformer);
95       c->Init(i,*fCrateManager);
96       TClonesArray& circuit = *fTriggerCircuit;
97       new(circuit[circuit.GetEntriesFast()])AliMUONTriggerCircuit(*c);
98       delete c;
99     }
100
101   
102 }
103
104 //_____________________________________________________________________________
105 AliMUONReconstructor::~AliMUONReconstructor()
106 {
107 /// Destructor
108
109   AliDebug(1,"");
110   delete fCalibrationData;
111   delete fDigitMaker;
112   delete fCrateManager;
113   delete fTriggerCircuit;
114   delete fTransformer;
115   delete fSegmentation;
116 }
117
118 //_____________________________________________________________________________
119 TTask* 
120 AliMUONReconstructor::GetCalibrationTask(AliMUONData* data) const
121 {
122 /// Create the calibration task(s). 
123   
124   //const AliRun* run = fRunLoader->GetAliRun();
125   //Int_t runNumber = run->GetRunNumber();     
126   Int_t runNumber = AliCDBManager::Instance()->GetRun();
127   AliInfo("Calibration will occur.");
128  
129   fCalibrationData = new AliMUONCalibrationData(runNumber);
130   if ( !fCalibrationData->IsValid() )
131     {
132       AliError("Could not retrieve calibrations !");
133       delete fCalibrationData;
134       fCalibrationData = 0x0;
135       return 0x0;
136     }    
137   TTask* calibration = new TTask("MUONCalibrator","MUON Digit calibrator");
138   
139   TString opt(GetOption());
140   opt.ToUpper();
141   Bool_t statusMap(kTRUE);
142   
143   if ( strstr(opt,"NOSTATUSMAP") )
144   {
145     AliWarning("Disconnecting status map : SHOULD BE USED FOR DEBUG ONLY. NOT FOR PRODUCTION !!!");
146     statusMap = kFALSE; 
147   }
148   calibration->Add(new AliMUONDigitCalibrator(data,fCalibrationData,statusMap));
149   return calibration;
150 }
151
152 //_____________________________________________________________________________
153 AliMUONClusterReconstructor*
154 AliMUONReconstructor::CreateClusterReconstructor(AliMUONData* data) const
155 {
156 /// Create cluster reconstructor
157
158   AliMUONVClusterFinder* clusterFinder(0x0);
159   
160   TString opt(GetOption());
161   opt.ToUpper();
162   
163   if ( strstr(opt,"PRECLUSTER") )
164   {
165     clusterFinder = new AliMUONPreClusterFinder;
166   }  
167   else if ( strstr(opt,"COG") )
168   {
169     clusterFinder = new AliMUONClusterFinderCOG;
170   }  
171   else if ( strstr(opt,"SIMPLEFIT") )
172   {
173     clusterFinder = new AliMUONClusterFinderSimpleFit;
174   }
175   else if ( strstr(opt,"MLEM:DRAW") )
176   {
177     clusterFinder = new AliMUONClusterFinderMLEM(kTRUE);
178   }
179   else if ( strstr(opt,"MLEM") )
180   {
181     clusterFinder = new AliMUONClusterFinderMLEM(kFALSE);
182   } 
183   
184   if ( clusterFinder) 
185   {
186     AliInfo(Form("Will use %s for clusterizing",clusterFinder->ClassName()));
187   }
188   
189   AliMUONClusterReconstructor* clusterReco = 
190     new AliMUONClusterReconstructor(data,clusterFinder,fTransformer);
191   return clusterReco;
192 }
193
194 //_____________________________________________________________________________
195 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
196 {
197 /// Reconstruct
198 /// \todo add more
199
200   AliLoader* loader = runLoader->GetLoader("MUONLoader");
201   Int_t nEvents     = runLoader->GetNumberOfEvents();
202   Int_t evtNumber   = runLoader->GetEventNumber();
203
204   AliMUONData* data = new AliMUONData(loader,"MUON","MUON");
205
206 // passing loader as argument.
207   AliMUONVTrackReconstructor* recoEvent;
208   if (strstr(GetOption(),"Original")) recoEvent = new AliMUONTrackReconstructor(data);
209   else if (strstr(GetOption(),"Combi")) recoEvent = new AliMUONTrackReconstructorK(data,"Combi");
210   else recoEvent = new AliMUONTrackReconstructorK(data,"Kalman");
211   
212   recoEvent->SetTriggerCircuit(fTriggerCircuit);
213
214   AliMUONClusterReconstructor* recoCluster = CreateClusterReconstructor(data);
215   
216   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
217
218   if (!strstr(GetOption(),"VS")) {
219     recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
220     recoCluster->SetRecoModel(recModel);
221   }
222   recModel->SetGhostChi2Cut(10);
223   recModel->SetEventNumber(evtNumber);
224
225   loader->LoadDigits("READ");
226   loader->LoadRecPoints("RECREATE");
227   loader->LoadTracks("RECREATE");
228   
229   TTask* calibration = GetCalibrationTask(data);
230   
231   Int_t chBeg = (strstr(GetOption(),"Combi") ? 6 : 0);
232
233   //   Loop over events              
234   for(Int_t ievent = 0; ievent < nEvents; ievent++) {
235
236     AliDebug(1,Form("Event %d",ievent));
237     
238     runLoader->GetEvent(ievent);
239
240     //----------------------- digit2cluster & Trigger2Trigger -------------------
241     if (!loader->TreeR()) loader->MakeRecPointsContainer();
242      
243     // tracking branch
244     if (!strstr(GetOption(),"Combi")) {
245       data->MakeBranch("RC");
246       data->SetTreeAddress("D,RC");
247     } else {
248       data->SetTreeAddress("D");
249       data->SetTreeAddress("RCC");
250     }
251     // Important for avoiding a memory leak when reading digits ( to be investigated more in detail)
252     // In any case the reading of GLT is needed for the Trigger2Tigger method below
253     data->SetTreeAddress("GLT");
254
255     data->GetDigits();
256     
257     if ( calibration ) 
258     {
259       calibration->ExecuteTask();
260     }
261     
262     recoCluster->Digits2Clusters(chBeg); 
263     
264     if (strstr(GetOption(),"Combi")) {
265       // Combined cluster / track finder
266       AliMUONEventRecoCombi::Instance()->FillEvent(data, (AliMUONClusterFinderAZ*)recModel);
267       ((AliMUONClusterFinderAZ*) recModel)->SetReco(2); 
268     }
269     else data->Fill("RC"); 
270
271     // trigger branch
272     data->MakeBranch("TC");
273     data->SetTreeAddress("TC");
274     recoCluster->Trigger2Trigger(); 
275     data->Fill("TC");
276
277     //AZ loader->WriteRecPoints("OVERWRITE");
278
279     //---------------------------- Track & TriggerTrack ---------------------
280     if (!loader->TreeT()) loader->MakeTracksContainer();
281
282     // trigger branch
283     data->MakeBranch("RL"); //trigger track
284     data->SetTreeAddress("RL");
285     recoEvent->EventReconstructTrigger();
286     data->Fill("RL");
287
288     // tracking branch
289     data->MakeBranch("RT"); //track
290     data->SetTreeAddress("RT");
291     recoEvent->EventReconstruct();
292     data->Fill("RT");
293
294     loader->WriteTracks("OVERWRITE"); 
295   
296     if (strstr(GetOption(),"Combi")) { 
297       // Combined cluster / track
298       ((AliMUONClusterFinderAZ*) recModel)->SetReco(1);
299       data->MakeBranch("RC");
300       data->SetTreeAddress("RC");
301       AliMUONEventRecoCombi::Instance()->FillRecP(data, (AliMUONTrackReconstructorK*)recoEvent); 
302       data->Fill("RC"); 
303     }
304     loader->WriteRecPoints("OVERWRITE"); 
305
306     //--------------------------- Resetting branches -----------------------
307     data->ResetDigits();
308     data->ResetRawClusters();
309     data->ResetTrigger();
310     data->ResetRecTracks();  
311     data->ResetRecTriggerTracks();
312
313   }
314   loader->UnloadDigits();
315   loader->UnloadRecPoints();
316   loader->UnloadTracks();
317
318   delete recoCluster;
319   delete recoEvent;
320   delete data;
321   delete calibration;
322 }
323
324 //_____________________________________________________________________________
325 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, 
326                                        AliRawReader* rawReader) const
327 {
328 /// Recontruct
329 /// \todo add more
330
331   //  AliLoader
332   AliLoader* loader = runLoader->GetLoader("MUONLoader");
333   Int_t evtNumber   = runLoader->GetEventNumber();
334  
335   AliMUONData data(loader,"MUON","MUON");
336
337   // passing loader as argument.
338   fDigitMaker->SetMUONData(&data);
339
340   // disable trigger rawdata reading
341   if (strstr(GetOption(),"TriggerDisable"))
342       fDigitMaker->DisableTrigger();
343
344   AliMUONClusterReconstructor* recoCluster = CreateClusterReconstructor(&data);
345
346   AliMUONVTrackReconstructor *recoEvent;
347   if (strstr(GetOption(),"Original")) recoEvent = new AliMUONTrackReconstructor(&data);
348   else if (strstr(GetOption(),"Combi")) recoEvent = new AliMUONTrackReconstructorK(&data,"Combi");
349   else recoEvent = new AliMUONTrackReconstructorK(&data,"Kalman");
350
351   recoEvent->SetTriggerCircuit(fTriggerCircuit);
352
353   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
354
355   if (!strstr(GetOption(),"VS")) 
356   {
357     recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
358     recoCluster->SetRecoModel(recModel);
359   }
360   recModel->SetGhostChi2Cut(10);
361   recModel->SetEventNumber(evtNumber);
362
363   TTask* calibration = GetCalibrationTask(&data);
364   
365   loader->LoadRecPoints("RECREATE");
366   loader->LoadTracks("RECREATE");
367  
368   //  Digits are not stored on disk and created on flight from rawdata.
369   // In order to write digits on disk the following line should be uncommented
370   // loader->LoadDigits("RECREATE"); 
371
372   //   Loop over events  
373   Int_t iEvent = 0;
374            
375   TStopwatch totalTimer;
376   TStopwatch rawTimer;
377   TStopwatch calibTimer;
378   TStopwatch clusterTimer;
379   TStopwatch trackingTimer;
380   
381   rawTimer.Start(kTRUE); rawTimer.Stop();
382   calibTimer.Start(kTRUE); calibTimer.Stop();
383   clusterTimer.Start(kTRUE); clusterTimer.Stop();
384   trackingTimer.Start(kTRUE); trackingTimer.Stop();
385   
386   totalTimer.Start(kTRUE);
387   
388   while (rawReader->NextEvent()) 
389   {
390     AliDebug(1,Form("Event %d",iEvent));
391     
392     runLoader->GetEvent(iEvent++);
393
394     //----------------------- raw2digits & raw2trigger-------------------
395 //     if (!loader->TreeD()) 
396 //     {
397 //       AliDebug(1,Form("Making Digit Container for event %d",iEvent));
398 //       loader->MakeDigitsContainer();
399 //     }
400   //  Digits are not stored on disk and created on flight from rawdata.
401   //  In order to write digits on disk the following lines should be uncommented
402   //  data.MakeBranch("D,GLT");
403   //  data.SetTreeAddress("D,GLT");
404     data.SetDataContainer("D, GLT");
405     rawTimer.Start(kFALSE);
406     fDigitMaker->Raw2Digits(rawReader);
407     rawTimer.Stop();
408     
409     if ( calibration )
410     {
411       calibTimer.Start(kFALSE);
412       calibration->ExecuteTask();
413       calibTimer.Stop();
414     }
415     // Digits are not stored on disk and created on flight from rawdata.
416     // In order to write digits on disk the following lines should be uncommented
417     // data.Fill("D,GLT");
418     // loader->WriteDigits("OVERWRITE");
419     //----------------------- digit2cluster & Trigger2Trigger -------------------
420     clusterTimer.Start(kFALSE);
421
422     if (!loader->TreeR()) loader->MakeRecPointsContainer();
423      
424     // tracking branch
425     data.MakeBranch("RC");
426     data.SetTreeAddress("RC");
427     recoCluster->Digits2Clusters(); 
428     data.Fill("RC"); 
429
430     // trigger branch
431     data.MakeBranch("TC");
432     data.SetTreeAddress("TC");
433     data.Fill("TC");
434     
435     loader->WriteRecPoints("OVERWRITE");
436
437     clusterTimer.Stop();
438
439     //---------------------------- Track & TriggerTrack ---------------------
440     trackingTimer.Start(kFALSE);
441     if (!loader->TreeT()) loader->MakeTracksContainer();
442
443     // trigger branch
444     data.MakeBranch("RL"); //trigger track
445     data.SetTreeAddress("RL");
446     recoEvent->EventReconstructTrigger();
447     data.Fill("RL");
448
449     // tracking branch
450     data.MakeBranch("RT"); //track
451     data.SetTreeAddress("RT");
452     recoEvent->EventReconstruct();
453     data.Fill("RT");
454
455     loader->WriteTracks("OVERWRITE");  
456     trackingTimer.Stop();
457     
458     //--------------------------- Resetting branches -----------------------
459     data.ResetDigits();
460     data.ResetRawClusters();
461     data.ResetTrigger();
462     data.ResetRecTracks();
463     data.ResetRecTriggerTracks();
464   
465   }
466   
467   totalTimer.Stop();
468   
469   loader->UnloadRecPoints();
470   loader->UnloadTracks();
471   loader->UnloadDigits();
472   
473   delete recoEvent;
474
475   delete recoCluster;
476   
477   AliInfo(Form("Execution time for converting RAW data to digits in MUON : R:%.2fs C:%.2fs",
478                rawTimer.RealTime(),rawTimer.CpuTime()));
479   AliInfo(Form("Execution time for calibrating MUON : R:%.2fs C:%.2fs",
480                calibTimer.RealTime(),calibTimer.CpuTime()));
481   AliInfo(Form("Execution time for clusterizing MUON : R:%.2fs C:%.2fs",
482                clusterTimer.RealTime(),clusterTimer.CpuTime()));
483   AliInfo(Form("Execution time for tracking MUON : R:%.2fs C:%.2fs",
484                trackingTimer.RealTime(),trackingTimer.CpuTime()));
485   AliInfo(Form("Total Execution time for Reconstruct(from raw) MUON : R:%.2fs C:%.2fs",
486                totalTimer.RealTime(),totalTimer.CpuTime()));
487 }
488
489 //_____________________________________________________________________________
490 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
491 {
492 /// Fill ESD
493 /// \todo add more
494
495   TClonesArray* recTracksArray = 0;
496   TClonesArray* recTrigTracksArray = 0;
497   
498   AliLoader* loader = runLoader->GetLoader("MUONLoader");
499   loader->LoadTracks("READ");
500   AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");
501
502    // declaration  
503   Int_t iEvent;// nPart;
504   Int_t nTrackHits;// nPrimary;
505   Double_t fitFmin;
506
507   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
508   Double_t xRec, yRec, zRec, chi2MatchTrigger;
509   Bool_t matchTrigger;
510
511   // setting pointer for tracks, triggertracks & trackparam at vertex
512   AliMUONTrack* recTrack = 0;
513   AliMUONTrackParam* trackParam = 0;
514   AliMUONTriggerTrack* recTriggerTrack = 0;
515
516   iEvent = runLoader->GetEventNumber(); 
517   runLoader->GetEvent(iEvent);
518
519   // Get vertex 
520   Double_t vertex[3] = {0};
521   const AliESDVertex *esdVert = esd->GetVertex(); 
522   if (esdVert->GetNContributors()) {
523     esdVert->GetXYZ(vertex);
524     printf("find vertex\n");
525   }
526   // setting ESD MUON class
527   AliESDMuonTrack* theESDTrack = new  AliESDMuonTrack() ;
528
529   //-------------------- trigger tracks-------------
530   Long_t trigPat = 0;
531   muonData->SetTreeAddress("RL");
532   muonData->GetRecTriggerTracks();
533   recTrigTracksArray = muonData->RecTriggerTracks();
534
535   // ready global trigger pattern from first track
536   if (recTrigTracksArray) 
537     recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
538   if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
539
540   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
541  
542   // -------------------- tracks-------------
543   muonData->SetTreeAddress("RT");
544   muonData->GetRecTracks();
545   recTracksArray = muonData->RecTracks();
546         
547   Int_t nRecTracks = 0;
548   if (recTracksArray)
549     nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
550   
551   // loop over tracks
552   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
553
554     // reading info from tracks
555     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
556
557     trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
558    
559     if (esdVert->GetNContributors())
560       AliMUONTrackExtrap::ExtrapToVertex(trackParam, vertex[0],vertex[1],vertex[2]);
561
562     bendingSlope            = trackParam->GetBendingSlope();
563     nonBendingSlope         = trackParam->GetNonBendingSlope();
564     inverseBendingMomentum  = trackParam->GetInverseBendingMomentum();
565     xRec  = trackParam->GetNonBendingCoor();
566     yRec  = trackParam->GetBendingCoor();
567     zRec  = trackParam->GetZ();
568
569     nTrackHits       = recTrack->GetNTrackHits();
570     fitFmin          = recTrack->GetFitFMin();
571     matchTrigger     = recTrack->GetMatchTrigger();
572     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
573
574     // setting data member of ESD MUON
575     theESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
576     theESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
577     theESDTrack->SetThetaY(TMath::ATan(bendingSlope));
578     theESDTrack->SetZ(zRec);
579     theESDTrack->SetBendingCoor(yRec); // calculate vertex at ESD or Tracking level ?
580     theESDTrack->SetNonBendingCoor(xRec);
581     theESDTrack->SetChi2(fitFmin);
582     theESDTrack->SetNHit(nTrackHits);
583     theESDTrack->SetMatchTrigger(matchTrigger);
584     theESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
585
586     // storing ESD MUON Track into ESD Event 
587     if (nRecTracks != 0)  
588       esd->AddMuonTrack(theESDTrack);
589   } // end loop tracks
590
591   // reset muondata
592   muonData->ResetRecTracks();
593   muonData->ResetRecTriggerTracks();
594
595   //} // end loop on event  
596   loader->UnloadTracks(); 
597
598   delete theESDTrack;
599   delete muonData;
600 }//_____________________________________________________________________________
601 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliRawReader* /*rawReader*/, AliESD* esd) const
602 {
603 /// Fill ESD
604 /// \todo add more
605
606   // don't need rawReader ???
607   FillESD(runLoader, esd);
608 }