]>
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 | /* History of cvs commits: | |
19 | * | |
20 | * $Log$ | |
21 | * Revision 1.15 2007/03/06 06:57:05 kharlov | |
22 | * DP:calculation of distance to CPV done in TSM | |
23 | * | |
24 | * Revision 1.14 2006/09/07 18:31:08 kharlov | |
25 | * Effective c++ corrections (T.Pocheptsov) | |
26 | * | |
27 | * Revision 1.13 2005/05/28 14:19:04 schutz | |
28 | * Compilation warnings fixed by T.P. | |
29 | * | |
30 | */ | |
31 | ||
32 | //_________________________________________________________________________ | |
33 | // Implementation version v0 of the PHOS particle identifier | |
34 | // Particle identification based on the | |
35 | // - CPV information, | |
36 | // - Preshower information (in MIXT or GPS2 geometries) | |
37 | // - shower width. | |
38 | // | |
39 | // CPV or Preshower clusters should be closer in PHOS plane than fCpvEmcDistance (in cm). | |
40 | // This parameter can be set by method SetCpvtoEmcDistanceCut(Float_t cut) | |
41 | // | |
42 | // One can set desirable ID method by the function SetIdentificationMethod(option). | |
43 | // Presently the following options can be used together or separately : | |
44 | // - "disp": use dispersion cut on shower width | |
45 | // (width can be set by method SetDispersionCut(Float_t cut) | |
46 | // - "ell" : use cut on the axis of the ellipse, drawn around shower | |
47 | // (this cut can be changed by SetShowerProfileCut(char* formula), | |
48 | // where formula - any function of two variables f(lambda[0],lambda[1]). | |
49 | // Shower is considered as EM if f() > 0 ) | |
50 | // One can visualize current cuts calling method PlotDispersionCuts(). | |
51 | // | |
52 | // use case: | |
53 | // root [0] AliPHOSPIDv0 * p1 = new AliPHOSPIDv0("galice.root") | |
54 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated | |
55 | // root [1] p1->SetIdentificationMethod("disp ellipse") | |
56 | // root [2] p1->ExecuteTask() | |
57 | // root [3] AliPHOSPIDv0 * p2 = new AliPHOSPIDv0("galice1.root","ts1") | |
58 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated | |
59 | // // reading headers from file galice1.root and TrackSegments | |
60 | // // with title "ts1" | |
61 | // root [4] p2->SetRecParticlesBranch("rp1") | |
62 | // // set file name for the branch RecParticles | |
63 | // root [5] p2->ExecuteTask("deb all time") | |
64 | // // available options | |
65 | // // "deb" - prints # of reconstructed particles | |
66 | // // "deb all" - prints # and list of RecParticles | |
67 | // // "time" - prints benchmarking results | |
68 | // | |
69 | //*-- Author: Yves Schutz (SUBATECH) & Gines Martinez (SUBATECH) & | |
70 | // Dmitri Peressounko (SUBATECH & Kurchatov Institute) | |
71 | // Completely redesined by Dmitri Peressounko, March 2001 | |
72 | ||
73 | // --- ROOT system --- | |
74 | #include "TTree.h" | |
75 | #include "TF2.h" | |
76 | #include "TFormula.h" | |
77 | #include "TCanvas.h" | |
78 | #include "TClonesArray.h" | |
79 | ||
80 | #include "TBenchmark.h" | |
81 | // --- Standard library --- | |
82 | ||
83 | // --- AliRoot header files --- | |
84 | #include "AliLog.h" | |
85 | #include "AliPHOSPIDv0.h" | |
86 | #include "AliPHOSEmcRecPoint.h" | |
87 | #include "AliPHOSTrackSegment.h" | |
88 | #include "AliPHOSRecParticle.h" | |
89 | #include "AliPHOSGeometry.h" | |
90 | ||
91 | ClassImp( AliPHOSPIDv0) | |
92 | ||
93 | //____________________________________________________________________________ | |
94 | AliPHOSPIDv0::AliPHOSPIDv0(): | |
95 | fIDOptions("dis time"), | |
96 | fClusterizer(0), | |
97 | fTSMaker(0), | |
98 | fFormula(0), | |
99 | fDispersion(0.f), | |
100 | fCpvEmcDistance(0.f), | |
101 | fTimeGate(2.e-9f) | |
102 | { | |
103 | // default ctor | |
104 | } | |
105 | ||
106 | //____________________________________________________________________________ | |
107 | AliPHOSPIDv0::AliPHOSPIDv0(AliPHOSGeometry *geom) : | |
108 | AliPHOSPID(geom), | |
109 | fIDOptions("dis time"), | |
110 | fClusterizer(0), | |
111 | fTSMaker(0), | |
112 | fFormula(new TFormula("LambdaCuts","(x>1)*(x<2.5)*(y>0)*(y<x)")), | |
113 | fDispersion(2.f), | |
114 | fCpvEmcDistance(3.f), | |
115 | fTimeGate(2.e-9f) | |
116 | { | |
117 | //ctor | |
118 | } | |
119 | ||
120 | //____________________________________________________________________________ | |
121 | AliPHOSPIDv0::AliPHOSPIDv0(const AliPHOSPIDv0 & rhs) : | |
122 | AliPHOSPID(rhs), | |
123 | fIDOptions(rhs.fIDOptions), | |
124 | fClusterizer(rhs.fClusterizer), | |
125 | fTSMaker(rhs.fTSMaker), | |
126 | fFormula(rhs.fFormula), | |
127 | fDispersion(rhs.fDispersion), | |
128 | fCpvEmcDistance(rhs.fCpvEmcDistance), | |
129 | fTimeGate(rhs.fTimeGate) | |
130 | { | |
131 | //Copy ctor, the same as compiler-generated, possibly wrong if | |
132 | //someone implements dtor correctly. | |
133 | } | |
134 | ||
135 | //____________________________________________________________________________ | |
136 | AliPHOSPIDv0 & AliPHOSPIDv0::operator = (const AliPHOSPIDv0 & rhs) | |
137 | { | |
138 | //Copy-assignment, emulates compiler generated, possibly wrong. | |
139 | AliPHOSPID::operator = (rhs); | |
140 | fIDOptions = rhs.fIDOptions; | |
141 | fClusterizer = rhs.fClusterizer; | |
142 | fTSMaker = rhs.fTSMaker; | |
143 | fFormula = rhs.fFormula; | |
144 | fDispersion = rhs.fDispersion; | |
145 | fCpvEmcDistance = rhs.fCpvEmcDistance; | |
146 | fTimeGate = rhs.fTimeGate; | |
147 | ||
148 | return *this; | |
149 | } | |
150 | ||
151 | //____________________________________________________________________________ | |
152 | AliPHOSPIDv0::~AliPHOSPIDv0() | |
153 | { | |
154 | //Empty dtor, fFormula leaks | |
155 | } | |
156 | ||
157 | //DP | |
158 | ////____________________________________________________________________________ | |
159 | //Float_t AliPHOSPIDv0::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * Axis)const | |
160 | //{ | |
161 | // // Calculates the distance between the EMC RecPoint and the PPSD RecPoint | |
162 | // | |
163 | // const AliPHOSGeometry * geom = AliPHOSLoader::GetPHOSGeometry() ; | |
164 | // TVector3 vecEmc ; | |
165 | // TVector3 vecCpv ; | |
166 | // | |
167 | // emc->GetLocalPosition(vecEmc) ; | |
168 | // cpv->GetLocalPosition(vecCpv) ; | |
169 | // if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ | |
170 | // | |
171 | // // Correct to difference in CPV and EMC position due to different distance to center. | |
172 | // // we assume, that particle moves from center | |
173 | // Float_t dCPV = geom->GetIPtoOuterCoverDistance(); | |
174 | // Float_t dEMC = geom->GetIPtoCrystalSurface() ; | |
175 | // dEMC = dEMC / dCPV ; | |
176 | // vecCpv = dEMC * vecCpv - vecEmc ; | |
177 | // if (Axis == "X") return vecCpv.X(); | |
178 | // if (Axis == "Y") return vecCpv.Y(); | |
179 | // if (Axis == "Z") return vecCpv.Z(); | |
180 | // if (Axis == "R") return vecCpv.Mag(); | |
181 | // } | |
182 | // | |
183 | // return 100000000 ; | |
184 | //} | |
185 | ||
186 | //____________________________________________________________________________ | |
187 | void AliPHOSPIDv0::TrackSegments2RecParticles(Option_t * option) | |
188 | { | |
189 | //Steering method | |
190 | ||
191 | if(strstr(option,"tim")) | |
192 | gBenchmark->Start("PHOSPID"); | |
193 | ||
194 | if(strstr(option,"print")) { | |
195 | Print() ; | |
196 | return ; | |
197 | } | |
198 | ||
199 | AliInfo(Form("%d emc clusters, %d track segments", | |
200 | fEMCRecPoints->GetEntriesFast(), | |
201 | fTrackSegments->GetEntriesFast())) ; | |
202 | ||
203 | MakeRecParticles() ; | |
204 | ||
205 | if(strstr(option,"deb")) | |
206 | PrintRecParticles(option) ; | |
207 | ||
208 | if(strstr(option,"tim")){ | |
209 | gBenchmark->Stop("PHOSPID"); | |
210 | AliInfo(Form("took %f seconds for PID", | |
211 | gBenchmark->GetCpuTime("PHOSPID"))); | |
212 | } | |
213 | ||
214 | } | |
215 | ||
216 | //____________________________________________________________________________ | |
217 | void AliPHOSPIDv0::MakeRecParticles() | |
218 | { | |
219 | // Reconstructs the particles from the tracksegments | |
220 | ||
221 | fRecParticles->Clear(); | |
222 | ||
223 | TIter next(fTrackSegments) ; | |
224 | AliPHOSTrackSegment * ts ; | |
225 | Int_t index = 0 ; | |
226 | AliPHOSRecParticle * rp ; | |
227 | ||
228 | Bool_t ellips = fIDOptions.Contains("ell",TString::kIgnoreCase ) ; | |
229 | Bool_t disp = fIDOptions.Contains("dis",TString::kIgnoreCase ) ; | |
230 | Bool_t time = fIDOptions.Contains("tim",TString::kIgnoreCase ) ; | |
231 | ||
232 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { | |
233 | ||
234 | new( (*fRecParticles)[index] ) AliPHOSRecParticle() ; | |
235 | rp = (AliPHOSRecParticle *)fRecParticles->At(index) ; | |
236 | rp->SetTrackSegment(index) ; | |
237 | rp->SetIndexInList(index) ; | |
238 | ||
239 | AliPHOSEmcRecPoint * emc = 0 ; | |
240 | if(ts->GetEmcIndex()>=0) | |
241 | emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ; | |
242 | ||
243 | AliPHOSRecPoint * cpv = 0 ; | |
244 | if(ts->GetCpvIndex()>=0) | |
245 | cpv = (AliPHOSRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ; | |
246 | ||
247 | //set momentum and energy first | |
248 | Float_t e = emc->GetEnergy() ; | |
249 | TVector3 dir = GetMomentumDirection(emc,cpv) ; | |
250 | dir.SetMag(e) ; | |
251 | ||
252 | rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ; | |
253 | rp->SetCalcMass(0); | |
254 | ||
255 | //now set type (reconstructed) of the particle | |
256 | Int_t showerprofile = 0; // 0 narrow and 1 wide | |
257 | ||
258 | if(ellips){ | |
259 | Float_t lambda[2] ; | |
260 | emc->GetElipsAxis(lambda) ; | |
261 | if(fFormula->Eval(lambda[0],lambda[1]) <= 0 ) | |
262 | showerprofile = 1 ; // not narrow | |
263 | } | |
264 | ||
265 | if(disp) | |
266 | if(emc->GetDispersion() > fDispersion ) | |
267 | showerprofile = 1 ; // not narrow | |
268 | ||
269 | Int_t slow = 0 ; | |
270 | if(time) | |
271 | if(emc->GetTime() > fTimeGate ) | |
272 | slow = 0 ; | |
273 | ||
274 | // Looking at the CPV detector | |
275 | Int_t cpvdetector= 0 ; //1 hit and 0 no hit | |
276 | if(cpv) | |
277 | if(ts->GetCpvDistance("R") < fCpvEmcDistance) | |
278 | cpvdetector = 1 ; | |
279 | ||
280 | Int_t type = showerprofile + 2 * slow + 4 * cpvdetector ; | |
281 | rp->SetType(type) ; | |
282 | rp->SetProductionVertex(0,0,0,0); | |
283 | rp->SetFirstMother(-1); | |
284 | rp->SetLastMother(-1); | |
285 | rp->SetFirstDaughter(-1); | |
286 | rp->SetLastDaughter(-1); | |
287 | rp->SetPolarisation(0,0,0); | |
288 | index++ ; | |
289 | } | |
290 | ||
291 | } | |
292 | ||
293 | //____________________________________________________________________________ | |
294 | void AliPHOSPIDv0:: Print(const Option_t *) const | |
295 | { | |
296 | // Print the parameters used for the particle type identification | |
297 | TString message ; | |
298 | message = "=============== AliPHOSPIDv0 ================\n" ; | |
299 | message += "Making PID\n" ; | |
300 | message += "with parameters:\n" ; | |
301 | message += " Maximal EMC - CPV distance (cm) %f\n" ; | |
302 | AliInfo(Form( message.Data(), | |
303 | fCpvEmcDistance )); | |
304 | ||
305 | if(fIDOptions.Contains("dis",TString::kIgnoreCase )) | |
306 | AliInfo(Form(" dispersion cut %f", fDispersion )) ; | |
307 | if(fIDOptions.Contains("ell",TString::kIgnoreCase )) | |
308 | AliInfo(Form(" Eliptic cuts function: %s", | |
309 | fFormula->GetTitle() )) ; | |
310 | if(fIDOptions.Contains("tim",TString::kIgnoreCase )) | |
311 | AliInfo(Form(" Time Gate used: %f", fTimeGate)) ; | |
312 | } | |
313 | ||
314 | //____________________________________________________________________________ | |
315 | void AliPHOSPIDv0::SetShowerProfileCut(const char * formula) | |
316 | { | |
317 | //set shape of the cut on the axis of ellipce, drown around shouer | |
318 | //shower considered "narrow" if Formula(lambda[0],lambda[1]) > 0. | |
319 | if(fFormula) | |
320 | delete fFormula; | |
321 | fFormula = new TFormula("Lambda Cut",formula) ; | |
322 | } | |
323 | ||
324 | //____________________________________________________________________________ | |
325 | void AliPHOSPIDv0::PlotDispersionCuts()const | |
326 | { | |
327 | // produces a plot of the dispersion cut | |
328 | TCanvas* lambdas = new TCanvas("lambdas","Cuts on the ellipse axis",200,10,700,500); | |
329 | ||
330 | if(fIDOptions.Contains("ell",TString::kIgnoreCase ) ){ | |
331 | TF2 * ell = new TF2("Elliptic Cuts",fFormula->GetName(),0,3,0,3) ; | |
332 | ell->SetMinimum(0.0000001) ; | |
333 | ell->SetMaximum(0.001) ; | |
334 | ell->SetLineStyle(1) ; | |
335 | ell->SetLineWidth(2) ; | |
336 | ell->Draw() ; | |
337 | } | |
338 | ||
339 | if( fIDOptions.Contains("dis",TString::kIgnoreCase ) ){ | |
340 | TF2 * dsp = new TF2("dispersion","(y<x)*(x*x+y*y < [0]*[0])",0,3,0,3) ; | |
341 | dsp->SetParameter(0,fDispersion) ; | |
342 | dsp->SetMinimum(0.0000001) ; | |
343 | dsp->SetMaximum(0.001) ; | |
344 | dsp->SetLineStyle(1) ; | |
345 | dsp->SetLineColor(2) ; | |
346 | dsp->SetLineWidth(2) ; | |
347 | dsp->SetNpx(200) ; | |
348 | dsp->SetNpy(200) ; | |
349 | if(fIDOptions.Contains("ell",TString::kIgnoreCase ) ) | |
350 | dsp->Draw("same") ; | |
351 | else | |
352 | dsp->Draw() ; | |
353 | } | |
354 | lambdas->Update(); | |
355 | } | |
356 | ||
357 | //____________________________________________________________________________ | |
358 | TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * )const | |
359 | { | |
360 | // Calculates the momentum direction: | |
361 | // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint | |
362 | // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints | |
363 | // However because of the poor position resolution of PPSD the direction is always taken as if we were | |
364 | // in case 1. | |
365 | ||
366 | TVector3 dir(0,0,0) ; | |
367 | ||
368 | TVector3 emcglobalpos ; | |
369 | TMatrix dummy ; | |
370 | ||
371 | emc->GetGlobalPosition(emcglobalpos, dummy) ; | |
372 | ||
373 | ||
374 | // The following commented code becomes valid once the PPSD provides | |
375 | // a reasonable position resolution, at least as good as EMC ! | |
376 | // TVector3 ppsdlglobalpos ; | |
377 | // TVector3 ppsduglobalpos ; | |
378 | // if( fPpsdLowRecPoint ){ // certainly a photon that has concerted | |
379 | // fPpsdLowRecPoint->GetGlobalPosition(ppsdlglobalpos, mdummy) ; | |
380 | // dir = emcglobalpos - ppsdlglobalpos ; | |
381 | // if( fPpsdUpRecPoint ){ // not looks like a charged | |
382 | // fPpsdUpRecPoint->GetGlobalPosition(ppsduglobalpos, mdummy) ; | |
383 | // dir = ( dir + emcglobalpos - ppsduglobalpos ) * 0.5 ; | |
384 | // } | |
385 | // } | |
386 | // else { // looks like a neutral | |
387 | // dir = emcglobalpos ; | |
388 | // } | |
389 | ||
390 | dir = emcglobalpos ; | |
391 | dir.SetZ( -dir.Z() ) ; // why ? | |
392 | dir.SetMag(1.) ; | |
393 | ||
394 | // One can not access MC information in the reconstruction!! | |
395 | // PLEASE FIT IT, EITHER BY TAKING 0,0,0 OR ACCESSING THE | |
396 | // VERTEX DIAMOND FROM CDB GRP FOLDER. | |
397 | //account correction to the position of IP | |
398 | // Float_t xo,yo,zo ; //Coordinates of the origin | |
399 | // gAlice->Generator()->GetOrigin(xo,yo,zo) ; | |
400 | // TVector3 origin(xo,yo,zo); | |
401 | // dir = dir - origin ; | |
402 | ||
403 | return dir ; | |
404 | } | |
405 | //____________________________________________________________________________ | |
406 | void AliPHOSPIDv0::PrintRecParticles(Option_t * option) | |
407 | { | |
408 | // Print table of reconstructed particles | |
409 | ||
410 | TString message ; | |
411 | message = "Found %d RecParticles\n" ; | |
412 | AliInfo(Form(message.Data(), | |
413 | fRecParticles->GetEntriesFast() )) ; | |
414 | ||
415 | if(strstr(option,"all")) { // printing found TS | |
416 | AliInfo(" PARTICLE Index" ) ; | |
417 | ||
418 | Int_t index ; | |
419 | for (index = 0 ; index < fRecParticles->GetEntries() ; index++) { | |
420 | AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ; | |
421 | ||
422 | Text_t particle[11]; | |
423 | switch(rp->GetType()) { | |
424 | case AliPHOSFastRecParticle::kNEUTRALEMFAST: | |
425 | strcpy( particle, "NEUTRAL EM FAST"); | |
426 | break; | |
427 | case AliPHOSFastRecParticle::kNEUTRALHAFAST: | |
428 | strcpy(particle, "NEUTRAL HA FAST"); | |
429 | break; | |
430 | case AliPHOSFastRecParticle::kNEUTRALEMSLOW: | |
431 | strcpy(particle, "NEUTRAL EM SLOW"); | |
432 | break ; | |
433 | case AliPHOSFastRecParticle::kNEUTRALHASLOW: | |
434 | strcpy(particle, "NEUTRAL HA SLOW"); | |
435 | break ; | |
436 | case AliPHOSFastRecParticle::kCHARGEDEMFAST: | |
437 | strcpy(particle, "CHARGED EM FAST") ; | |
438 | break ; | |
439 | case AliPHOSFastRecParticle::kCHARGEDHAFAST: | |
440 | strcpy(particle, "CHARGED HA FAST") ; | |
441 | break ; | |
442 | case AliPHOSFastRecParticle::kCHARGEDEMSLOW: | |
443 | strcpy(particle, "CHARGEDEMSLOW") ; | |
444 | break ; | |
445 | case AliPHOSFastRecParticle::kCHARGEDHASLOW: | |
446 | strcpy(particle, "CHARGED HA SLOW") ; | |
447 | break ; | |
448 | } | |
449 | ||
450 | // Int_t * primaries; | |
451 | // Int_t nprimaries; | |
452 | // primaries = rp->GetPrimaries(nprimaries); | |
453 | ||
454 | AliInfo(Form(" %s %d", | |
455 | particle, rp->GetIndexInList())) ; | |
456 | } | |
457 | } | |
458 | } |