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