Transition to NewIO
[u/mrichter/AliRoot.git] / 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     // track parameters at vertex of best compatible generated track:
145     // in fact muon with the right charge
146  //     for (int iPart = 0; iPart < gAlice->Particles()->GetEntriesFast(); iPart++) {
147 //        TParticle *particle = (TParticle*) gAlice->Particles()->UncheckedAt(iPart);
148 //        if ((particle->GetPdgCode() * recTrackNt->fCharge) == -13) {
149 //      recTrackNt->fPxGen = particle->Px();
150 //      recTrackNt->fPyGen = particle->Py();
151 //      recTrackNt->fPzGen = particle->Pz();
152 //        }  // commented by Gines
153     //    }
154   } // for (trackIndex = 0;...
155
156   printf(">>> Filling Ntuple of reconstructed tracks\n");
157   ntuple->Fill();
158
159   return;
160 }
161
162 void MUONrecoNtuple (Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t RecGeantHits = 0, Text_t *FileName = "galice.root", Text_t *BkgGeantFileName = "")
163 {
164   //
165   cout << "MUON_recoNtuple" << endl;
166   cout << "FirstEvent " << FirstEvent << endl;
167   cout << "LastEvent " << LastEvent << endl;
168   cout << "RecGeantHits " << RecGeantHits << endl;
169   cout << "FileName ``" << FileName << "''" << endl;
170   cout << "BkgGeantFileName ``" << BkgGeantFileName << "''" << endl;
171 //   // Dynamically link some shared libs                    
172 //   if (gClassTable->GetID("AliRun") < 0) {
173 //     gROOT->LoadMacro("loadlibs.C");
174 //     loadlibs();
175 //   }
176
177
178   // Creating Run Loader and openning file containing Hits, Digits and RecPoints
179   AliRunLoader * RunLoader = AliRunLoader::Open(FileName,"Event","UPDATE");
180   if (RunLoader ==0x0) {
181     printf(">>> Error : Error Opening %s file \n",FileName);
182     return;
183   }
184   // Loading AliRun master
185   RunLoader->LoadgAlice();
186   gAlice = RunLoader->GetAliRun();
187
188   // Loading MUON subsystem
189   AliMUON * MUON = (AliMUON *) gAlice->GetDetector("MUON");
190   AliLoader * MUONLoader = RunLoader->GetLoader("MUONLoader");
191   MUONLoader->LoadHits("READ");
192   MUONLoader->LoadRecPoints("READ");
193
194   Int_t ievent, nevents;
195   nevents = RunLoader->GetNumberOfEvents();
196
197   // Initializations
198   // AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
199   AliMUONEventReconstructor *Reco = new AliMUONEventReconstructor();
200
201   Reco->SetRecGeantHits(RecGeantHits);
202
203   // The right place for changing AliMUONEventReconstructor parameters
204   // with respect to the default ones
205 //   Reco->SetMaxSigma2Distance(100.0);
206   Reco->SetPrintLevel(20);
207 //   Reco->SetPrintLevel(1);
208 //   Reco->SetBendingResolution(0.0);
209 //   Reco->SetNonBendingResolution(0.0);
210   cout << "AliMUONEventReconstructor: actual parameters" << endl;
211   Reco->Dump();
212 //   gObjectTable->Print();
213
214   // Loop over events
215   for (Int_t event = FirstEvent; event <= LastEvent; event++) {
216     MUON->ResetRawClusters();
217     cout << "Event: " << event << endl;
218     RunLoader->GetEvent(event);
219     // Addressing
220     MUON->SetTreeAddress(); 
221
222
223
224  //     Int_t nparticles = gAlice->GetEvent(event);
225 //      cout << "nparticles: " << nparticles << endl;
226     // prepare background file and/or event if necessary
227     if (RecGeantHits == 1) {
228       if (event == FirstEvent) Reco->SetBkgGeantFile(BkgGeantFileName);
229       if (Reco->GetBkgGeantFile())Reco->NextBkgGeantEvent();
230     }
231     Reco->EventReconstruct();
232     // Dump current event
233     Reco->EventDump();
234     // Fill Ntuple
235     AliMUONEventRecNtupleFill(Reco, 0);
236 //     gObjectTable->Print();
237
238   } // Event loop
239     // Write Ntuple
240     AliMUONEventRecNtupleFill(Reco, -1);
241 }