1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /* History of cvs commits:
21 * Revision 1.15 2007/03/06 06:57:05 kharlov
22 * DP:calculation of distance to CPV done in TSM
24 * Revision 1.14 2006/09/07 18:31:08 kharlov
25 * Effective c++ corrections (T.Pocheptsov)
27 * Revision 1.13 2005/05/28 14:19:04 schutz
28 * Compilation warnings fixed by T.P.
32 //_________________________________________________________________________
33 // Implementation version v0 of the PHOS particle identifier
34 // Particle identification based on the
36 // - Preshower information (in MIXT or GPS2 geometries)
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)
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().
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
69 //*-- Author: Yves Schutz (SUBATECH) & Gines Martinez (SUBATECH) &
70 // Dmitri Peressounko (SUBATECH & Kurchatov Institute)
71 // Completely redesined by Dmitri Peressounko, March 2001
73 // --- ROOT system ---
78 #include "TClonesArray.h"
80 #include "TBenchmark.h"
81 // --- Standard library ---
83 // --- AliRoot header files ---
85 #include "AliPHOSPIDv0.h"
86 #include "AliPHOSEmcRecPoint.h"
87 #include "AliPHOSTrackSegment.h"
88 #include "AliPHOSRecParticle.h"
89 #include "AliPHOSGeometry.h"
91 ClassImp( AliPHOSPIDv0)
93 //____________________________________________________________________________
94 AliPHOSPIDv0::AliPHOSPIDv0():
95 fIDOptions("dis time"),
100 fCpvEmcDistance(0.f),
106 //____________________________________________________________________________
107 AliPHOSPIDv0::AliPHOSPIDv0(AliPHOSGeometry *geom) :
109 fIDOptions("dis time"),
112 fFormula(new TFormula("LambdaCuts","(x>1)*(x<2.5)*(y>0)*(y<x)")),
114 fCpvEmcDistance(3.f),
120 //____________________________________________________________________________
121 AliPHOSPIDv0::AliPHOSPIDv0(const AliPHOSPIDv0 & 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)
131 //Copy ctor, the same as compiler-generated, possibly wrong if
132 //someone implements dtor correctly.
135 //____________________________________________________________________________
136 AliPHOSPIDv0 & AliPHOSPIDv0::operator = (const AliPHOSPIDv0 & rhs)
138 //Copy-assignment, emulates compiler generated, possibly wrong.
139 AliPHOSPID::operator = (rhs);
141 fIDOptions = rhs.fIDOptions;
142 fClusterizer = rhs.fClusterizer;
143 fTSMaker = rhs.fTSMaker;
144 fFormula = rhs.fFormula;
145 fDispersion = rhs.fDispersion;
146 fCpvEmcDistance = rhs.fCpvEmcDistance;
147 fTimeGate = rhs.fTimeGate;
150 AliFatal("Self assignment!");
155 //____________________________________________________________________________
156 AliPHOSPIDv0::~AliPHOSPIDv0()
158 //Empty dtor, fFormula leaks
162 ////____________________________________________________________________________
163 //Float_t AliPHOSPIDv0::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * Axis)const
165 // // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
167 // const AliPHOSGeometry * geom = AliPHOSLoader::GetPHOSGeometry() ;
171 // emc->GetLocalPosition(vecEmc) ;
172 // cpv->GetLocalPosition(vecCpv) ;
173 // if(emc->GetPHOSMod() == cpv->GetPHOSMod()){
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();
187 // return 100000000 ;
190 //____________________________________________________________________________
191 void AliPHOSPIDv0::TrackSegments2RecParticles(Option_t * option)
195 if(strstr(option,"tim"))
196 gBenchmark->Start("PHOSPID");
198 if(strstr(option,"print")) {
203 AliInfo(Form("%d emc clusters, %d track segments",
204 fEMCRecPoints->GetEntriesFast(),
205 fTrackSegments->GetEntriesFast())) ;
209 if(strstr(option,"deb"))
210 PrintRecParticles(option) ;
212 if(strstr(option,"tim")){
213 gBenchmark->Stop("PHOSPID");
214 AliInfo(Form("took %f seconds for PID",
215 gBenchmark->GetCpuTime("PHOSPID")));
220 //____________________________________________________________________________
221 void AliPHOSPIDv0::MakeRecParticles()
223 // Reconstructs the particles from the tracksegments
225 fRecParticles->Clear();
227 TIter next(fTrackSegments) ;
228 AliPHOSTrackSegment * ts ;
230 AliPHOSRecParticle * rp ;
232 Bool_t ellips = fIDOptions.Contains("ell",TString::kIgnoreCase ) ;
233 Bool_t disp = fIDOptions.Contains("dis",TString::kIgnoreCase ) ;
234 Bool_t time = fIDOptions.Contains("tim",TString::kIgnoreCase ) ;
236 while ( (ts = (AliPHOSTrackSegment *)next()) ) {
238 new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
239 rp = (AliPHOSRecParticle *)fRecParticles->At(index) ;
240 rp->SetTrackSegment(index) ;
241 rp->SetIndexInList(index) ;
243 AliPHOSEmcRecPoint * emc = 0 ;
244 if(ts->GetEmcIndex()>=0)
245 emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
247 AliFatal(Form("ts->GetEmcIndex() return illegal index %d",ts->GetEmcIndex()));
249 AliPHOSRecPoint * cpv = 0 ;
250 if(ts->GetCpvIndex()>=0)
251 cpv = (AliPHOSRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ;
253 //set momentum and energy first
254 Float_t e = emc->GetEnergy() ;
255 TVector3 dir = GetMomentumDirection(emc,cpv) ;
258 rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
261 //now set type (reconstructed) of the particle
262 Int_t showerprofile = 0; // 0 narrow and 1 wide
266 emc->GetElipsAxis(lambda) ;
267 if(fFormula->Eval(lambda[0],lambda[1]) <= 0 )
268 showerprofile = 1 ; // not narrow
272 if(emc->GetDispersion() > fDispersion )
273 showerprofile = 1 ; // not narrow
277 if(emc->GetTime() > fTimeGate )
280 // Looking at the CPV detector
281 Int_t cpvdetector= 0 ; //1 hit and 0 no hit
283 if(ts->GetCpvDistance("R") < fCpvEmcDistance)
286 Int_t type = showerprofile + 2 * slow + 4 * cpvdetector ;
288 rp->SetProductionVertex(0,0,0,0);
289 rp->SetFirstMother(-1);
290 rp->SetLastMother(-1);
291 rp->SetFirstDaughter(-1);
292 rp->SetLastDaughter(-1);
293 rp->SetPolarisation(0,0,0);
299 //____________________________________________________________________________
300 void AliPHOSPIDv0:: Print(const Option_t *) const
302 // Print the parameters used for the particle type identification
304 message = "=============== AliPHOSPIDv0 ================\n" ;
305 message += "Making PID\n" ;
306 message += "with parameters:\n" ;
307 message += " Maximal EMC - CPV distance (cm) %f\n" ;
308 AliInfo(Form( message.Data(),
311 if(fIDOptions.Contains("dis",TString::kIgnoreCase ))
312 AliInfo(Form(" dispersion cut %f", fDispersion )) ;
313 if(fIDOptions.Contains("ell",TString::kIgnoreCase ))
314 AliInfo(Form(" Eliptic cuts function: %s",
315 fFormula->GetTitle() )) ;
316 if(fIDOptions.Contains("tim",TString::kIgnoreCase ))
317 AliInfo(Form(" Time Gate used: %f", fTimeGate)) ;
320 //____________________________________________________________________________
321 void AliPHOSPIDv0::SetShowerProfileCut(const char * formula)
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.
327 fFormula = new TFormula("Lambda Cut",formula) ;
330 //____________________________________________________________________________
331 void AliPHOSPIDv0::PlotDispersionCuts()const
333 // produces a plot of the dispersion cut
334 TCanvas* lambdas = new TCanvas("lambdas","Cuts on the ellipse axis",200,10,700,500);
336 if(fIDOptions.Contains("ell",TString::kIgnoreCase ) ){
337 TF2 * ell = new TF2("Elliptic Cuts",fFormula->GetName(),0,3,0,3) ;
338 ell->SetMinimum(0.0000001) ;
339 ell->SetMaximum(0.001) ;
340 ell->SetLineStyle(1) ;
341 ell->SetLineWidth(2) ;
345 if( fIDOptions.Contains("dis",TString::kIgnoreCase ) ){
346 TF2 * dsp = new TF2("dispersion","(y<x)*(x*x+y*y < [0]*[0])",0,3,0,3) ;
347 dsp->SetParameter(0,fDispersion) ;
348 dsp->SetMinimum(0.0000001) ;
349 dsp->SetMaximum(0.001) ;
350 dsp->SetLineStyle(1) ;
351 dsp->SetLineColor(2) ;
352 dsp->SetLineWidth(2) ;
355 if(fIDOptions.Contains("ell",TString::kIgnoreCase ) )
363 //____________________________________________________________________________
364 TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * )const
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
372 TVector3 dir(0,0,0) ;
374 TVector3 emcglobalpos ;
377 emc->GetGlobalPosition(emcglobalpos, dummy) ;
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 ;
392 // else { // looks like a neutral
393 // dir = emcglobalpos ;
397 dir.SetZ( -dir.Z() ) ; // why ?
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 ;
411 //____________________________________________________________________________
412 void AliPHOSPIDv0::PrintRecParticles(Option_t * option)
414 // Print table of reconstructed particles
417 message = "Found %d RecParticles\n" ;
418 AliInfo(Form(message.Data(),
419 fRecParticles->GetEntriesFast() )) ;
421 if(strstr(option,"all")) { // printing found TS
422 AliInfo(" PARTICLE Index" ) ;
425 for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
426 AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;
428 Text_t particle[100];
429 switch(rp->GetType()) {
430 case AliPHOSFastRecParticle::kNEUTRALEMFAST:
431 strncpy( particle, "NEUTRAL EM FAST",100);
433 case AliPHOSFastRecParticle::kNEUTRALHAFAST:
434 strncpy(particle, "NEUTRAL HA FAST",100);
436 case AliPHOSFastRecParticle::kNEUTRALEMSLOW:
437 strncpy(particle, "NEUTRAL EM SLOW",100);
439 case AliPHOSFastRecParticle::kNEUTRALHASLOW:
440 strncpy(particle, "NEUTRAL HA SLOW",100);
442 case AliPHOSFastRecParticle::kCHARGEDEMFAST:
443 strncpy(particle, "CHARGED EM FAST",100);
445 case AliPHOSFastRecParticle::kCHARGEDHAFAST:
446 strncpy(particle, "CHARGED HA FAST",100);
448 case AliPHOSFastRecParticle::kCHARGEDEMSLOW:
449 strncpy(particle, "CHARGEDEMSLOW",100);
451 case AliPHOSFastRecParticle::kCHARGEDHASLOW:
452 strncpy(particle, "CHARGED HA SLOW",100);
456 // Int_t * primaries;
458 // primaries = rp->GetPrimaries(nprimaries);
460 AliInfo(Form(" %s %d",
461 particle, rp->GetIndexInList())) ;