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 "AliGenEventHeader.h"
  28 #include "AliHeader.h"
  29 #include "AliMUON.h"
  30 #include "AliMUONCalibrationData.h"
  31 #include "AliMUONClusterFinderAZ.h"
  32 #include "AliMUONClusterFinderVS.h"
  33 #include "AliMUONClusterReconstructor.h"
  34 #include "AliMUONData.h"
  35 #include "AliMUONDigitCalibrator.h"
  36 #include "AliMUONEventRecoCombi.h"
  37 #include "AliMUONRawReader.h"
  38 #include "AliMUONTrack.h"
  39 #include "AliMUONTrackParam.h"
  40 #include "AliMUONTrackReconstructor.h"
  41 #include "AliMUONTriggerTrack.h"
  42 #include "AliRawReader.h"
  43 #include "AliRun.h"
  44 #include "AliRunLoader.h"
  45 #include "TTask.h"
  46 #include <TArrayF.h>
  47 #include <TParticle.h>
  48
  49 ClassImp(AliMUONReconstructor)
  50
  51 //_____________________________________________________________________________
  52 AliMUONReconstructor::AliMUONReconstructor()
  53   : AliReconstructor(), fCalibrationData(0x0)
  54 {
  55 }
  56 //_____________________________________________________________________________
  57 AliMUONReconstructor::~AliMUONReconstructor()
  58 {
  59   delete fCalibrationData;
  60 }
  61
  62 //_____________________________________________________________________________
  63 TTask*
  64 AliMUONReconstructor::GetCalibrationTask(AliMUONData* data) const
  65 {
  66   //
  67   // Create the calibration task(s).
  68   //
  69
  70   const AliRun* run = fRunLoader->GetAliRun();
  71
  72   // Not really clean, but for the moment we must check whether the
  73   // simulation has decalibrated the data or not...
  74   const AliMUON* muon = static_cast<AliMUON*>(run->GetModule("MUON"));
  75   if ( muon->DigitizerType().Contains("NewDigitizer") )
  76   {
  77     AliInfo("Calibration will occur.");
  78     Int_t runNumber = run->GetRunNumber();
  79     fCalibrationData = new AliMUONCalibrationData(runNumber);
  80     if ( !fCalibrationData->IsValid() )
  81     {
  82       AliError("Could not retrieve calibrations !");
  83       delete fCalibrationData;
  84       fCalibrationData = 0x0;
  85       return 0x0;
  86     }
  87     TTask* calibration = new TTask("MUONCalibrator","MUON Digit calibrator");
  88     calibration->Add(new AliMUONDigitCalibrator(data,fCalibrationData));
  89     //FIXME: calibration->Add(something about dead channels should go here).
  90     return calibration;
  91   }
  92   else
  93   {
  94     AliInfo("Detected the usage of old digitizer (w/o decalibration). "
  95             "Will not calibrate then.");
  96     return 0x0;
  97   }
  98 }
  99
 100 //_____________________________________________________________________________
 101 void
 102 AliMUONReconstructor::Init(AliRunLoader* runLoader)
 103 {
 104   fRunLoader = runLoader;
 105 }
 106
 107 //_____________________________________________________________________________
 108 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
 109 {
 110 //  AliLoader
 111   AliLoader* loader = runLoader->GetLoader("MUONLoader");
 112   Int_t nEvents = runLoader->GetNumberOfEvents();
 113
 114   AliMUONData* data = new AliMUONData(loader,"MUON","MUON");
 115
 116 // passing loader as argument.
 117   AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader, data);
 118
 119   if (strstr(GetOption(),"Original"))
 120     recoEvent->SetTrackMethod(1); // Original tracking
 121   else if (strstr(GetOption(),"Combi"))
 122     recoEvent->SetTrackMethod(3); // Combined cluster / track
 123   else
 124     recoEvent->SetTrackMethod(2); // Kalman
 125
 126   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader, data);
 127
 128   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
 129
 130   if (!strstr(GetOption(),"VS")) {
 131     recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
 132     recoCluster->SetRecoModel(recModel);
 133   }
 134   recModel->SetGhostChi2Cut(10);
 135
 136   loader->LoadDigits("READ");
 137   loader->LoadRecPoints("RECREATE");
 138   loader->LoadTracks("RECREATE");
 139
 140   TTask* calibration = GetCalibrationTask(data);
 141
 142   Int_t chBeg = recoEvent->GetTrackMethod() == 3 ? 6 : 0;
 143   //   Loop over events
 144   for(Int_t ievent = 0; ievent < nEvents; ievent++) {
 145     printf("Event %d\n",ievent);
 146     runLoader->GetEvent(ievent);
 147
 148     //----------------------- digit2cluster & Trigger2Trigger -------------------
 149     if (!loader->TreeR()) loader->MakeRecPointsContainer();
 150
 151     // tracking branch
 152     if (recoEvent->GetTrackMethod() != 3) {
 153       data->MakeBranch("RC");
 154       data->SetTreeAddress("D,RC");
 155     } else {
 156       data->SetTreeAddress("D");
 157       data->SetTreeAddress("RCC");
 158     }
 159     // Important for avoiding a memory leak when reading digits ( to be investigated more in detail)
 160     // In any case the reading of GLT is needed for the Trigger2Tigger method below
 161     data->SetTreeAddress("GLT");
 162
 163     data->GetDigits();
 164
 165     if ( calibration )
 166     {
 167       calibration->ExecuteTask();
 168     }
 169
 170     recoCluster->Digits2Clusters(chBeg);
 171
 172     if (recoEvent->GetTrackMethod() == 3) {
 173       // Combined cluster / track finder
 174       AliMUONEventRecoCombi::Instance()->FillEvent(data, (AliMUONClusterFinderAZ*)recModel);
 175       ((AliMUONClusterFinderAZ*) recModel)->SetReco(2);
 176     }
 177     else data->Fill("RC");
 178
 179     // trigger branch
 180     data->MakeBranch("TC");
 181     data->SetTreeAddress("TC");
 182     recoCluster->Trigger2Trigger();
 183     data->Fill("TC");
 184
 185     //AZ loader->WriteRecPoints("OVERWRITE");
 186
 187     //---------------------------- Track & TriggerTrack ---------------------
 188     if (!loader->TreeT()) loader->MakeTracksContainer();
 189
 190     // trigger branch
 191     data->MakeBranch("RL"); //trigger track
 192     data->SetTreeAddress("RL");
 193     recoEvent->EventReconstructTrigger();
 194     data->Fill("RL");
 195
 196     // tracking branch
 197     data->MakeBranch("RT"); //track
 198     data->SetTreeAddress("RT");
 199     recoEvent->EventReconstruct();
 200     data->Fill("RT");
 201
 202     loader->WriteTracks("OVERWRITE");
 203
 204     if (recoEvent->GetTrackMethod() == 3) {
 205       // Combined cluster / track
 206       ((AliMUONClusterFinderAZ*) recModel)->SetReco(1);
 207       data->MakeBranch("RC");
 208       data->SetTreeAddress("RC");
 209       AliMUONEventRecoCombi::Instance()->FillRecP(data, recoEvent);
 210       data->Fill("RC");
 211     }
 212     loader->WriteRecPoints("OVERWRITE");
 213
 214     //--------------------------- Resetting branches -----------------------
 215     data->ResetDigits();
 216     data->ResetRawClusters();
 217     data->ResetTrigger();
 218
 219     data->ResetRawClusters();
 220     data->ResetTrigger();
 221     data->ResetRecTracks();
 222     data->ResetRecTriggerTracks();
 223
 224   }
 225   loader->UnloadDigits();
 226   loader->UnloadRecPoints();
 227   loader->UnloadTracks();
 228
 229   delete recoCluster;
 230   delete recoEvent;
 231   delete data;
 232   delete calibration;
 233 }
 234
 235 //_____________________________________________________________________________
 236 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawReader) const
 237 {
 238 //  AliLoader
 239   AliLoader* loader = runLoader->GetLoader("MUONLoader");
 240   AliMUONData* data = new AliMUONData(loader,"MUON","MUON");
 241
 242 // passing loader as argument.
 243   AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader, data);
 244
 245   AliMUONRawReader* rawData = new AliMUONRawReader(loader, data);
 246
 247   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader, data);
 248
 249   if (strstr(GetOption(),"Original"))
 250     recoEvent->SetTrackMethod(1); // Original tracking
 251 //   else if (strstr(GetOption(),"Combi"))
 252 //     recoEvent->SetTrackMethod(3); // Combined cluster / track
 253   else
 254     recoEvent->SetTrackMethod(2); // Kalman
 255
 256   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
 257   if (!strstr(GetOption(),"VS")) {
 258     recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
 259     recoCluster->SetRecoModel(recModel);
 260   }
 261   recModel->SetGhostChi2Cut(10);
 262
 263   TTask* calibration = GetCalibrationTask(data);
 264
 265   loader->LoadRecPoints("RECREATE");
 266   loader->LoadTracks("RECREATE");
 267   loader->LoadDigits("RECREATE");
 268
 269   //   Loop over events
 270   Int_t iEvent = 0;
 271
 272   while (rawReader->NextEvent()) {
 273     printf("Event %d\n",iEvent);
 274     runLoader->GetEvent(iEvent++);
 275
 276     //----------------------- raw2digits & raw2trigger-------------------
 277     if (!loader->TreeD()) loader->MakeDigitsContainer();
 278
 279     // tracking branch
 280     data->MakeBranch("D");
 281     data->SetTreeAddress("D");
 282     rawData->ReadTrackerDDL(rawReader);
 283     data->Fill("D");
 284
 285     // trigger branch
 286     data->MakeBranch("GLT");
 287     data->SetTreeAddress("GLT");
 288     rawData->ReadTriggerDDL(rawReader);
 289     data->Fill("GLT");
 290
 291     if ( calibration )
 292     {
 293       calibration->ExecuteTask();
 294     }
 295
 296     loader->WriteDigits("OVERWRITE");
 297
 298     //----------------------- digit2cluster & Trigger2Trigger -------------------
 299     if (!loader->TreeR()) loader->MakeRecPointsContainer();
 300
 301     // tracking branch
 302     data->MakeBranch("RC");
 303     data->SetTreeAddress("RC");
 304     recoCluster->Digits2Clusters();
 305     data->Fill("RC");
 306
 307     // trigger branch
 308     data->MakeBranch("TC");
 309     data->SetTreeAddress("TC");
 310     recoCluster->Trigger2Trigger();
 311     data->Fill("TC");
 312
 313     loader->WriteRecPoints("OVERWRITE");
 314
 315     //---------------------------- Track & TriggerTrack ---------------------
 316     if (!loader->TreeT()) loader->MakeTracksContainer();
 317
 318     // trigger branch
 319     data->MakeBranch("RL"); //trigger track
 320     data->SetTreeAddress("RL");
 321     recoEvent->EventReconstructTrigger();
 322     data->Fill("RL");
 323
 324     // tracking branch
 325     data->MakeBranch("RT"); //track
 326     data->SetTreeAddress("RT");
 327     recoEvent->EventReconstruct();
 328     data->Fill("RT");
 329
 330     loader->WriteTracks("OVERWRITE");
 331
 332     //--------------------------- Resetting branches -----------------------
 333     data->ResetDigits();
 334     data->ResetRawClusters();
 335     data->ResetTrigger();
 336
 337     data->ResetRawClusters();
 338     data->ResetTrigger();
 339     data->ResetRecTracks();
 340     data->ResetRecTriggerTracks();
 341
 342   }
 343   loader->UnloadRecPoints();
 344   loader->UnloadTracks();
 345   loader->UnloadDigits();
 346
 347   delete recoCluster;
 348   delete recoEvent;
 349   delete data;
 350 }
 351
 352 //_____________________________________________________________________________
 353 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
 354 {
 355   TClonesArray* recTracksArray = 0;
 356   TClonesArray* recTrigTracksArray = 0;
 357
 358   AliLoader* loader = runLoader->GetLoader("MUONLoader");
 359   loader->LoadTracks("READ");
 360   AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");
 361
 362    // declaration
 363   Int_t iEvent;// nPart;
 364   Int_t nTrackHits;// nPrimary;
 365   Double_t fitFmin;
 366   TArrayF vertex(3);
 367
 368   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
 369   Double_t xRec, yRec, zRec, chi2MatchTrigger;
 370   Bool_t matchTrigger;
 371
 372   // setting pointer for tracks, triggertracks & trackparam at vertex
 373   AliMUONTrack* recTrack = 0;
 374   AliMUONTrackParam* trackParam = 0;
 375   AliMUONTriggerTrack* recTriggerTrack = 0;
 376 //   TParticle* particle = new TParticle();
 377 //   AliGenEventHeader* header = 0;
 378   iEvent = runLoader->GetEventNumber();
 379   runLoader->GetEvent(iEvent);
 380
 381   // vertex calculation (maybe it exists already somewhere else)
 382   vertex[0] = vertex[1] = vertex[2] = 0.;
 383  //  nPrimary = 0;
 384 //   if ( (header = runLoader->GetHeader()->GenEventHeader()) ) {
 385 //     header->PrimaryVertex(vertex);
 386 //   } else {
 387 //     runLoader->LoadKinematics("READ");
 388 //     runLoader->TreeK()->GetBranch("Particles")->SetAddress(&particle);
 389 //     nPart = (Int_t)runLoader->TreeK()->GetEntries();
 390 //     for(Int_t iPart = 0; iPart < nPart; iPart++) {
 391 //       runLoader->TreeK()->GetEvent(iPart);
 392 //       if (particle->GetFirstMother() == -1) {
 393 //      vertex[0] += particle->Vx();
 394 //      vertex[1] += particle->Vy();
 395 //      vertex[2] += particle->Vz();
 396 //      nPrimary++;
 397 //       }
 398 //       if (nPrimary) {
 399 //      vertex[0] /= (double)nPrimary;
 400 //      vertex[1] /= (double)nPrimary;
 401 //      vertex[2] /= (double)nPrimary;
 402 //       }
 403 //     }
 404 //   }
 405   // setting ESD MUON class
 406   AliESDMuonTrack* theESDTrack = new  AliESDMuonTrack() ;
 407
 408   //-------------------- trigger tracks-------------
 409   Long_t trigPat = 0;
 410   muonData->SetTreeAddress("RL");
 411   muonData->GetRecTriggerTracks();
 412   recTrigTracksArray = muonData->RecTriggerTracks();
 413
 414   // ready global trigger pattern from first track
 415   if (recTrigTracksArray)
 416     recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
 417   if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
 418
 419   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
 420
 421   // -------------------- tracks-------------
 422   muonData->SetTreeAddress("RT");
 423   muonData->GetRecTracks();
 424   recTracksArray = muonData->RecTracks();
 425
 426   Int_t nRecTracks = 0;
 427   if (recTracksArray)
 428     nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
 429
 430   // loop over tracks
 431   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
 432
 433     // reading info from tracks
 434     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
 435
 436     trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
 437     trackParam->ExtrapToVertex(vertex[0],vertex[1],vertex[2]);
 438
 439     bendingSlope            = trackParam->GetBendingSlope();
 440     nonBendingSlope         = trackParam->GetNonBendingSlope();
 441     inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
 442     xRec  = trackParam->GetNonBendingCoor();
 443     yRec  = trackParam->GetBendingCoor();
 444     zRec  = trackParam->GetZ();
 445
 446     nTrackHits       = recTrack->GetNTrackHits();
 447     fitFmin          = recTrack->GetFitFMin();
 448     matchTrigger     = recTrack->GetMatchTrigger();
 449     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
 450
 451     // setting data member of ESD MUON
 452     theESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
 453     theESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
 454     theESDTrack->SetThetaY(TMath::ATan(bendingSlope));
 455     theESDTrack->SetZ(zRec);
 456     theESDTrack->SetBendingCoor(yRec); // calculate vertex at ESD or Tracking level ?
 457     theESDTrack->SetNonBendingCoor(xRec);
 458     theESDTrack->SetChi2(fitFmin);
 459     theESDTrack->SetNHit(nTrackHits);
 460     theESDTrack->SetMatchTrigger(matchTrigger);
 461     theESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
 462
 463     // storing ESD MUON Track into ESD Event
 464     if (nRecTracks != 0)
 465       esd->AddMuonTrack(theESDTrack);
 466   } // end loop tracks
 467
 468   // add global trigger pattern
 469   if (nRecTracks != 0)
 470     esd->SetTrigger(trigPat);
 471
 472   // reset muondata
 473   muonData->ResetRecTracks();
 474   muonData->ResetRecTriggerTracks();
 475
 476   //} // end loop on event
 477   loader->UnloadTracks();
 478  //  if (!header)
 479 //     runLoader->UnloadKinematics();
 480   delete theESDTrack;
 481   delete muonData;
 482   // delete particle;
 483 }//_____________________________________________________________________________
 484 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliRawReader* /*rawReader*/, AliESD* esd) const
 485 {
 486   // don't need rawReader ???
 487   FillESD(runLoader, esd);
 488 }