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 #include <TParticle.h>
  23 #include <TArrayF.h>
  24
  25 #include "AliRunLoader.h"
  26 #include "AliHeader.h"
  27 #include "AliGenEventHeader.h"
  28 #include "AliESD.h"
  29 #include "AliMUONReconstructor.h"
  30
  31 #include "AliMUONData.h"
  32 #include "AliMUONTrackReconstructor.h"
  33 #include "AliMUONClusterReconstructor.h"
  34 #include "AliMUONClusterFinderVS.h"
  35 #include "AliMUONClusterFinderAZ.h"
  36 #include "AliMUONEventRecoCombi.h"
  37 #include "AliMUONTrack.h"
  38 #include "AliMUONTrackParam.h"
  39 #include "AliMUONTriggerTrack.h"
  40 #include "AliESDMuonTrack.h"
  41 #include "AliMUONRawData.h"
  42
  43 #include "AliRawReader.h"
  44
  45
  46 ClassImp(AliMUONReconstructor)
  47 //_____________________________________________________________________________
  48 AliMUONReconstructor::AliMUONReconstructor()
  49   : AliReconstructor()
  50 {
  51 }
  52 //_____________________________________________________________________________
  53 AliMUONReconstructor::~AliMUONReconstructor()
  54 {
  55 }
  56 //_____________________________________________________________________________
  57 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
  58 {
  59 //  AliLoader
  60   AliLoader* loader = runLoader->GetLoader("MUONLoader");
  61   Int_t nEvents = runLoader->GetNumberOfEvents();
  62
  63 // used local container for each method
  64 // passing fLoader as argument, could be avoided ???
  65   AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader);
  66   AliMUONData* dataEvent = recoEvent->GetMUONData();
  67   if (strstr(GetOption(),"Original")) recoEvent->SetTrackMethod(1); // Original tracking
  68   else if (strstr(GetOption(),"Combi")) recoEvent->SetTrackMethod(3); // Combined cluster / track
  69   else recoEvent->SetTrackMethod(2); // Kalman
  70
  71   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader);
  72   AliMUONData* dataCluster = recoCluster->GetMUONData();
  73   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
  74   if (!strstr(GetOption(),"VS")) {
  75     recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
  76     recoCluster->SetRecoModel(recModel);
  77   }
  78   recModel->SetGhostChi2Cut(10);
  79
  80   loader->LoadDigits("READ");
  81   loader->LoadRecPoints("RECREATE");
  82   loader->LoadTracks("RECREATE");
  83
  84   Int_t chBeg = recoEvent->GetTrackMethod() == 3 ? 6 : 0;
  85   //   Loop over events
  86   for(Int_t ievent = 0; ievent < nEvents; ievent++) {
  87     printf("Event %d\n",ievent);
  88     runLoader->GetEvent(ievent);
  89
  90     //----------------------- digit2cluster & Trigger2Trigger -------------------
  91     if (!loader->TreeR()) loader->MakeRecPointsContainer();
  92
  93     // tracking branch
  94     if (recoEvent->GetTrackMethod() != 3) {
  95       dataCluster->MakeBranch("RC");
  96       dataCluster->SetTreeAddress("D,RC");
  97     } else {
  98       dataCluster->SetTreeAddress("D");
  99       dataCluster->SetTreeAddress("RCC");
 100     }
 101     // Important for avoiding a memory leak when reading digits ( to be investigated more in detail)
 102     // In any case the reading of GLT is needed for the Trigger2Tigger method below
 103     dataCluster->SetTreeAddress("GLT");
 104
 105     recoCluster->Digits2Clusters(chBeg);
 106     if (recoEvent->GetTrackMethod() == 3) {
 107       // Combined cluster / track finder
 108       AliMUONEventRecoCombi::Instance()->FillEvent(dataCluster, dataEvent, (AliMUONClusterFinderAZ*)recModel);
 109       ((AliMUONClusterFinderAZ*) recModel)->SetReco(2);
 110     }
 111     else dataCluster->Fill("RC");
 112
 113     // trigger branch
 114     dataCluster->MakeBranch("TC");
 115     dataCluster->SetTreeAddress("TC");
 116     recoCluster->Trigger2Trigger();
 117     dataCluster->Fill("TC");
 118
 119     //AZ loader->WriteRecPoints("OVERWRITE");
 120
 121     //---------------------------- Track & TriggerTrack ---------------------
 122     if (!loader->TreeT()) loader->MakeTracksContainer();
 123
 124     // trigger branch
 125     dataEvent->MakeBranch("RL"); //trigger track
 126     dataEvent->SetTreeAddress("RL");
 127     recoEvent->EventReconstructTrigger();
 128     dataEvent->Fill("RL");
 129
 130     // tracking branch
 131     dataEvent->MakeBranch("RT"); //track
 132     dataEvent->SetTreeAddress("RT");
 133     recoEvent->EventReconstruct();
 134     dataEvent->Fill("RT");
 135
 136     loader->WriteTracks("OVERWRITE");
 137
 138     if (recoEvent->GetTrackMethod() == 3) {
 139       // Combined cluster / track
 140       ((AliMUONClusterFinderAZ*) recModel)->SetReco(1);
 141       dataCluster->MakeBranch("RC");
 142       dataCluster->SetTreeAddress("RC");
 143       AliMUONEventRecoCombi::Instance()->FillRecP(dataCluster, recoEvent);
 144       dataCluster->Fill("RC");
 145     }
 146     loader->WriteRecPoints("OVERWRITE");
 147
 148     //--------------------------- Resetting branches -----------------------
 149     dataCluster->ResetDigits();
 150     dataCluster->ResetRawClusters();
 151     dataCluster->ResetTrigger();
 152
 153     dataEvent->ResetRawClusters();
 154     dataEvent->ResetTrigger();
 155     dataEvent->ResetRecTracks();
 156     dataEvent->ResetRecTriggerTracks();
 157
 158   }
 159   loader->UnloadDigits();
 160   loader->UnloadRecPoints();
 161   loader->UnloadTracks();
 162
 163   delete recoCluster;
 164   delete recoEvent;
 165 }
 166
 167 //_____________________________________________________________________________
 168 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawReader) const
 169 {
 170 //  AliLoader
 171   AliLoader* loader = runLoader->GetLoader("MUONLoader");
 172
 173 // used local container for each method
 174 // passing fLoader as argument, could be avoided ???
 175   AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader);
 176   AliMUONData* dataEvent = recoEvent->GetMUONData();
 177
 178   AliMUONRawData* rawData = new AliMUONRawData(loader);
 179   AliMUONData* dataCluster = rawData->GetMUONData();
 180
 181   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader, dataCluster);
 182   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
 183   recModel->SetGhostChi2Cut(10);
 184
 185   loader->LoadRecPoints("RECREATE");
 186   loader->LoadTracks("RECREATE");
 187   loader->LoadDigits("RECREATE");
 188
 189
 190   //   Loop over events
 191   Int_t iEvent = 0;
 192
 193   while (rawReader->NextEvent()) {
 194     printf("Event %d\n",iEvent);
 195     runLoader->GetEvent(iEvent++);
 196
 197     //----------------------- raw2digits & raw2trigger-------------------
 198     if (!loader->TreeD()) loader->MakeDigitsContainer();
 199
 200     // tracking branch
 201     dataCluster->MakeBranch("D");
 202     dataCluster->SetTreeAddress("D");
 203     rawData->ReadTrackerDDL(rawReader);
 204     dataCluster->Fill("D");
 205
 206     // trigger branch
 207     dataCluster->MakeBranch("GLT");
 208     dataCluster->SetTreeAddress("GLT");
 209     rawData->ReadTriggerDDL(rawReader);
 210     dataCluster->Fill("GLT");
 211
 212     loader->WriteDigits("OVERWRITE");
 213
 214     //----------------------- digit2cluster & Trigger2Trigger -------------------
 215     if (!loader->TreeR()) loader->MakeRecPointsContainer();
 216
 217     // tracking branch
 218     dataCluster->MakeBranch("RC");
 219     dataCluster->SetTreeAddress("RC");
 220     recoCluster->Digits2Clusters();
 221     dataCluster->Fill("RC");
 222
 223     // trigger branch
 224     dataCluster->MakeBranch("TC");
 225     dataCluster->SetTreeAddress("TC");
 226     recoCluster->Trigger2Trigger();
 227     dataCluster->Fill("TC");
 228
 229     loader->WriteRecPoints("OVERWRITE");
 230
 231     //---------------------------- Track & TriggerTrack ---------------------
 232     if (!loader->TreeT()) loader->MakeTracksContainer();
 233
 234     // trigger branch
 235     dataEvent->MakeBranch("RL"); //trigger track
 236     dataEvent->SetTreeAddress("RL");
 237     recoEvent->EventReconstructTrigger();
 238     dataEvent->Fill("RL");
 239
 240     // tracking branch
 241     dataEvent->MakeBranch("RT"); //track
 242     dataEvent->SetTreeAddress("RT");
 243     recoEvent->EventReconstruct();
 244     dataEvent->Fill("RT");
 245
 246     loader->WriteTracks("OVERWRITE");
 247
 248     //--------------------------- Resetting branches -----------------------
 249     dataCluster->ResetDigits();
 250     dataCluster->ResetRawClusters();
 251     dataCluster->ResetTrigger();
 252
 253     dataEvent->ResetRawClusters();
 254     dataEvent->ResetTrigger();
 255     dataEvent->ResetRecTracks();
 256     dataEvent->ResetRecTriggerTracks();
 257
 258   }
 259   loader->UnloadRecPoints();
 260   loader->UnloadTracks();
 261   loader->UnloadDigits();
 262
 263   delete recoCluster;
 264   delete recoEvent;
 265 }
 266
 267 //_____________________________________________________________________________
 268 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
 269 {
 270   TClonesArray* recTracksArray = 0;
 271   TClonesArray* recTrigTracksArray = 0;
 272
 273   AliLoader* loader = runLoader->GetLoader("MUONLoader");
 274   loader->LoadTracks("READ");
 275   AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");
 276
 277    // declaration
 278   Int_t iEvent;// nPart;
 279   Int_t nTrackHits;// nPrimary;
 280   Double_t fitFmin;
 281   TArrayF vertex(3);
 282
 283   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
 284   Double_t xRec, yRec, zRec, chi2MatchTrigger;
 285   Bool_t matchTrigger;
 286
 287   // setting pointer for tracks, triggertracks & trackparam at vertex
 288   AliMUONTrack* recTrack = 0;
 289   AliMUONTrackParam* trackParam = 0;
 290   AliMUONTriggerTrack* recTriggerTrack = 0;
 291 //   TParticle* particle = new TParticle();
 292 //   AliGenEventHeader* header = 0;
 293   iEvent = runLoader->GetEventNumber();
 294   runLoader->GetEvent(iEvent);
 295
 296   // vertex calculation (maybe it exists already somewhere else)
 297   vertex[0] = vertex[1] = vertex[2] = 0.;
 298  //  nPrimary = 0;
 299 //   if ( (header = runLoader->GetHeader()->GenEventHeader()) ) {
 300 //     header->PrimaryVertex(vertex);
 301 //   } else {
 302 //     runLoader->LoadKinematics("READ");
 303 //     runLoader->TreeK()->GetBranch("Particles")->SetAddress(&particle);
 304 //     nPart = (Int_t)runLoader->TreeK()->GetEntries();
 305 //     for(Int_t iPart = 0; iPart < nPart; iPart++) {
 306 //       runLoader->TreeK()->GetEvent(iPart);
 307 //       if (particle->GetFirstMother() == -1) {
 308 //      vertex[0] += particle->Vx();
 309 //      vertex[1] += particle->Vy();
 310 //      vertex[2] += particle->Vz();
 311 //      nPrimary++;
 312 //       }
 313 //       if (nPrimary) {
 314 //      vertex[0] /= (double)nPrimary;
 315 //      vertex[1] /= (double)nPrimary;
 316 //      vertex[2] /= (double)nPrimary;
 317 //       }
 318 //     }
 319 //   }
 320   // setting ESD MUON class
 321   AliESDMuonTrack* theESDTrack = new  AliESDMuonTrack() ;
 322
 323   //-------------------- trigger tracks-------------
 324   Long_t trigPat = 0;
 325   muonData->SetTreeAddress("RL");
 326   muonData->GetRecTriggerTracks();
 327   recTrigTracksArray = muonData->RecTriggerTracks();
 328
 329   // ready global trigger pattern from first track
 330   if (recTrigTracksArray)
 331     recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
 332   if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
 333
 334   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
 335
 336   // -------------------- tracks-------------
 337   muonData->SetTreeAddress("RT");
 338   muonData->GetRecTracks();
 339   recTracksArray = muonData->RecTracks();
 340
 341   Int_t nRecTracks = 0;
 342   if (recTracksArray)
 343     nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
 344
 345   // loop over tracks
 346   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
 347
 348     // reading info from tracks
 349     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
 350
 351     trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
 352     trackParam->ExtrapToVertex(vertex[0],vertex[1],vertex[2]);
 353
 354     bendingSlope            = trackParam->GetBendingSlope();
 355     nonBendingSlope         = trackParam->GetNonBendingSlope();
 356     inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
 357     xRec  = trackParam->GetNonBendingCoor();
 358     yRec  = trackParam->GetBendingCoor();
 359     zRec  = trackParam->GetZ();
 360
 361     nTrackHits       = recTrack->GetNTrackHits();
 362     fitFmin          = recTrack->GetFitFMin();
 363     matchTrigger     = recTrack->GetMatchTrigger();
 364     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
 365
 366     // setting data member of ESD MUON
 367     theESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
 368     theESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
 369     theESDTrack->SetThetaY(TMath::ATan(bendingSlope));
 370     theESDTrack->SetZ(zRec);
 371     theESDTrack->SetBendingCoor(yRec); // calculate vertex at ESD or Tracking level ?
 372     theESDTrack->SetNonBendingCoor(xRec);
 373     theESDTrack->SetChi2(fitFmin);
 374     theESDTrack->SetNHit(nTrackHits);
 375     theESDTrack->SetMatchTrigger(matchTrigger);
 376     theESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
 377
 378     // storing ESD MUON Track into ESD Event
 379     if (nRecTracks != 0)
 380       esd->AddMuonTrack(theESDTrack);
 381   } // end loop tracks
 382
 383   // add global trigger pattern
 384   if (nRecTracks != 0)
 385     esd->SetTrigger(trigPat);
 386
 387   // reset muondata
 388   muonData->ResetRecTracks();
 389   muonData->ResetRecTriggerTracks();
 390
 391   //} // end loop on event
 392   loader->UnloadTracks();
 393  //  if (!header)
 394 //     runLoader->UnloadKinematics();
 395   delete theESDTrack;
 396   delete muonData;
 397   // delete particle;
 398 }//_____________________________________________________________________________
 399 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliRawReader* /*rawReader*/, AliESD* esd) const
 400 {
 401   // don't need rawReader ???
 402   FillESD(runLoader, esd);
 403 }