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 "AliMUONEventReconstructor.h"
  34 #include "AliMUONTrack.h"
  35 #include "AliMUONTrackHit.h"
  36 #include "AliMUONTrackParam.h"
  37
  38 // Classes for Ntuple ////////////////////////////////////////////////////
  39
  40 class AliMUONTrackRecNtuple : public TObject {
  41  public:
  42   // for direct access to members
  43   Int_t fCharge; // charge of reconstructed track (+/- 1)
  44   Float_t fPxRec; // Px of reconstructed track at vertex (GeV/c)
  45   Float_t fPyRec; // Py of reconstructed track at vertex (GeV/c)
  46   Float_t fPzRec; // Pz of reconstructed track at vertex (GeV/c)
  47   Float_t fZRec; // Z of reconstructed track at vertex (cm)
  48   Float_t fZRec1; // Z of reconstructed track at first hit (cm)
  49   Int_t fNHits; // number of hits
  50   Float_t fChi2; // chi2 of fit
  51   Float_t fPxGen; // Px of best compatible generated track at vertex (GeV/c)
  52   Float_t fPyGen; // Py of best compatible generated track at vertex (GeV/c)
  53   Float_t fPzGen; // Pz of best compatible generated track at vertex (GeV/c)
  54   AliMUONTrackRecNtuple(){;} // Constructor
  55   virtual ~AliMUONTrackRecNtuple(){;} // Destructor
  56  protected:
  57  private:
  58   ClassDef(AliMUONTrackRecNtuple, 1) // AliMUONTrackRecNtuple
  59     };
  60
  61 class AliMUONHeaderRecNtuple : public TObject {
  62  public:
  63   // for direct access to members
  64   Int_t fEvent; // event number
  65   AliMUONHeaderRecNtuple(){;} // Constructor
  66   virtual ~AliMUONHeaderRecNtuple(){;} // Destructor
  67  protected:
  68  private:
  69   ClassDef(AliMUONHeaderRecNtuple, 1) // AliMUONHeaderRecNtuple
  70     };
  71
  72 ClassImp(AliMUONTrackRecNtuple) // Class implementation in ROOT context
  73 ClassImp(AliMUONHeaderRecNtuple) // Class implementation in ROOT context
  74
  75   //__________________________________________________________________________
  76 void AliMUONEventRecNtupleFill(AliMUONEventReconstructor *Reco, Int_t FillWrite = 0)
  77 {
  78   // Fill Ntuple for reconstructed event pointed to by "Reco"
  79   // if "FillWrite" different from -1.
  80   // Ntuple is created automatically before filling first event.
  81   // If "FillWrite" = -1, write and close the file
  82
  83   static Bool_t firstTime = kTRUE;
  84   static TTree *ntuple;
  85   static AliMUONHeaderRecNtuple *header = new AliMUONHeaderRecNtuple();
  86   static TClonesArray *recTracks = new TClonesArray("AliMUONTrackRecNtuple",5);
  87
  88   Int_t trackIndex;
  89   AliMUONTrack *track;
  90   AliMUONTrackParam *trackParam;
  91   AliMUONTrackRecNtuple *recTrackNt;
  92   Double_t bendingSlope, nonBendingSlope, pYZ;
  93
  94   if (FillWrite == -1) {
  95     // better to create the file before the Ntuple ????
  96     TFile *file = new TFile("MUONtrackReco.root","recreate");
  97     ntuple->Write();
  98     file->Write();
  99     file->Close();
 100   }
 101
 102   if (firstTime) {
 103     firstTime = kFALSE;
 104     // first call: create tree for Ntuple...
 105     ntuple = new TTree("MUONtrackReco", "MUONtrackReco");
 106     ntuple->Branch("Header","AliMUONHeaderRecNtuple", &header);
 107     ntuple->Branch("Tracks", &recTracks);
 108   }
 109
 110   // header
 111   header->fEvent = gAlice->GetHeader()->GetEvent();
 112
 113   TClonesArray *recoTracksPtr = Reco->GetRecTracksPtr();
 114   recoTracksPtr->Compress(); // for simple loop without "Next" since no hole
 115   recTracks->Clear(); // to reset the TClonesArray of tracks to be put in the ntuple
 116   // Loop over reconstructed tracks
 117   for (trackIndex = 0; trackIndex < Reco->GetNRecTracks(); trackIndex++) {
 118     track = (AliMUONTrack*) ((*recoTracksPtr)[trackIndex]);
 119     recTrackNt = (AliMUONTrackRecNtuple*)
 120       new ((*recTracks)[trackIndex]) AliMUONTrackRecNtuple();
 121     // track parameters at Vertex
 122     trackParam = track->GetTrackParamAtVertex();
 123     recTrackNt->fCharge =
 124       Int_t(TMath::Sign(1., trackParam->GetInverseBendingMomentum()));
 125     bendingSlope = trackParam->GetBendingSlope();
 126     nonBendingSlope = trackParam->GetNonBendingSlope();
 127     pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
 128     recTrackNt->fPzRec = pYZ / TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
 129     recTrackNt->fPxRec = recTrackNt->fPzRec * nonBendingSlope;
 130     recTrackNt->fPyRec = recTrackNt->fPzRec * bendingSlope;
 131     recTrackNt->fZRec = trackParam->GetZ();
 132     // track parameters at first hit
 133     trackParam = ((AliMUONTrackHit*)
 134                   (track->GetTrackHitsPtr()->First()))->GetTrackParam();
 135     recTrackNt->fZRec1 = trackParam->GetZ();
 136     // chi2
 137     recTrackNt->fChi2 = track->GetFitFMin();
 138     // number of hits
 139     recTrackNt->fNHits = track->GetNTrackHits();
 140     // track parameters at vertex of best compatible generated track:
 141     // in fact muon with the right charge
 142     for (int iPart = 0; iPart < gAlice->Particles()->GetEntriesFast(); iPart++) {
 143       TParticle *particle = (TParticle*) gAlice->Particles()->UncheckedAt(iPart);
 144       if ((particle->GetPdgCode() * recTrackNt->fCharge) == -13) {
 145         recTrackNt->fPxGen = particle->Px();
 146         recTrackNt->fPyGen = particle->Py();
 147         recTrackNt->fPzGen = particle->Pz();
 148       }
 149     }
 150   } // for (trackIndex = 0;...
 151
 152   ntuple->Fill();
 153
 154   return;
 155 }
 156
 157 void MUONrecoNtuple (Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t RecGeantHits = 0, Text_t *FileName = "galice.root", Text_t *BkgGeantFileName = "")
 158 {
 159   //
 160   cout << "MUON_recoNtuple" << endl;
 161   cout << "FirstEvent " << FirstEvent << endl;
 162   cout << "LastEvent " << LastEvent << endl;
 163   cout << "RecGeantHits " << RecGeantHits << endl;
 164   cout << "FileName ``" << FileName << "''" << endl;
 165   cout << "BkgGeantFileName ``" << BkgGeantFileName << "''" << endl;
 166 //   // Dynamically link some shared libs
 167 //   if (gClassTable->GetID("AliRun") < 0) {
 168 //     gROOT->LoadMacro("loadlibs.C");
 169 //     loadlibs();
 170 //   }
 171
 172   // Connect the Root Galice file containing Geometry, Kine, Hits
 173   // and eventually RawClusters
 174   TFile *file = (TFile*) gROOT->GetListOfFiles()->FindObject(FileName);
 175   if (!file) {
 176     printf("\n Creating file %s\n", FileName);
 177     file = new TFile(FileName);
 178   }
 179   else printf("\n File %s found in file list\n", FileName);
 180
 181   // Get AliRun object from file or create it if not on file
 182   if (!gAlice) {
 183     gAlice = (AliRun*) file->Get("gAlice");
 184     if (gAlice) printf("AliRun object found on file\n");
 185     if (!gAlice) {
 186       printf("\n Create new gAlice object");
 187       gAlice = new AliRun("gAlice","Alice test program");
 188     }
 189   }
 190
 191   // Initializations
 192   // AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 193   AliMUONEventReconstructor *Reco = new AliMUONEventReconstructor();
 194
 195   Reco->SetRecGeantHits(RecGeantHits);
 196
 197   // The right place for changing AliMUONEventReconstructor parameters
 198   // with respect to the default ones
 199 //   Reco->SetMaxSigma2Distance(100.0);
 200   Reco->SetPrintLevel(10);
 201 //   Reco->SetPrintLevel(1);
 202 //   Reco->SetBendingResolution(0.0);
 203 //   Reco->SetNonBendingResolution(0.0);
 204   cout << "AliMUONEventReconstructor: actual parameters" << endl;
 205   Reco->Dump();
 206 //   gObjectTable->Print();
 207
 208   // Loop over events
 209   for (Int_t event = FirstEvent; event <= LastEvent; event++) {
 210     cout << "Event: " << event << endl;
 211 //     AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 212     Int_t nparticles = gAlice->GetEvent(event);
 213     cout << "nparticles: " << nparticles << endl;
 214     // prepare background file and/or event if necessary
 215     if (RecGeantHits == 1) {
 216       if (event == FirstEvent) Reco->SetBkgGeantFile(BkgGeantFileName);
 217       if (Reco->GetBkgGeantFile())Reco->NextBkgGeantEvent();
 218     }
 219     Reco->EventReconstruct();
 220     // Dump current event
 221     Reco->EventDump();
 222     // Fill Ntuple
 223     AliMUONEventRecNtupleFill(Reco, 0);
 224 //     gObjectTable->Print();
 225
 226   } // Event loop
 227     // Write Ntuple
 228     AliMUONEventRecNtupleFill(Reco, -1);
 229 }