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