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