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