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