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