Modifications needed by the HBT analysis (P.Skowronski)
[u/mrichter/AliRoot.git] / 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 }