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