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