]>
Commit | Line | Data |
---|---|---|
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 | // Implementation version v1 of the EMCAL particle identifier | |
20 | ||
21 | //*-- Author: Yves Schutz (SUBATECH) | |
22 | // | |
23 | // --- ROOT system --- | |
24 | //#include "TROOT.h" | |
25 | #include "TTree.h" | |
26 | #include "TBenchmark.h" | |
27 | #include "TSystem.h" | |
28 | ||
29 | // --- Standard library --- | |
30 | ||
31 | // --- AliRoot header files --- | |
32 | #include "AliGenerator.h" | |
33 | #include "AliEMCALPIDv1.h" | |
34 | #include "AliEMCALRecParticle.h" | |
35 | #include "AliEMCALGetter.h" | |
36 | ||
37 | ClassImp( AliEMCALPIDv1) | |
38 | ||
39 | //____________________________________________________________________________ | |
40 | AliEMCALPIDv1::AliEMCALPIDv1():AliEMCALPID() | |
41 | { | |
42 | // default ctor | |
43 | ||
44 | InitParameters() ; | |
45 | fDefaultInit = kTRUE ; | |
46 | ||
47 | } | |
48 | ||
49 | //____________________________________________________________________________ | |
50 | AliEMCALPIDv1::AliEMCALPIDv1(const AliEMCALPIDv1 & pid ):AliEMCALPID(pid) | |
51 | { | |
52 | // ctor | |
53 | InitParameters() ; | |
54 | Init() ; | |
55 | ||
56 | } | |
57 | ||
58 | //____________________________________________________________________________ | |
59 | AliEMCALPIDv1::AliEMCALPIDv1(const TString alirunFileName, const TString eventFolderName):AliEMCALPID(alirunFileName, eventFolderName) | |
60 | { | |
61 | //ctor with the indication on where to look for the track segments | |
62 | ||
63 | InitParameters() ; | |
64 | Init() ; | |
65 | fDefaultInit = kFALSE ; | |
66 | ||
67 | } | |
68 | ||
69 | //____________________________________________________________________________ | |
70 | AliEMCALPIDv1::~AliEMCALPIDv1() | |
71 | { | |
72 | // dtor | |
73 | } | |
74 | ||
75 | //____________________________________________________________________________ | |
76 | const TString AliEMCALPIDv1::BranchName() const | |
77 | { | |
78 | ||
79 | return GetName() ; | |
80 | } | |
81 | ||
82 | //____________________________________________________________________________ | |
83 | Float_t AliEMCALPIDv1::GetCalibratedEnergy(Float_t e) const | |
84 | { | |
85 | // It calibrates Energy depending on the recpoint energy. | |
86 | // The energy of the reconstructed cluster is corrected with | |
87 | // the formula A + B* E + C* E^2, whose parameters where obtained | |
88 | // through the study of the reconstructed energy distribution of | |
89 | // monoenergetic photons. | |
90 | ||
91 | //Float_t p[]={0.,0.,0.}; | |
92 | //for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i); | |
93 | Float_t enerec = e ; // p[0] + p[1]*e + p[2]*e*e; | |
94 | return enerec ; | |
95 | ||
96 | } | |
97 | //____________________________________________________________________________ | |
98 | TVector3 AliEMCALPIDv1::GetMomentumDirection(AliEMCALRecPoint * emc)const | |
99 | { | |
100 | // Calculates the momentum direction: | |
101 | // direction is given by IP and this RecPoint | |
102 | ||
103 | ||
104 | TVector3 dir(0,0,0) ; | |
105 | TVector3 emcglobalpos ; | |
106 | // TMatrix dummy ; | |
107 | ||
108 | emc->GetGlobalPosition(emcglobalpos) ; | |
109 | ||
110 | ||
111 | dir = emcglobalpos ; | |
112 | // dir.SetMag(1.) ; Removed to avoid warings !!!!!!!!!!!!!! TO BE REVISED | |
113 | ||
114 | //account correction to the position of IP | |
115 | Float_t xo,yo,zo ; //Coordinates of the origin | |
116 | gAlice->Generator()->GetOrigin(xo,yo,zo) ; | |
117 | TVector3 origin(xo,yo,zo); | |
118 | dir = dir - origin ; | |
119 | ||
120 | return dir ; | |
121 | } | |
122 | ||
123 | //____________________________________________________________________________ | |
124 | void AliEMCALPIDv1::Init() | |
125 | { | |
126 | // Make all memory allocations that are not possible in default constructor | |
127 | // Add the PID task to the list of EMCAL tasks | |
128 | ||
129 | ||
130 | AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data()) ; | |
131 | ||
132 | if ( !gime->PID() ) | |
133 | gime->PostPID(this) ; | |
134 | } | |
135 | ||
136 | //____________________________________________________________________________ | |
137 | void AliEMCALPIDv1::InitParameters() | |
138 | { | |
139 | // Initialize the parameters | |
140 | fRecParticlesInRun = 0 ; | |
141 | fNEvent = 0 ; | |
142 | fRecParticlesInRun = 0 ; | |
143 | } | |
144 | ||
145 | //____________________________________________________________________________ | |
146 | ||
147 | void AliEMCALPIDv1::Exec(Option_t * option) | |
148 | { | |
149 | //Steering method | |
150 | ||
151 | if(strstr(option,"tim")) | |
152 | gBenchmark->Start("EMCALPID"); | |
153 | ||
154 | if(strstr(option,"print")) { | |
155 | Print("") ; | |
156 | return ; | |
157 | } | |
158 | AliEMCALGetter * gime = AliEMCALGetter::Instance() ; | |
159 | ||
160 | if (fLastEvent == -1) | |
161 | fLastEvent = gime->MaxEvent() - 1 ; | |
162 | else | |
163 | fLastEvent = TMath::Min(fLastEvent,gime->MaxEvent()); | |
164 | Int_t nEvents = fLastEvent - fFirstEvent + 1; | |
165 | ||
166 | Int_t ievent ; | |
167 | ||
168 | for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { | |
169 | gime->Event(ievent,"R") ; | |
170 | MakeRecParticles() ; | |
171 | WriteRecParticles(); | |
172 | if(strstr(option,"deb")) | |
173 | PrintRecParticles(option) ; | |
174 | //increment the total number of rec particles per run | |
175 | fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; | |
176 | } | |
177 | if(strstr(option,"tim")){ | |
178 | gBenchmark->Stop("EMCALPID"); | |
179 | printf("Exec: took %f seconds for PID %f seconds per event", | |
180 | gBenchmark->GetCpuTime("EMCALPID"), | |
181 | gBenchmark->GetCpuTime("EMCALPID")/nEvents) ; | |
182 | } | |
183 | ||
184 | Unload(); | |
185 | } | |
186 | ||
187 | //____________________________________________________________________________ | |
188 | void AliEMCALPIDv1::MakeRecParticles(){ | |
189 | ||
190 | // Makes a RecParticle out of a TrackSegment | |
191 | ||
192 | AliEMCALGetter * gime = AliEMCALGetter::Instance() ; | |
193 | TObjArray * aECARecPoints = gime->ECARecPoints() ; | |
194 | if ( !aECARecPoints ) { | |
195 | Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ; | |
196 | } | |
197 | TClonesArray * recParticles = gime->RecParticles() ; | |
198 | recParticles->Clear(); | |
199 | ||
200 | TIter next(aECARecPoints) ; | |
201 | AliEMCALRecPoint * eca ; | |
202 | Int_t index = 0 ; | |
203 | AliEMCALRecParticle * rp ; | |
204 | while ( (eca = (AliEMCALRecPoint *)next()) ) { | |
205 | ||
206 | new( (*recParticles)[index] ) AliEMCALRecParticle() ; | |
207 | rp = (AliEMCALRecParticle *)recParticles->At(index) ; | |
208 | rp->SetRecPoint(index) ; | |
209 | rp->SetIndexInList(index) ; | |
210 | ||
211 | // Now set type (reconstructed) of the particle | |
212 | ||
213 | // Choose the cluster energy range | |
214 | ||
215 | Float_t lambda[2] ; | |
216 | eca->GetElipsAxis(lambda) ; | |
217 | ||
218 | if((lambda[0]>0.01) && (lambda[1]>0.01)){ | |
219 | // Looking PCA. Define and calculate the data (X), | |
220 | // introduce in the function X2P that gives the components (P). | |
221 | ||
222 | Float_t spher = 0. ; | |
223 | Float_t emaxdtotal = 0. ; | |
224 | ||
225 | if((lambda[0]+lambda[1])!=0) | |
226 | spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); | |
227 | ||
228 | emaxdtotal=eca->GetMaximalEnergy()/eca->GetEnergy(); | |
229 | } | |
230 | ||
231 | // Float_t time = eca->GetTime() ; | |
232 | ||
233 | //Set momentum, energy and other parameters | |
234 | Float_t encal = GetCalibratedEnergy(eca->GetEnergy()); | |
235 | TVector3 dir = GetMomentumDirection(eca) ; | |
236 | // dir.SetMag(encal) ;Removed to avoid warings !!!!!!!!!!!!!! TO BE REVISED | |
237 | rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; | |
238 | rp->SetCalcMass(0); | |
239 | rp->Name(); //If photon sets the particle pdg name to gamma | |
240 | rp->SetProductionVertex(0,0,0,0); | |
241 | rp->SetFirstMother(-1); | |
242 | rp->SetLastMother(-1); | |
243 | rp->SetFirstDaughter(-1); | |
244 | rp->SetLastDaughter(-1); | |
245 | rp->SetPolarisation(0,0,0); | |
246 | //Set the position in global coordinate system from the RecPoint | |
247 | //AliEMCALGeometry * geom = gime->EMCALGeometry() ; | |
248 | //AliEMCALTowerRecPoint * erp = gime->ECARecPoint(rp->GetEMCALRPIndex()) ; | |
249 | TVector3 pos ; | |
250 | //geom->GetGlobal(erp, pos) ; !!!!!!!!!! to check | |
251 | rp->SetPos(pos); | |
252 | ||
253 | ||
254 | index++ ; | |
255 | } | |
256 | ||
257 | } | |
258 | ||
259 | //____________________________________________________________________________ | |
260 | void AliEMCALPIDv1:: Print(Option_t * /*option*/) const | |
261 | { | |
262 | // Print the parameters used for the particle type identification | |
263 | ||
264 | printf("Print: =============== AliEMCALPID1 ================") ; | |
265 | printf("Making PID\n"); | |
266 | printf(" Pricipal analysis file from 0.5 to 100 %s\n", fFileName.Data() ) ; | |
267 | printf(" Name of parameters file %s\n", fFileNamePar.Data() ) ; | |
268 | } | |
269 | ||
270 | //____________________________________________________________________________ | |
271 | void AliEMCALPIDv1::Print() const | |
272 | { | |
273 | // Print the parameters used for the particle type identification | |
274 | ||
275 | Info("Print", "=============== AliEMCALPIDv1 ================") ; | |
276 | } | |
277 | ||
278 | //____________________________________________________________________________ | |
279 | void AliEMCALPIDv1::PrintRecParticles(Option_t * option) | |
280 | { | |
281 | // Print table of reconstructed particles | |
282 | ||
283 | AliEMCALGetter *gime = AliEMCALGetter::Instance() ; | |
284 | ||
285 | TClonesArray * recParticles = gime->RecParticles() ; | |
286 | ||
287 | printf("\nevent %i", gAlice->GetEvNumber()); | |
288 | printf(" found %i", recParticles->GetEntriesFast()); | |
289 | printf(" RecParticles\n"); | |
290 | ||
291 | if(strstr(option,"all")) { // printing found TS | |
292 | printf("\n PARTICLE Index \n"); | |
293 | ||
294 | Int_t index ; | |
295 | for (index = 0 ; index < recParticles->GetEntries() ; index++) { | |
296 | AliEMCALRecParticle * rp = (AliEMCALRecParticle * ) recParticles->At(index) ; | |
297 | printf("\n"); | |
298 | printf(rp->Name().Data()); | |
299 | printf(" %i", rp->GetIndexInList()); | |
300 | printf(" %i", rp->GetType()); | |
301 | } | |
302 | } | |
303 | } | |
304 | ||
305 | //____________________________________________________________________________ | |
306 | void AliEMCALPIDv1::Unload() | |
307 | { | |
308 | // Unloads RecPoints, TrackSegments and RecParticles from the folder | |
309 | AliEMCALGetter * gime = AliEMCALGetter::Instance() ; | |
310 | gime->EmcalLoader()->UnloadRecPoints() ; | |
311 | gime->EmcalLoader()->UnloadTracks() ; | |
312 | gime->EmcalLoader()->UnloadRecParticles() ; | |
313 | } | |
314 | ||
315 | //____________________________________________________________________________ | |
316 | void AliEMCALPIDv1::WriteRecParticles() | |
317 | { | |
318 | // Write RecParticles array to file | |
319 | AliEMCALGetter *gime = AliEMCALGetter::Instance() ; | |
320 | ||
321 | TClonesArray * recParticles = gime->RecParticles() ; | |
322 | recParticles->Expand(recParticles->GetEntriesFast() ) ; | |
323 | TTree * treeP = gime->TreeP() ; | |
324 | ||
325 | ||
326 | ||
327 | //First rp | |
328 | Int_t bufferSize = 32000 ; | |
329 | TBranch * rpBranch = treeP->Branch("EMCALRP",&recParticles,bufferSize); | |
330 | rpBranch->SetTitle(BranchName()); | |
331 | ||
332 | rpBranch->Fill() ; | |
333 | ||
334 | gime->WriteRecParticles("OVERWRITE"); | |
335 | gime->WritePID("OVERWRITE"); | |
336 | ||
337 | } | |
338 |