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