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