MUON/MUONrecoNtuple.C

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