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