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