1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
19 // class for MUON reconstruction //
21 ///////////////////////////////////////////////////////////////////////////////
23 #include "AliMUONReconstructor.h"
26 #include "AliESDMuonTrack.h"
29 #include "AliMUONCalibrationData.h"
30 #include "AliMUONClusterFinderAZ.h"
31 #include "AliMUONClusterFinderVS.h"
32 #include "AliMUONClusterReconstructor.h"
33 #include "AliMUONData.h"
34 #include "AliMUONDigitCalibrator.h"
35 #include "AliMUONEventRecoCombi.h"
36 #include "AliMUONRawReader.h"
37 #include "AliMUONTrack.h"
38 #include "AliMUONTrackParam.h"
39 #include "AliMUONTrackReconstructor.h"
40 #include "AliMUONTriggerTrack.h"
41 #include "AliRawReader.h"
43 #include "AliRunLoader.h"
45 #include "TStopwatch.h"
47 ClassImp(AliMUONReconstructor)
49 //_____________________________________________________________________________
50 AliMUONReconstructor::AliMUONReconstructor()
51 : AliReconstructor(), fCalibrationData(0x0)
53 /// Default constructor
58 //______________________________________________________________________________
59 AliMUONReconstructor::AliMUONReconstructor(const AliMUONReconstructor& right)
60 : AliReconstructor(right)
62 /// Protected copy constructor (not implemented)
64 AliFatal("Copy constructor not provided.");
67 //_____________________________________________________________________________
68 AliMUONReconstructor::~AliMUONReconstructor()
73 delete fCalibrationData;
76 //______________________________________________________________________________
78 AliMUONReconstructor::operator=(const AliMUONReconstructor& right)
80 /// Protected assignement operator (not implemented)
82 // check assignement to self
83 if (this == &right) return *this;
85 AliFatal("Assignement operator not provided.");
90 //_____________________________________________________________________________
92 AliMUONReconstructor::GetCalibrationTask(AliMUONData* data) const
94 /// Create the calibration task(s).
96 const AliRun* run = fRunLoader->GetAliRun();
98 // Not really clean, but for the moment we must check whether the
99 // simulation has decalibrated the data or not...
100 const AliMUON* muon = static_cast<AliMUON*>(run->GetModule("MUON"));
101 if ( muon->DigitizerType().Contains("NewDigitizer") )
103 AliInfo("Calibration will occur.");
104 Int_t runNumber = run->GetRunNumber();
105 fCalibrationData = new AliMUONCalibrationData(runNumber);
106 if ( !fCalibrationData->IsValid() )
108 AliError("Could not retrieve calibrations !");
109 delete fCalibrationData;
110 fCalibrationData = 0x0;
113 TTask* calibration = new TTask("MUONCalibrator","MUON Digit calibrator");
114 calibration->Add(new AliMUONDigitCalibrator(data,fCalibrationData));
115 //FIXME: calibration->Add(something about dead channels should go here).
120 AliInfo("Detected the usage of old digitizer (w/o decalibration). "
121 "Will not calibrate then.");
126 //_____________________________________________________________________________
128 AliMUONReconstructor::Init(AliRunLoader* runLoader)
132 fRunLoader = runLoader;
135 //_____________________________________________________________________________
136 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
141 AliLoader* loader = runLoader->GetLoader("MUONLoader");
142 Int_t nEvents = runLoader->GetNumberOfEvents();
144 AliMUONData* data = new AliMUONData(loader,"MUON","MUON");
146 // passing loader as argument.
147 AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader, data);
149 if (strstr(GetOption(),"Original"))
150 recoEvent->SetTrackMethod(1); // Original tracking
151 else if (strstr(GetOption(),"Combi"))
152 recoEvent->SetTrackMethod(3); // Combined cluster / track
154 recoEvent->SetTrackMethod(2); // Kalman
156 AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(data);
158 AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
160 if (!strstr(GetOption(),"VS")) {
161 recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
162 recoCluster->SetRecoModel(recModel);
164 recModel->SetGhostChi2Cut(10);
166 loader->LoadDigits("READ");
167 loader->LoadRecPoints("RECREATE");
168 loader->LoadTracks("RECREATE");
170 TTask* calibration = GetCalibrationTask(data);
172 Int_t chBeg = recoEvent->GetTrackMethod() == 3 ? 6 : 0;
174 for(Int_t ievent = 0; ievent < nEvents; ievent++) {
176 AliDebug(1,Form("Event %d",ievent));
178 runLoader->GetEvent(ievent);
180 //----------------------- digit2cluster & Trigger2Trigger -------------------
181 if (!loader->TreeR()) loader->MakeRecPointsContainer();
184 if (recoEvent->GetTrackMethod() != 3) {
185 data->MakeBranch("RC");
186 data->SetTreeAddress("D,RC");
188 data->SetTreeAddress("D");
189 data->SetTreeAddress("RCC");
191 // Important for avoiding a memory leak when reading digits ( to be investigated more in detail)
192 // In any case the reading of GLT is needed for the Trigger2Tigger method below
193 data->SetTreeAddress("GLT");
199 calibration->ExecuteTask();
202 recoCluster->Digits2Clusters(chBeg);
204 if (recoEvent->GetTrackMethod() == 3) {
205 // Combined cluster / track finder
206 AliMUONEventRecoCombi::Instance()->FillEvent(data, (AliMUONClusterFinderAZ*)recModel);
207 ((AliMUONClusterFinderAZ*) recModel)->SetReco(2);
209 else data->Fill("RC");
212 data->MakeBranch("TC");
213 data->SetTreeAddress("TC");
214 recoCluster->Trigger2Trigger();
217 //AZ loader->WriteRecPoints("OVERWRITE");
219 //---------------------------- Track & TriggerTrack ---------------------
220 if (!loader->TreeT()) loader->MakeTracksContainer();
223 data->MakeBranch("RL"); //trigger track
224 data->SetTreeAddress("RL");
225 recoEvent->EventReconstructTrigger();
229 data->MakeBranch("RT"); //track
230 data->SetTreeAddress("RT");
231 recoEvent->EventReconstruct();
234 loader->WriteTracks("OVERWRITE");
236 if (recoEvent->GetTrackMethod() == 3) {
237 // Combined cluster / track
238 ((AliMUONClusterFinderAZ*) recModel)->SetReco(1);
239 data->MakeBranch("RC");
240 data->SetTreeAddress("RC");
241 AliMUONEventRecoCombi::Instance()->FillRecP(data, recoEvent);
244 loader->WriteRecPoints("OVERWRITE");
246 //--------------------------- Resetting branches -----------------------
248 data->ResetRawClusters();
249 data->ResetTrigger();
251 data->ResetRawClusters();
252 data->ResetTrigger();
253 data->ResetRecTracks();
254 data->ResetRecTriggerTracks();
257 loader->UnloadDigits();
258 loader->UnloadRecPoints();
259 loader->UnloadTracks();
267 //_____________________________________________________________________________
268 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawReader) const
274 AliLoader* loader = runLoader->GetLoader("MUONLoader");
275 AliMUONData data(loader,"MUON","MUON");
277 // passing loader as argument.
278 AliMUONTrackReconstructor recoEvent(loader, &data);
280 AliMUONRawReader rawData(&data);
282 AliMUONClusterReconstructor recoCluster(&data);
284 if (strstr(GetOption(),"Original"))
286 recoEvent.SetTrackMethod(1); // Original tracking
290 recoEvent.SetTrackMethod(2); // Kalman
293 AliMUONClusterFinderVS *recModel = recoCluster.GetRecoModel();
294 if (!strstr(GetOption(),"VS"))
296 recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
297 recoCluster.SetRecoModel(recModel);
299 recModel->SetGhostChi2Cut(10);
301 TTask* calibration = GetCalibrationTask(&data);
303 loader->LoadRecPoints("RECREATE");
304 loader->LoadTracks("RECREATE");
305 loader->LoadDigits("READ");
310 TStopwatch totalTimer;
312 TStopwatch calibTimer;
313 TStopwatch clusterTimer;
314 TStopwatch trackingTimer;
316 rawTimer.Start(kTRUE); rawTimer.Stop();
317 calibTimer.Start(kTRUE); calibTimer.Stop();
318 clusterTimer.Start(kTRUE); clusterTimer.Stop();
319 trackingTimer.Start(kTRUE); trackingTimer.Stop();
321 totalTimer.Start(kTRUE);
323 while (rawReader->NextEvent())
325 AliDebug(1,Form("Event %d",iEvent));
327 runLoader->GetEvent(iEvent++);
329 //----------------------- raw2digits & raw2trigger-------------------
330 if (!loader->TreeD())
332 AliDebug(1,Form("Making Digit Container for event %d",iEvent));
333 loader->MakeDigitsContainer();
336 data.SetTreeAddress("D,GLT");
337 rawTimer.Start(kFALSE);
338 rawData.Raw2Digits(rawReader);
343 calibTimer.Start(kFALSE);
344 calibration->ExecuteTask();
348 //----------------------- digit2cluster & Trigger2Trigger -------------------
349 clusterTimer.Start(kFALSE);
351 if (!loader->TreeR()) loader->MakeRecPointsContainer();
354 data.MakeBranch("RC");
355 data.SetTreeAddress("RC");
356 recoCluster.Digits2Clusters();
360 data.MakeBranch("TC");
361 data.SetTreeAddress("TC");
362 // recoCluster.Trigger2Trigger();
365 loader->WriteRecPoints("OVERWRITE");
369 //---------------------------- Track & TriggerTrack ---------------------
370 trackingTimer.Start(kFALSE);
371 if (!loader->TreeT()) loader->MakeTracksContainer();
374 data.MakeBranch("RL"); //trigger track
375 data.SetTreeAddress("RL");
376 recoEvent.EventReconstructTrigger();
380 data.MakeBranch("RT"); //track
381 data.SetTreeAddress("RT");
382 recoEvent.EventReconstruct();
385 loader->WriteTracks("OVERWRITE");
386 trackingTimer.Stop();
388 //--------------------------- Resetting branches -----------------------
390 data.ResetRawClusters();
393 data.ResetRawClusters();
395 data.ResetRecTracks();
396 data.ResetRecTriggerTracks();
402 loader->UnloadRecPoints();
403 loader->UnloadTracks();
404 loader->UnloadDigits();
406 AliInfo(Form("Execution time for converting RAW data to digits in MUON : R:%.2fs C:%.2fs",
407 rawTimer.RealTime(),rawTimer.CpuTime()));
408 AliInfo(Form("Execution time for calibrating MUON : R:%.2fs C:%.2fs",
409 calibTimer.RealTime(),calibTimer.CpuTime()));
410 AliInfo(Form("Execution time for clusterizing MUON : R:%.2fs C:%.2fs",
411 clusterTimer.RealTime(),clusterTimer.CpuTime()));
412 AliInfo(Form("Execution time for tracking MUON : R:%.2fs C:%.2fs",
413 trackingTimer.RealTime(),trackingTimer.CpuTime()));
414 AliInfo(Form("Total Execution time for Reconstruct(from raw) MUON : R:%.2fs C:%.2fs",
415 totalTimer.RealTime(),totalTimer.CpuTime()));
418 //_____________________________________________________________________________
419 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
424 TClonesArray* recTracksArray = 0;
425 TClonesArray* recTrigTracksArray = 0;
427 AliLoader* loader = runLoader->GetLoader("MUONLoader");
428 loader->LoadTracks("READ");
429 AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");
432 Int_t iEvent;// nPart;
433 Int_t nTrackHits;// nPrimary;
436 Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
437 Double_t xRec, yRec, zRec, chi2MatchTrigger;
440 // setting pointer for tracks, triggertracks & trackparam at vertex
441 AliMUONTrack* recTrack = 0;
442 AliMUONTrackParam* trackParam = 0;
443 AliMUONTriggerTrack* recTriggerTrack = 0;
444 // TParticle* particle = new TParticle();
445 // AliGenEventHeader* header = 0;
446 iEvent = runLoader->GetEventNumber();
447 runLoader->GetEvent(iEvent);
450 Double_t vertex[3] = {0};
451 const AliESDVertex *esdVert = esd->GetVertex();
452 if (esdVert) esdVert->GetXYZ(vertex);
455 // if ( (header = runLoader->GetHeader()->GenEventHeader()) ) {
456 // header->PrimaryVertex(vertex);
458 // runLoader->LoadKinematics("READ");
459 // runLoader->TreeK()->GetBranch("Particles")->SetAddress(&particle);
460 // nPart = (Int_t)runLoader->TreeK()->GetEntries();
461 // for(Int_t iPart = 0; iPart < nPart; iPart++) {
462 // runLoader->TreeK()->GetEvent(iPart);
463 // if (particle->GetFirstMother() == -1) {
464 // vertex[0] += particle->Vx();
465 // vertex[1] += particle->Vy();
466 // vertex[2] += particle->Vz();
470 // vertex[0] /= (double)nPrimary;
471 // vertex[1] /= (double)nPrimary;
472 // vertex[2] /= (double)nPrimary;
476 // setting ESD MUON class
477 AliESDMuonTrack* theESDTrack = new AliESDMuonTrack() ;
479 //-------------------- trigger tracks-------------
481 muonData->SetTreeAddress("RL");
482 muonData->GetRecTriggerTracks();
483 recTrigTracksArray = muonData->RecTriggerTracks();
485 // ready global trigger pattern from first track
486 if (recTrigTracksArray)
487 recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
488 if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
490 //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
492 // -------------------- tracks-------------
493 muonData->SetTreeAddress("RT");
494 muonData->GetRecTracks();
495 recTracksArray = muonData->RecTracks();
497 Int_t nRecTracks = 0;
499 nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
502 for (Int_t iRecTracks = 0; iRecTracks < nRecTracks; iRecTracks++) {
504 // reading info from tracks
505 recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
507 trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
508 trackParam->ExtrapToVertex(vertex[0],vertex[1],vertex[2]);
510 bendingSlope = trackParam->GetBendingSlope();
511 nonBendingSlope = trackParam->GetNonBendingSlope();
512 inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
513 xRec = trackParam->GetNonBendingCoor();
514 yRec = trackParam->GetBendingCoor();
515 zRec = trackParam->GetZ();
517 nTrackHits = recTrack->GetNTrackHits();
518 fitFmin = recTrack->GetFitFMin();
519 matchTrigger = recTrack->GetMatchTrigger();
520 chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
522 // setting data member of ESD MUON
523 theESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
524 theESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
525 theESDTrack->SetThetaY(TMath::ATan(bendingSlope));
526 theESDTrack->SetZ(zRec);
527 theESDTrack->SetBendingCoor(yRec); // calculate vertex at ESD or Tracking level ?
528 theESDTrack->SetNonBendingCoor(xRec);
529 theESDTrack->SetChi2(fitFmin);
530 theESDTrack->SetNHit(nTrackHits);
531 theESDTrack->SetMatchTrigger(matchTrigger);
532 theESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
534 // storing ESD MUON Track into ESD Event
536 esd->AddMuonTrack(theESDTrack);
539 // add global trigger pattern
541 esd->SetTrigger(trigPat);
544 muonData->ResetRecTracks();
545 muonData->ResetRecTriggerTracks();
547 //} // end loop on event
548 loader->UnloadTracks();
550 // runLoader->UnloadKinematics();
554 }//_____________________________________________________________________________
555 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliRawReader* /*rawReader*/, AliESD* esd) const
560 // don't need rawReader ???
561 FillESD(runLoader, esd);