]>
Commit | Line | Data |
---|---|---|
c2a319d4 | 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 | ||
5ac84131 | 30 | #include "AliHeader.h" |
c2a319d4 | 31 | #include "AliRun.h" |
32 | ||
88cb7938 | 33 | #include "AliMUON.h" |
c2a319d4 | 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) { | |
88cb7938 | 96 | printf(">>> Writing Ntuple of reconstructed tracks\n"); |
c2a319d4 | 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... | |
88cb7938 | 107 | printf(">>> Creating Ntuple of reconstructed tracks\n"); |
c2a319d4 | 108 | ntuple = new TTree("MUONtrackReco", "MUONtrackReco"); |
109 | ntuple->Branch("Header","AliMUONHeaderRecNtuple", &header); | |
110 | ntuple->Branch("Tracks", &recTracks); | |
111 | } | |
112 | ||
113 | // header | |
88cb7938 | 114 | |
c2a319d4 | 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(); | |
04c865aa | 144 | printf("test> Px %f Py %f Pz %f \n", recTrackNt->fPxRec, recTrackNt->fPyRec, recTrackNt->fPzRec); |
c2a319d4 | 145 | // track parameters at vertex of best compatible generated track: |
146 | // in fact muon with the right charge | |
04c865aa | 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 | } | |
c2a319d4 | 162 | } // for (trackIndex = 0;... |
163 | ||
88cb7938 | 164 | printf(">>> Filling Ntuple of reconstructed tracks\n"); |
c2a319d4 | 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 | ||
88cb7938 | 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; | |
c2a319d4 | 191 | } |
88cb7938 | 192 | // Loading AliRun master |
193 | RunLoader->LoadgAlice(); | |
194 | gAlice = RunLoader->GetAliRun(); | |
04c865aa | 195 | RunLoader->LoadKinematics("READ"); |
196 | ||
88cb7938 | 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(); | |
c2a319d4 | 205 | |
206 | // Initializations | |
207 | // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? | |
896c4449 | 208 | MUON->SetTreeAddress(); |
c2a319d4 | 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); | |
04c865aa | 216 | // Reco->SetPrintLevel(20); |
217 | Reco->SetPrintLevel(1); | |
c2a319d4 | 218 | // Reco->SetBendingResolution(0.0); |
219 | // Reco->SetNonBendingResolution(0.0); | |
220 | cout << "AliMUONEventReconstructor: actual parameters" << endl; | |
221 | Reco->Dump(); | |
222 | // gObjectTable->Print(); | |
c2a319d4 | 223 | // Loop over events |
896c4449 | 224 | for (Int_t event = FirstEvent; event < LastEvent; event++) { |
c2a319d4 | 225 | cout << "Event: " << event << endl; |
896c4449 | 226 | RunLoader->GetEvent(event); |
227 | // Int_t nparticles = gAlice->GetEvent(event); | |
228 | // cout << "nparticles: " << nparticles << endl; | |
c2a319d4 | 229 | // prepare background file and/or event if necessary |
230 | if (RecGeantHits == 1) { | |
231 | if (event == FirstEvent) Reco->SetBkgGeantFile(BkgGeantFileName); | |
232 | if (Reco->GetBkgGeantFile())Reco->NextBkgGeantEvent(); | |
233 | } | |
234 | Reco->EventReconstruct(); | |
235 | // Dump current event | |
236 | Reco->EventDump(); | |
237 | // Fill Ntuple | |
238 | AliMUONEventRecNtupleFill(Reco, 0); | |
239 | // gObjectTable->Print(); | |
896c4449 | 240 | MUON->ResetRawClusters(); |
c2a319d4 | 241 | } // Event loop |
242 | // Write Ntuple | |
243 | AliMUONEventRecNtupleFill(Reco, -1); | |
244 | } |