Line data Source code
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 20 : ClassImp( AliPHOSPIDv0)
92 :
93 : //____________________________________________________________________________
94 0 : AliPHOSPIDv0::AliPHOSPIDv0():
95 0 : fIDOptions("dis time"),
96 0 : fClusterizer(0),
97 0 : fTSMaker(0),
98 0 : fFormula(0),
99 0 : fDispersion(0.f),
100 0 : fCpvEmcDistance(0.f),
101 0 : fTimeGate(2.e-9f)
102 0 : {
103 : // default ctor
104 0 : }
105 :
106 : //____________________________________________________________________________
107 : AliPHOSPIDv0::AliPHOSPIDv0(AliPHOSGeometry *geom) :
108 0 : AliPHOSPID(geom),
109 0 : fIDOptions("dis time"),
110 0 : fClusterizer(0),
111 0 : fTSMaker(0),
112 0 : fFormula(new TFormula("LambdaCuts","(x>1)*(x<2.5)*(y>0)*(y<x)")),
113 0 : fDispersion(2.f),
114 0 : fCpvEmcDistance(3.f),
115 0 : fTimeGate(2.e-9f)
116 0 : {
117 : //ctor
118 0 : }
119 :
120 : //____________________________________________________________________________
121 : AliPHOSPIDv0::AliPHOSPIDv0(const AliPHOSPIDv0 & rhs) :
122 0 : AliPHOSPID(rhs),
123 0 : fIDOptions(rhs.fIDOptions),
124 0 : fClusterizer(rhs.fClusterizer),
125 0 : fTSMaker(rhs.fTSMaker),
126 0 : fFormula(rhs.fFormula),
127 0 : fDispersion(rhs.fDispersion),
128 0 : fCpvEmcDistance(rhs.fCpvEmcDistance),
129 0 : fTimeGate(rhs.fTimeGate)
130 0 : {
131 : //Copy ctor, the same as compiler-generated, possibly wrong if
132 : //someone implements dtor correctly.
133 0 : }
134 :
135 : //____________________________________________________________________________
136 : AliPHOSPIDv0 & AliPHOSPIDv0::operator = (const AliPHOSPIDv0 & rhs)
137 : {
138 : //Copy-assignment, emulates compiler generated, possibly wrong.
139 0 : AliPHOSPID::operator = (rhs);
140 0 : if (this != &rhs) {
141 0 : fIDOptions = rhs.fIDOptions;
142 0 : fClusterizer = rhs.fClusterizer;
143 0 : fTSMaker = rhs.fTSMaker;
144 0 : fFormula = rhs.fFormula;
145 0 : fDispersion = rhs.fDispersion;
146 0 : fCpvEmcDistance = rhs.fCpvEmcDistance;
147 0 : fTimeGate = rhs.fTimeGate;
148 0 : }
149 : else {
150 0 : AliFatal("Self assignment!");
151 : }
152 0 : return *this;
153 : }
154 :
155 : //____________________________________________________________________________
156 : AliPHOSPIDv0::~AliPHOSPIDv0()
157 0 : {
158 : //Empty dtor, fFormula leaks
159 0 : }
160 :
161 : //DP
162 : ////____________________________________________________________________________
163 : //Float_t AliPHOSPIDv0::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * Axis)const
164 : //{
165 : // // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
166 : //
167 : // const AliPHOSGeometry * geom = AliPHOSLoader::GetPHOSGeometry() ;
168 : // TVector3 vecEmc ;
169 : // TVector3 vecCpv ;
170 : //
171 : // emc->GetLocalPosition(vecEmc) ;
172 : // cpv->GetLocalPosition(vecCpv) ;
173 : // if(emc->GetPHOSMod() == cpv->GetPHOSMod()){
174 : //
175 : // // Correct to difference in CPV and EMC position due to different distance to center.
176 : // // we assume, that particle moves from center
177 : // Float_t dCPV = geom->GetIPtoOuterCoverDistance();
178 : // Float_t dEMC = geom->GetIPtoCrystalSurface() ;
179 : // dEMC = dEMC / dCPV ;
180 : // vecCpv = dEMC * vecCpv - vecEmc ;
181 : // if (Axis == "X") return vecCpv.X();
182 : // if (Axis == "Y") return vecCpv.Y();
183 : // if (Axis == "Z") return vecCpv.Z();
184 : // if (Axis == "R") return vecCpv.Mag();
185 : // }
186 : //
187 : // return 100000000 ;
188 : //}
189 :
190 : //____________________________________________________________________________
191 : void AliPHOSPIDv0::TrackSegments2RecParticles(Option_t * option)
192 : {
193 : //Steering method
194 :
195 0 : if(strstr(option,"tim"))
196 0 : gBenchmark->Start("PHOSPID");
197 :
198 0 : if(strstr(option,"print")) {
199 0 : Print() ;
200 0 : return ;
201 : }
202 :
203 0 : AliInfo(Form("%d emc clusters, %d track segments",
204 : fEMCRecPoints->GetEntriesFast(),
205 : fTrackSegments->GetEntriesFast())) ;
206 :
207 0 : MakeRecParticles() ;
208 :
209 0 : if(strstr(option,"deb"))
210 0 : PrintRecParticles(option) ;
211 :
212 0 : if(strstr(option,"tim")){
213 0 : gBenchmark->Stop("PHOSPID");
214 0 : AliInfo(Form("took %f seconds for PID",
215 : gBenchmark->GetCpuTime("PHOSPID")));
216 0 : }
217 :
218 0 : }
219 :
220 : //____________________________________________________________________________
221 : void AliPHOSPIDv0::MakeRecParticles()
222 : {
223 : // Reconstructs the particles from the tracksegments
224 :
225 0 : fRecParticles->Clear();
226 :
227 0 : TIter next(fTrackSegments) ;
228 : AliPHOSTrackSegment * ts ;
229 : Int_t index = 0 ;
230 : AliPHOSRecParticle * rp ;
231 :
232 0 : Bool_t ellips = fIDOptions.Contains("ell",TString::kIgnoreCase ) ;
233 0 : Bool_t disp = fIDOptions.Contains("dis",TString::kIgnoreCase ) ;
234 0 : Bool_t time = fIDOptions.Contains("tim",TString::kIgnoreCase ) ;
235 :
236 0 : while ( (ts = (AliPHOSTrackSegment *)next()) ) {
237 :
238 0 : new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
239 0 : rp = (AliPHOSRecParticle *)fRecParticles->At(index) ;
240 0 : rp->SetTrackSegment(index) ;
241 0 : rp->SetIndexInList(index) ;
242 :
243 : AliPHOSEmcRecPoint * emc = 0 ;
244 0 : if(ts->GetEmcIndex()>=0)
245 0 : emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
246 : else
247 0 : AliFatal(Form("ts->GetEmcIndex() return illegal index %d",ts->GetEmcIndex()));
248 :
249 : AliPHOSRecPoint * cpv = 0 ;
250 0 : if(ts->GetCpvIndex()>=0)
251 0 : cpv = (AliPHOSRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ;
252 :
253 : //set momentum and energy first
254 0 : Float_t e = emc->GetEnergy() ;
255 0 : TVector3 dir = GetMomentumDirection(emc,cpv) ;
256 0 : dir.SetMag(e) ;
257 :
258 0 : rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
259 0 : rp->SetCalcMass(0);
260 :
261 : //now set type (reconstructed) of the particle
262 : Int_t showerprofile = 0; // 0 narrow and 1 wide
263 :
264 0 : if(ellips){
265 0 : Float_t lambda[2] ;
266 0 : emc->GetElipsAxis(lambda) ;
267 0 : if(fFormula->Eval(lambda[0],lambda[1]) <= 0 )
268 0 : showerprofile = 1 ; // not narrow
269 0 : }
270 :
271 0 : if(disp)
272 0 : if(emc->GetDispersion() > fDispersion )
273 0 : showerprofile = 1 ; // not narrow
274 :
275 : Int_t slow = 0 ;
276 0 : if(time)
277 0 : if(emc->GetTime() > fTimeGate )
278 0 : slow = 0 ;
279 :
280 : // Looking at the CPV detector
281 : Int_t cpvdetector= 0 ; //1 hit and 0 no hit
282 0 : if(cpv)
283 0 : if(ts->GetCpvDistance("R") < fCpvEmcDistance)
284 0 : cpvdetector = 1 ;
285 :
286 0 : Int_t type = showerprofile + 2 * slow + 4 * cpvdetector ;
287 0 : rp->SetType(type) ;
288 0 : rp->SetProductionVertex(0,0,0,0);
289 0 : rp->SetFirstMother(-1);
290 0 : rp->SetLastMother(-1);
291 0 : rp->SetFirstDaughter(-1);
292 0 : rp->SetLastDaughter(-1);
293 0 : rp->SetPolarisation(0,0,0);
294 0 : index++ ;
295 0 : }
296 :
297 0 : }
298 :
299 : //____________________________________________________________________________
300 : void AliPHOSPIDv0:: Print(const Option_t *) const
301 : {
302 : // Print the parameters used for the particle type identification
303 0 : TString message ;
304 0 : message = "=============== AliPHOSPIDv0 ================\n" ;
305 0 : message += "Making PID\n" ;
306 0 : message += "with parameters:\n" ;
307 0 : message += " Maximal EMC - CPV distance (cm) %f\n" ;
308 0 : AliInfo(Form( message.Data(),
309 : fCpvEmcDistance ));
310 :
311 0 : if(fIDOptions.Contains("dis",TString::kIgnoreCase ))
312 0 : AliInfo(Form(" dispersion cut %f", fDispersion )) ;
313 0 : if(fIDOptions.Contains("ell",TString::kIgnoreCase ))
314 0 : AliInfo(Form(" Eliptic cuts function: %s",
315 : fFormula->GetTitle() )) ;
316 0 : if(fIDOptions.Contains("tim",TString::kIgnoreCase ))
317 0 : AliInfo(Form(" Time Gate used: %f", fTimeGate)) ;
318 0 : }
319 :
320 : //____________________________________________________________________________
321 : void AliPHOSPIDv0::SetShowerProfileCut(const char * formula)
322 : {
323 : //set shape of the cut on the axis of ellipce, drown around shouer
324 : //shower considered "narrow" if Formula(lambda[0],lambda[1]) > 0.
325 0 : if(fFormula)
326 0 : delete fFormula;
327 0 : fFormula = new TFormula("Lambda Cut",formula) ;
328 0 : }
329 :
330 : //____________________________________________________________________________
331 : void AliPHOSPIDv0::PlotDispersionCuts()const
332 : {
333 : // produces a plot of the dispersion cut
334 0 : TCanvas* lambdas = new TCanvas("lambdas","Cuts on the ellipse axis",200,10,700,500);
335 :
336 0 : if(fIDOptions.Contains("ell",TString::kIgnoreCase ) ){
337 0 : TF2 * ell = new TF2("Elliptic Cuts",fFormula->GetName(),0,3,0,3) ;
338 0 : ell->SetMinimum(0.0000001) ;
339 0 : ell->SetMaximum(0.001) ;
340 0 : ell->SetLineStyle(1) ;
341 0 : ell->SetLineWidth(2) ;
342 0 : ell->Draw() ;
343 0 : }
344 :
345 0 : if( fIDOptions.Contains("dis",TString::kIgnoreCase ) ){
346 0 : TF2 * dsp = new TF2("dispersion","(y<x)*(x*x+y*y < [0]*[0])",0,3,0,3) ;
347 0 : dsp->SetParameter(0,fDispersion) ;
348 0 : dsp->SetMinimum(0.0000001) ;
349 0 : dsp->SetMaximum(0.001) ;
350 0 : dsp->SetLineStyle(1) ;
351 0 : dsp->SetLineColor(2) ;
352 0 : dsp->SetLineWidth(2) ;
353 0 : dsp->SetNpx(200) ;
354 0 : dsp->SetNpy(200) ;
355 0 : if(fIDOptions.Contains("ell",TString::kIgnoreCase ) )
356 0 : dsp->Draw("same") ;
357 : else
358 0 : dsp->Draw() ;
359 0 : }
360 0 : lambdas->Update();
361 0 : }
362 :
363 : //____________________________________________________________________________
364 : TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * )const
365 : {
366 : // Calculates the momentum direction:
367 : // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint
368 : // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints
369 : // However because of the poor position resolution of PPSD the direction is always taken as if we were
370 : // in case 1.
371 :
372 0 : TVector3 dir(0,0,0) ;
373 :
374 0 : TVector3 emcglobalpos ;
375 0 : TMatrix dummy ;
376 :
377 0 : emc->GetGlobalPosition(emcglobalpos, dummy) ;
378 :
379 :
380 : // The following commented code becomes valid once the PPSD provides
381 : // a reasonable position resolution, at least as good as EMC !
382 : // TVector3 ppsdlglobalpos ;
383 : // TVector3 ppsduglobalpos ;
384 : // if( fPpsdLowRecPoint ){ // certainly a photon that has concerted
385 : // fPpsdLowRecPoint->GetGlobalPosition(ppsdlglobalpos, mdummy) ;
386 : // dir = emcglobalpos - ppsdlglobalpos ;
387 : // if( fPpsdUpRecPoint ){ // not looks like a charged
388 : // fPpsdUpRecPoint->GetGlobalPosition(ppsduglobalpos, mdummy) ;
389 : // dir = ( dir + emcglobalpos - ppsduglobalpos ) * 0.5 ;
390 : // }
391 : // }
392 : // else { // looks like a neutral
393 : // dir = emcglobalpos ;
394 : // }
395 :
396 0 : dir = emcglobalpos ;
397 0 : dir.SetZ( -dir.Z() ) ; // why ?
398 0 : dir.SetMag(1.) ;
399 :
400 : // One can not access MC information in the reconstruction!!
401 : // PLEASE FIT IT, EITHER BY TAKING 0,0,0 OR ACCESSING THE
402 : // VERTEX DIAMOND FROM CDB GRP FOLDER.
403 : //account correction to the position of IP
404 : // Float_t xo,yo,zo ; //Coordinates of the origin
405 : // gAlice->Generator()->GetOrigin(xo,yo,zo) ;
406 : // TVector3 origin(xo,yo,zo);
407 : // dir = dir - origin ;
408 :
409 : return dir ;
410 0 : }
411 : //____________________________________________________________________________
412 : void AliPHOSPIDv0::PrintRecParticles(Option_t * option)
413 : {
414 : // Print table of reconstructed particles
415 :
416 0 : TString message ;
417 0 : message = "Found %d RecParticles\n" ;
418 0 : AliInfo(Form(message.Data(),
419 : fRecParticles->GetEntriesFast() )) ;
420 :
421 0 : if(strstr(option,"all")) { // printing found TS
422 0 : AliInfo(" PARTICLE Index" ) ;
423 :
424 : Int_t index ;
425 0 : for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
426 0 : AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;
427 :
428 0 : Text_t particle[100];
429 0 : switch(rp->GetType()) {
430 : case AliPHOSFastRecParticle::kNEUTRALEMFAST:
431 0 : strncpy( particle, "NEUTRAL EM FAST",100);
432 : break;
433 : case AliPHOSFastRecParticle::kNEUTRALHAFAST:
434 0 : strncpy(particle, "NEUTRAL HA FAST",100);
435 : break;
436 : case AliPHOSFastRecParticle::kNEUTRALEMSLOW:
437 0 : strncpy(particle, "NEUTRAL EM SLOW",100);
438 : break ;
439 : case AliPHOSFastRecParticle::kNEUTRALHASLOW:
440 0 : strncpy(particle, "NEUTRAL HA SLOW",100);
441 : break ;
442 : case AliPHOSFastRecParticle::kCHARGEDEMFAST:
443 0 : strncpy(particle, "CHARGED EM FAST",100);
444 : break ;
445 : case AliPHOSFastRecParticle::kCHARGEDHAFAST:
446 0 : strncpy(particle, "CHARGED HA FAST",100);
447 : break ;
448 : case AliPHOSFastRecParticle::kCHARGEDEMSLOW:
449 0 : strncpy(particle, "CHARGEDEMSLOW",100);
450 : break ;
451 : case AliPHOSFastRecParticle::kCHARGEDHASLOW:
452 0 : strncpy(particle, "CHARGED HA SLOW",100);
453 : break ;
454 : }
455 :
456 : // Int_t * primaries;
457 : // Int_t nprimaries;
458 : // primaries = rp->GetPrimaries(nprimaries);
459 :
460 0 : AliInfo(Form(" %s %d",
461 : particle, rp->GetIndexInList())) ;
462 0 : }
463 0 : }
464 0 : }
|
