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