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