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