]> git.uio.no Git - u/mrichter/AliRoot.git/blob - MUON/AliMUONReconstructor.cxx
Check we get everything from OCDB even before starting, and correctly deleting calibr...
[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   // Check that we get all the calibrations we'll need
140   if ( !fCalibrationData->Pedestals() ||
141        !fCalibrationData->Gains() ||
142        !fCalibrationData->HV() )
143   {
144     AliFatal("Could not access all required calibration data");
145   }
146   TTask* calibration = new TTask("MUONCalibrator","MUON Digit calibrator");
147   
148   TString opt(GetOption());
149   opt.ToUpper();
150   Bool_t statusMap(kTRUE);
151   
152   if ( strstr(opt,"NOSTATUSMAP") )
153   {
154     AliWarning("Disconnecting status map : SHOULD BE USED FOR DEBUG ONLY. NOT FOR PRODUCTION !!!");
155     statusMap = kFALSE; 
156   }
157   calibration->Add(new AliMUONDigitCalibrator(fMUONData,fCalibrationData,statusMap));
158   return calibration;
159 }
160
161 //_____________________________________________________________________________
162 AliMUONClusterReconstructor*
163 AliMUONReconstructor::CreateClusterReconstructor() const
164 {
165 /// Create cluster reconstructor
166
167   AliMUONVClusterFinder* clusterFinder(0x0);
168   
169   TString opt(GetOption());
170   opt.ToUpper();
171   
172   if ( strstr(opt,"PRECLUSTER") )
173   {
174     clusterFinder = new AliMUONPreClusterFinder;
175   }  
176   else if ( strstr(opt,"COG") )
177   {
178     clusterFinder = new AliMUONClusterFinderCOG;
179   }  
180   else if ( strstr(opt,"SIMPLEFIT") )
181   {
182     clusterFinder = new AliMUONClusterFinderSimpleFit;
183   }
184   else if ( strstr(opt,"MLEM:DRAW") )
185   {
186     clusterFinder = new AliMUONClusterFinderMLEM(kTRUE);
187   }
188   else if ( strstr(opt,"MLEM") )
189   {
190     clusterFinder = new AliMUONClusterFinderMLEM(kFALSE);
191   } 
192   
193   if ( clusterFinder) 
194   {
195     AliInfo(Form("Will use %s for clusterizing",clusterFinder->ClassName()));
196   }
197   
198   AliMUONClusterReconstructor* clusterReco = 
199     new AliMUONClusterReconstructor(fMUONData,clusterFinder,fTransformer);
200   return clusterReco;
201 }
202
203 //_____________________________________________________________________________
204 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
205 {
206 /// Reconstruct
207 /// \todo add more
208
209   AliLoader* loader = runLoader->GetLoader("MUONLoader");
210   Int_t nEvents     = runLoader->GetNumberOfEvents();
211   Int_t evtNumber   = runLoader->GetEventNumber();
212
213   fMUONData->SetLoader(loader);
214
215 // passing loader as argument.
216   AliMUONVTrackReconstructor* recoEvent;
217   if (strstr(GetOption(),"Original")) recoEvent = new AliMUONTrackReconstructor(fMUONData);
218   else if (strstr(GetOption(),"Combi")) recoEvent = new AliMUONTrackReconstructorK(fMUONData,"Combi");
219   else recoEvent = new AliMUONTrackReconstructorK(fMUONData,"Kalman");
220   
221   recoEvent->SetTriggerCircuit(fTriggerCircuit);
222
223   AliMUONClusterReconstructor* recoCluster = CreateClusterReconstructor();
224   
225   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
226
227   if (!strstr(GetOption(),"VS")) {
228     recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
229     recoCluster->SetRecoModel(recModel);
230   }
231   recModel->SetGhostChi2Cut(10);
232   recModel->SetEventNumber(evtNumber);
233
234   loader->LoadDigits("READ");
235   loader->LoadRecPoints("RECREATE");
236   loader->LoadTracks("RECREATE");
237   
238   TTask* calibration = GetCalibrationTask();
239   
240   Int_t chBeg = (strstr(GetOption(),"Combi") ? 6 : 0);
241
242   //   Loop over events              
243   for(Int_t ievent = 0; ievent < nEvents; ievent++) {
244
245     AliDebug(1,Form("Event %d",ievent));
246     
247     runLoader->GetEvent(ievent);
248
249     //----------------------- digit2cluster & Trigger2Trigger -------------------
250     if (!loader->TreeR()) loader->MakeRecPointsContainer();
251      
252     // tracking branch
253     if (!strstr(GetOption(),"Combi")) {
254       fMUONData->MakeBranch("RC");
255       fMUONData->SetTreeAddress("D,RC");
256     } else {
257       fMUONData->SetTreeAddress("D");
258       fMUONData->SetTreeAddress("RCC");
259     }
260     // Important for avoiding a memory leak when reading digits ( to be investigated more in detail)
261     // In any case the reading of GLT is needed for the Trigger2Tigger method below
262     fMUONData->SetTreeAddress("GLT");
263
264     fMUONData->GetDigits();
265     
266     if ( calibration ) 
267     {
268       calibration->ExecuteTask();
269     }
270     
271     recoCluster->Digits2Clusters(chBeg); 
272     
273     if (strstr(GetOption(),"Combi")) {
274       // Combined cluster / track finder
275       AliMUONEventRecoCombi::Instance()->FillEvent(fMUONData, (AliMUONClusterFinderAZ*)recModel);
276       ((AliMUONClusterFinderAZ*) recModel)->SetReco(2); 
277     }
278     else fMUONData->Fill("RC"); 
279
280     // trigger branch
281     fMUONData->MakeBranch("TC");
282     fMUONData->SetTreeAddress("TC");
283     recoCluster->Trigger2Trigger(); 
284     fMUONData->Fill("TC");
285
286     //AZ loader->WriteRecPoints("OVERWRITE");
287
288     //---------------------------- Track & TriggerTrack ---------------------
289     if (!loader->TreeT()) loader->MakeTracksContainer();
290
291     // trigger branch
292     fMUONData->MakeBranch("RL"); //trigger track
293     fMUONData->SetTreeAddress("RL");
294     recoEvent->EventReconstructTrigger();
295     fMUONData->Fill("RL");
296
297     // tracking branch
298     fMUONData->MakeBranch("RT"); //track
299     fMUONData->SetTreeAddress("RT");
300     recoEvent->EventReconstruct();
301     fMUONData->Fill("RT");
302
303     loader->WriteTracks("OVERWRITE"); 
304   
305     if (strstr(GetOption(),"Combi")) { 
306       // Combined cluster / track
307       ((AliMUONClusterFinderAZ*) recModel)->SetReco(1);
308       fMUONData->MakeBranch("RC");
309       fMUONData->SetTreeAddress("RC");
310       AliMUONEventRecoCombi::Instance()->FillRecP(fMUONData, (AliMUONTrackReconstructorK*)recoEvent); 
311       fMUONData->Fill("RC"); 
312     }
313     loader->WriteRecPoints("OVERWRITE"); 
314
315     //--------------------------- Resetting branches -----------------------
316     fMUONData->ResetDigits();
317     fMUONData->ResetRawClusters();
318     fMUONData->ResetTrigger();
319     fMUONData->ResetRecTracks();  
320     fMUONData->ResetRecTriggerTracks();
321
322   }
323   loader->UnloadDigits();
324   loader->UnloadRecPoints();
325   loader->UnloadTracks();
326
327   delete recoCluster;
328   delete recoEvent;
329   delete calibration;
330 }
331
332 //_____________________________________________________________________________
333 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, 
334                                        AliRawReader* rawReader) const
335 {
336 /// Recontruct
337 /// \todo add more
338
339   //  AliLoader
340   AliLoader* loader = runLoader->GetLoader("MUONLoader");
341   Int_t evtNumber   = runLoader->GetEventNumber();
342  
343   fMUONData->SetLoader(loader);
344
345   // passing loader as argument.
346   fDigitMaker->SetMUONData(fMUONData);
347
348   // disable trigger rawdata reading
349   if (strstr(GetOption(),"TriggerDisable"))
350       fDigitMaker->DisableTrigger();
351
352   AliMUONClusterReconstructor* recoCluster = CreateClusterReconstructor();
353
354   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
355
356   if (!strstr(GetOption(),"VS")) 
357   {
358     recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
359     recoCluster->SetRecoModel(recModel);
360   }
361   recModel->SetGhostChi2Cut(10);
362   recModel->SetEventNumber(evtNumber);
363   
364   TTask* calibration = GetCalibrationTask();
365   
366   loader->LoadRecPoints("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   
380   rawTimer.Start(kTRUE); rawTimer.Stop();
381   calibTimer.Start(kTRUE); calibTimer.Stop();
382   clusterTimer.Start(kTRUE); clusterTimer.Stop();
383   
384   totalTimer.Start(kTRUE);
385   
386   while (rawReader->NextEvent()) 
387   {
388     AliDebug(1,Form("Event %d",iEvent));
389     
390     runLoader->GetEvent(iEvent++);
391     
392     //----------------------- raw2digits & raw2trigger-------------------
393     //  if (!loader->TreeD()) 
394     //  {
395     //    AliDebug(1,Form("Making Digit Container for event %d",iEvent));
396     //    loader->MakeDigitsContainer();
397     //  }
398     //  Digits are not stored on disk and created on flight from rawdata.
399     //  In order to write digits on disk the following lines should be uncommented
400     //  fMUONData->MakeBranch("D,GLT");
401     //  fMUONData->SetTreeAddress("D,GLT");
402     fMUONData->SetDataContainer("D, GLT");
403     rawTimer.Start(kFALSE);
404     fDigitMaker->Raw2Digits(rawReader);
405     rawTimer.Stop();
406     
407     if ( calibration )
408     {
409       calibTimer.Start(kFALSE);
410       calibration->ExecuteTask();
411       calibTimer.Stop();
412     }
413     // Digits are not stored on disk and created on flight from rawdata.
414     // In order to write digits on disk the following lines should be uncommented
415     // fMUONData->Fill("D,GLT");
416     // loader->WriteDigits("OVERWRITE");
417     //----------------------- digit2cluster & Trigger2Trigger -------------------
418     clusterTimer.Start(kFALSE);
419
420     if (!loader->TreeR()) loader->MakeRecPointsContainer();
421     
422     // tracking branch
423     fMUONData->MakeBranch("RC");
424     fMUONData->SetTreeAddress("RC");
425     recoCluster->Digits2Clusters(); 
426     fMUONData->Fill("RC"); 
427
428     // trigger branch
429     fMUONData->MakeBranch("TC");
430     fMUONData->SetTreeAddress("TC");
431     fMUONData->Fill("TC");
432     
433     loader->WriteRecPoints("OVERWRITE");
434
435     clusterTimer.Stop();
436     
437     
438     //--------------------------- Resetting branches -----------------------
439     fMUONData->ResetDigits();
440     fMUONData->ResetRawClusters();
441     fMUONData->ResetTrigger();
442   }
443   
444   totalTimer.Stop();
445   
446   loader->UnloadRecPoints();
447   loader->UnloadDigits();
448   
449   delete recoCluster;
450   delete calibration;
451   
452   AliInfo(Form("Execution time for converting RAW data to digits in MUON : R:%.2fs C:%.2fs",
453                rawTimer.RealTime(),rawTimer.CpuTime()));
454   AliInfo(Form("Execution time for calibrating MUON : R:%.2fs C:%.2fs",
455                calibTimer.RealTime(),calibTimer.CpuTime()));
456   AliInfo(Form("Execution time for clusterizing MUON : R:%.2fs C:%.2fs",
457                clusterTimer.RealTime(),clusterTimer.CpuTime()));
458   AliInfo(Form("Total Execution time for Reconstruct(from raw) MUON : R:%.2fs C:%.2fs",
459                totalTimer.RealTime(),totalTimer.CpuTime()));
460 }
461
462 //_____________________________________________________________________________
463 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
464 {
465 /// Fill ESD
466 /// \todo add more
467
468   TClonesArray* recTracksArray = 0;
469   TClonesArray* recTrigTracksArray = 0;
470   
471   AliLoader* loader = runLoader->GetLoader("MUONLoader");
472   loader->LoadTracks("READ");
473   fMUONData->SetLoader(loader);
474
475    // declaration  
476   Int_t iEvent;// nPart;
477   Int_t nTrackHits;// nPrimary;
478   Double_t fitFmin;
479
480   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
481   Double_t xRec, yRec, zRec, chi2MatchTrigger;
482   Bool_t matchTrigger;
483
484   // setting pointer for tracks, triggertracks & trackparam at vertex
485   AliMUONTrack* recTrack = 0;
486   AliMUONTrackParam trackParam;
487   AliMUONTriggerTrack* recTriggerTrack = 0;
488
489   iEvent = runLoader->GetEventNumber(); 
490   runLoader->GetEvent(iEvent);
491
492   // Get vertex 
493   Double_t vertex[3] = {0};
494   const AliESDVertex *esdVert = esd->GetVertex(); 
495   if (esdVert->GetNContributors()) {
496     esdVert->GetXYZ(vertex);
497     AliDebug(1, "find vertex\n");
498   }
499   // setting ESD MUON class
500   AliESDMuonTrack* theESDTrack = new  AliESDMuonTrack() ;
501
502   //-------------------- trigger tracks-------------
503   Long_t trigPat = 0;
504   fMUONData->SetTreeAddress("RL");
505   fMUONData->GetRecTriggerTracks();
506   recTrigTracksArray = fMUONData->RecTriggerTracks();
507
508   // ready global trigger pattern from first track
509   if (recTrigTracksArray) 
510     recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
511   if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
512
513   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
514  
515   // -------------------- tracks-------------
516   fMUONData->SetTreeAddress("RT");
517   fMUONData->GetRecTracks();
518   recTracksArray = fMUONData->RecTracks();
519         
520   Int_t nRecTracks = 0;
521   if (recTracksArray)
522     nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
523   
524   // loop over tracks
525   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
526
527     // reading info from tracks
528     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
529
530     trackParam = *((AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First());
531    
532     // extrapolate to the vertex if required
533     //   if the vertex is not available, extrapolate to (0,0,0)
534     if (!strstr(GetOption(),"NoExtrapToVtx")) {
535       if (esdVert->GetNContributors())
536         AliMUONTrackExtrap::ExtrapToVertex(&trackParam, vertex[0],vertex[1],vertex[2]);
537       else
538         AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0.,0.,0.);
539     }
540     
541     bendingSlope            = trackParam.GetBendingSlope();
542     nonBendingSlope         = trackParam.GetNonBendingSlope();
543     inverseBendingMomentum  = trackParam.GetInverseBendingMomentum();
544     xRec  = trackParam.GetNonBendingCoor();
545     yRec  = trackParam.GetBendingCoor();
546     zRec  = trackParam.GetZ();
547
548     nTrackHits       = recTrack->GetNTrackHits();
549     fitFmin          = recTrack->GetFitFMin();
550     matchTrigger     = recTrack->GetMatchTrigger();
551     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
552
553     // setting data member of ESD MUON
554     theESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
555     theESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
556     theESDTrack->SetThetaY(TMath::ATan(bendingSlope));
557     theESDTrack->SetZ(zRec);
558     theESDTrack->SetBendingCoor(yRec); // calculate vertex at ESD or Tracking level ?
559     theESDTrack->SetNonBendingCoor(xRec);
560     theESDTrack->SetChi2(fitFmin);
561     theESDTrack->SetNHit(nTrackHits);
562     theESDTrack->SetMatchTrigger(matchTrigger);
563     theESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
564
565     // storing ESD MUON Track into ESD Event 
566     if (nRecTracks != 0)  
567       esd->AddMuonTrack(theESDTrack);
568   } // end loop tracks
569
570   // reset muondata
571   fMUONData->ResetRecTracks();
572   fMUONData->ResetRecTriggerTracks();
573
574   //} // end loop on event  
575   loader->UnloadTracks(); 
576
577   delete theESDTrack;
578 }
579
580 //_____________________________________________________________________________
581 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliRawReader* /*rawReader*/, AliESD* esd) const
582 {
583 /// Fill ESD
584 /// \todo add more
585
586   // don't need rawReader ???
587   FillESD(runLoader, esd);
588 }
589
590 //_____________________________________________________________________________
591 AliTracker* AliMUONReconstructor::CreateTracker(AliRunLoader* runLoader) const
592 {
593   /// create tracker for MUON
594   /// go into the tracking framework finally (Ch.F)
595  
596   AliLoader* loader = runLoader->GetLoader("MUONLoader");
597
598   fMUONData->SetLoader(loader);
599
600   AliMUONTracker* tracker = new AliMUONTracker();
601   tracker->SetMUONData(fMUONData);
602   tracker->SetTriggerCircuit(fTriggerCircuit);
603   tracker->SetOption(GetOption());
604
605   return tracker;
606
607 }