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