MUON/MUONrecoNtuple.C

   1 // Macro MUONrecoNtuple.C (TO BE COMPILED)
   2 // for testing the C++ reconstruction code
   3 // and producing an Ntuple with reconstructed tracks
   4 // in output file "MUONtrackReco.root".
   5 // An example for using the Ntuple is in the macro MUONmassPlot.C
   6
   7 // Arguments:
   8 //   FirstEvent (default 0)
   9 //   LastEvent (default 0)
  10 //   RecGeantHits (1 to reconstruct GEANT hits) (default 0)
  11 //   FileName (for signal) (default "galice.root")
  12 //   BkgGeantFileName (for background),
  13 //      needed only if RecGeantHits = 1 and background to be added
  14
  15 // IMPORTANT NOTICE FOR USERS:
  16 // under "root" or "root.exe", execute the following commands:
  17 // 1. "gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER -I$ROOTSYS/include")" to get the right path at compilation time
  18 // 2. ".x loadlibs.C" to load the shared libraries
  19 // 3. ".x MUONrecoNtuple.C+()" with the right arguments, without forgetting the "+" which implies the compilation of the macro before its execution
  20
  21 #include <iostream.h>
  22
  23 #include <TClassTable.h>
  24 #include <TClonesArray.h>
  25 #include <TFile.h>
  26 #include <TParticle.h>
  27 #include <TROOT.h>
  28 #include <TTree.h>
  29
  30 #include "AliHeader.h"
  31 #include "AliRun.h"
  32
  33 #include "AliMUON.h"
  34 #include "AliMUONEventReconstructor.h"
  35 #include "AliMUONTrack.h"
  36 #include "AliMUONTrackHit.h"
  37 #include "AliMUONTrackParam.h"
  38
  39 // Classes for Ntuple ////////////////////////////////////////////////////
  40
  41 class AliMUONTrackRecNtuple : public TObject {
  42  public:
  43   // for direct access to members
  44   Int_t fCharge; // charge of reconstructed track (+/- 1)
  45   Float_t fPxRec; // Px of reconstructed track at vertex (GeV/c)
  46   Float_t fPyRec; // Py of reconstructed track at vertex (GeV/c)
  47   Float_t fPzRec; // Pz of reconstructed track at vertex (GeV/c)
  48   Float_t fZRec; // Z of reconstructed track at vertex (cm)
  49   Float_t fZRec1; // Z of reconstructed track at first hit (cm)
  50   Int_t fNHits; // number of hits
  51   Float_t fChi2; // chi2 of fit
  52   Float_t fPxGen; // Px of best compatible generated track at vertex (GeV/c)
  53   Float_t fPyGen; // Py of best compatible generated track at vertex (GeV/c)
  54   Float_t fPzGen; // Pz of best compatible generated track at vertex (GeV/c)
  55   AliMUONTrackRecNtuple(){;} // Constructor
  56   virtual ~AliMUONTrackRecNtuple(){;} // Destructor
  57  protected:
  58  private:
  59   ClassDef(AliMUONTrackRecNtuple, 1) // AliMUONTrackRecNtuple
  60     };
  61
  62 class AliMUONHeaderRecNtuple : public TObject {
  63  public:
  64   // for direct access to members
  65   Int_t fEvent; // event number
  66   AliMUONHeaderRecNtuple(){;} // Constructor
  67   virtual ~AliMUONHeaderRecNtuple(){;} // Destructor
  68  protected:
  69  private:
  70   ClassDef(AliMUONHeaderRecNtuple, 1) // AliMUONHeaderRecNtuple
  71     };
  72
  73 ClassImp(AliMUONTrackRecNtuple) // Class implementation in ROOT context
  74 ClassImp(AliMUONHeaderRecNtuple) // Class implementation in ROOT context
  75
  76   //__________________________________________________________________________
  77 void AliMUONEventRecNtupleFill(AliMUONEventReconstructor *Reco, Int_t FillWrite = 0)
  78 {
  79   // Fill Ntuple for reconstructed event pointed to by "Reco"
  80   // if "FillWrite" different from -1.
  81   // Ntuple is created automatically before filling first event.
  82   // If "FillWrite" = -1, write and close the file
  83
  84   static Bool_t firstTime = kTRUE;
  85   static TTree *ntuple;
  86   static AliMUONHeaderRecNtuple *header = new AliMUONHeaderRecNtuple();
  87   static TClonesArray *recTracks = new TClonesArray("AliMUONTrackRecNtuple",5);
  88
  89   Int_t trackIndex;
  90   AliMUONTrack *track;
  91   AliMUONTrackParam *trackParam;
  92   AliMUONTrackRecNtuple *recTrackNt;
  93   Double_t bendingSlope, nonBendingSlope, pYZ;
  94
  95   if (FillWrite == -1) {
  96     printf(">>> Writing Ntuple of reconstructed tracks\n");
  97     // better to create the file before the Ntuple ????
  98     TFile *file = new TFile("MUONtrackReco.root","recreate");
  99     ntuple->Write();
 100     file->Write();
 101     file->Close();
 102   }
 103
 104   if (firstTime) {
 105     firstTime = kFALSE;
 106     // first call: create tree for Ntuple...
 107     printf(">>> Creating Ntuple of reconstructed tracks\n");
 108     ntuple = new TTree("MUONtrackReco", "MUONtrackReco");
 109     ntuple->Branch("Header","AliMUONHeaderRecNtuple", &header);
 110     ntuple->Branch("Tracks", &recTracks);
 111   }
 112
 113   // header
 114
 115   header->fEvent = gAlice->GetHeader()->GetEvent();
 116
 117   TClonesArray *recoTracksPtr = Reco->GetRecTracksPtr();
 118   recoTracksPtr->Compress(); // for simple loop without "Next" since no hole
 119   recTracks->Clear(); // to reset the TClonesArray of tracks to be put in the ntuple
 120   // Loop over reconstructed tracks
 121   for (trackIndex = 0; trackIndex < Reco->GetNRecTracks(); trackIndex++) {
 122     track = (AliMUONTrack*) ((*recoTracksPtr)[trackIndex]);
 123     recTrackNt = (AliMUONTrackRecNtuple*)
 124       new ((*recTracks)[trackIndex]) AliMUONTrackRecNtuple();
 125     // track parameters at Vertex
 126     trackParam = track->GetTrackParamAtVertex();
 127     recTrackNt->fCharge =
 128       Int_t(TMath::Sign(1., trackParam->GetInverseBendingMomentum()));
 129     bendingSlope = trackParam->GetBendingSlope();
 130     nonBendingSlope = trackParam->GetNonBendingSlope();
 131     pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
 132     recTrackNt->fPzRec = pYZ / TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
 133     recTrackNt->fPxRec = recTrackNt->fPzRec * nonBendingSlope;
 134     recTrackNt->fPyRec = recTrackNt->fPzRec * bendingSlope;
 135     recTrackNt->fZRec = trackParam->GetZ();
 136     // track parameters at first hit
 137     trackParam = ((AliMUONTrackHit*)
 138                   (track->GetTrackHitsPtr()->First()))->GetTrackParam();
 139     recTrackNt->fZRec1 = trackParam->GetZ();
 140     // chi2
 141     recTrackNt->fChi2 = track->GetFitFMin();
 142     // number of hits
 143     recTrackNt->fNHits = track->GetNTrackHits();
 144     printf("test> Px %f Py %f Pz %f \n", recTrackNt->fPxRec,  recTrackNt->fPyRec,  recTrackNt->fPzRec);
 145     // track parameters at vertex of best compatible generated track:
 146     // in fact muon with the right charge
 147     TTree* mtreeK=gAlice->TreeK();
 148     TBranch *brparticle = mtreeK->GetBranch("Particles");
 149     Int_t nPart = brparticle->GetEntries();
 150     TParticle *particle = new TParticle();
 151     mtreeK->SetBranchAddress("Particles",&particle);
 152     for (Int_t iPart = 0; iPart < nPart; iPart++) {
 153        brparticle->GetEntry(iPart);
 154        //cout << "Code Particle: " << particle->GetPdgCode() << "\n";
 155        if ((particle->GetPdgCode() * recTrackNt->fCharge) == -13) {
 156         recTrackNt->fPxGen = particle->Px();
 157         recTrackNt->fPyGen = particle->Py();
 158         recTrackNt->fPzGen = particle->Pz();
 159         printf("Gen: Px %f Py %f Pz %f \n", recTrackNt->fPxGen, recTrackNt->fPyGen, recTrackNt->fPzGen);
 160        }
 161    }
 162   } // for (trackIndex = 0;...
 163
 164   printf(">>> Filling Ntuple of reconstructed tracks\n");
 165   ntuple->Fill();
 166
 167   return;
 168 }
 169
 170 void MUONrecoNtuple (Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t RecGeantHits = 0, Text_t *FileName = "galice.root", Text_t *BkgGeantFileName = "")
 171 {
 172   //
 173   cout << "MUON_recoNtuple" << endl;
 174   cout << "FirstEvent " << FirstEvent << endl;
 175   cout << "LastEvent " << LastEvent << endl;
 176   cout << "RecGeantHits " << RecGeantHits << endl;
 177   cout << "FileName ``" << FileName << "''" << endl;
 178   cout << "BkgGeantFileName ``" << BkgGeantFileName << "''" << endl;
 179 //   // Dynamically link some shared libs
 180 //   if (gClassTable->GetID("AliRun") < 0) {
 181 //     gROOT->LoadMacro("loadlibs.C");
 182 //     loadlibs();
 183 //   }
 184
 185
 186   // Creating Run Loader and openning file containing Hits, Digits and RecPoints
 187   AliRunLoader * RunLoader = AliRunLoader::Open(FileName,"Event","UPDATE");
 188   if (RunLoader ==0x0) {
 189     printf(">>> Error : Error Opening %s file \n",FileName);
 190     return;
 191   }
 192   // Loading AliRun master
 193   RunLoader->LoadgAlice();
 194   gAlice = RunLoader->GetAliRun();
 195   RunLoader->LoadKinematics("READ");
 196
 197   // Loading MUON subsystem
 198   AliMUON * MUON = (AliMUON *) gAlice->GetDetector("MUON");
 199   AliLoader * MUONLoader = RunLoader->GetLoader("MUONLoader");
 200   MUONLoader->LoadHits("READ");
 201   MUONLoader->LoadRecPoints("READ");
 202
 203   Int_t ievent, nevents;
 204   nevents = RunLoader->GetNumberOfEvents();
 205
 206   // Initializations
 207   // AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 208   MUON->SetTreeAddress();
 209   AliMUONEventReconstructor *Reco = new AliMUONEventReconstructor();
 210
 211   Reco->SetRecGeantHits(RecGeantHits);
 212
 213   // The right place for changing AliMUONEventReconstructor parameters
 214   // with respect to the default ones
 215 //   Reco->SetMaxSigma2Distance(100.0);
 216 //  Reco->SetPrintLevel(20);
 217    Reco->SetPrintLevel(1);
 218 //   Reco->SetBendingResolution(0.0);
 219 //   Reco->SetNonBendingResolution(0.0);
 220   cout << "AliMUONEventReconstructor: actual parameters" << endl;
 221   Reco->Dump();
 222 //   gObjectTable->Print();
 223   // Loop over events
 224   if (LastEvent>nevents) LastEvent = nevents;
 225   for (Int_t event = FirstEvent; event < LastEvent; event++) {
 226     cout << "Event: " << event << endl;
 227     RunLoader->GetEvent(event);
 228     //     Int_t nparticles = gAlice->GetEvent(event);
 229     //      cout << "nparticles: " << nparticles << endl;
 230     // prepare background file and/or event if necessary
 231     if (RecGeantHits == 1) {
 232       if (event == FirstEvent) Reco->SetBkgGeantFile(BkgGeantFileName);
 233       if (Reco->GetBkgGeantFile())Reco->NextBkgGeantEvent();
 234     }
 235     Reco->EventReconstruct();
 236     // Dump current event
 237     Reco->EventDump();
 238     // Fill Ntuple
 239     AliMUONEventRecNtupleFill(Reco, 0);
 240 //     gObjectTable->Print();
 241     MUON->ResetRawClusters();
 242   } // Event loop
 243     // Write Ntuple
 244     AliMUONEventRecNtupleFill(Reco, -1);
 245 }