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);
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;
151 //____________________________________________________________________________
152 AliPHOSPIDv0::~AliPHOSPIDv0()
154 //Empty dtor, fFormula leaks
158 ////____________________________________________________________________________
159 //Float_t AliPHOSPIDv0::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * Axis)const
161 // // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
163 // const AliPHOSGeometry * geom = AliPHOSLoader::GetPHOSGeometry() ;
167 // emc->GetLocalPosition(vecEmc) ;
168 // cpv->GetLocalPosition(vecCpv) ;
169 // if(emc->GetPHOSMod() == cpv->GetPHOSMod()){
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();
183 // return 100000000 ;
186 //____________________________________________________________________________
187 void AliPHOSPIDv0::TrackSegments2RecParticles(Option_t * option)
191 if(strstr(option,"tim"))
192 gBenchmark->Start("PHOSPID");
194 if(strstr(option,"print")) {
199 AliInfo(Form("%d emc clusters, %d track segments",
200 fEMCRecPoints->GetEntriesFast(),
201 fTrackSegments->GetEntriesFast())) ;
205 if(strstr(option,"deb"))
206 PrintRecParticles(option) ;
208 if(strstr(option,"tim")){
209 gBenchmark->Stop("PHOSPID");
210 AliInfo(Form("took %f seconds for PID",
211 gBenchmark->GetCpuTime("PHOSPID")));
216 //____________________________________________________________________________
217 void AliPHOSPIDv0::MakeRecParticles()
219 // Reconstructs the particles from the tracksegments
221 fRecParticles->Clear();
223 TIter next(fTrackSegments) ;
224 AliPHOSTrackSegment * ts ;
226 AliPHOSRecParticle * rp ;
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 ) ;
232 while ( (ts = (AliPHOSTrackSegment *)next()) ) {
234 new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
235 rp = (AliPHOSRecParticle *)fRecParticles->At(index) ;
236 rp->SetTrackSegment(index) ;
237 rp->SetIndexInList(index) ;
239 AliPHOSEmcRecPoint * emc = 0 ;
240 if(ts->GetEmcIndex()>=0)
241 emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
243 AliFatal(Form("ts->GetEmcIndex() return illegal index %d",ts->GetEmcIndex()));
245 AliPHOSRecPoint * cpv = 0 ;
246 if(ts->GetCpvIndex()>=0)
247 cpv = (AliPHOSRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ;
249 //set momentum and energy first
250 Float_t e = emc->GetEnergy() ;
251 TVector3 dir = GetMomentumDirection(emc,cpv) ;
254 rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
257 //now set type (reconstructed) of the particle
258 Int_t showerprofile = 0; // 0 narrow and 1 wide
262 emc->GetElipsAxis(lambda) ;
263 if(fFormula->Eval(lambda[0],lambda[1]) <= 0 )
264 showerprofile = 1 ; // not narrow
268 if(emc->GetDispersion() > fDispersion )
269 showerprofile = 1 ; // not narrow
273 if(emc->GetTime() > fTimeGate )
276 // Looking at the CPV detector
277 Int_t cpvdetector= 0 ; //1 hit and 0 no hit
279 if(ts->GetCpvDistance("R") < fCpvEmcDistance)
282 Int_t type = showerprofile + 2 * slow + 4 * cpvdetector ;
284 rp->SetProductionVertex(0,0,0,0);
285 rp->SetFirstMother(-1);
286 rp->SetLastMother(-1);
287 rp->SetFirstDaughter(-1);
288 rp->SetLastDaughter(-1);
289 rp->SetPolarisation(0,0,0);
295 //____________________________________________________________________________
296 void AliPHOSPIDv0:: Print(const Option_t *) const
298 // Print the parameters used for the particle type identification
300 message = "=============== AliPHOSPIDv0 ================\n" ;
301 message += "Making PID\n" ;
302 message += "with parameters:\n" ;
303 message += " Maximal EMC - CPV distance (cm) %f\n" ;
304 AliInfo(Form( message.Data(),
307 if(fIDOptions.Contains("dis",TString::kIgnoreCase ))
308 AliInfo(Form(" dispersion cut %f", fDispersion )) ;
309 if(fIDOptions.Contains("ell",TString::kIgnoreCase ))
310 AliInfo(Form(" Eliptic cuts function: %s",
311 fFormula->GetTitle() )) ;
312 if(fIDOptions.Contains("tim",TString::kIgnoreCase ))
313 AliInfo(Form(" Time Gate used: %f", fTimeGate)) ;
316 //____________________________________________________________________________
317 void AliPHOSPIDv0::SetShowerProfileCut(const char * formula)
319 //set shape of the cut on the axis of ellipce, drown around shouer
320 //shower considered "narrow" if Formula(lambda[0],lambda[1]) > 0.
323 fFormula = new TFormula("Lambda Cut",formula) ;
326 //____________________________________________________________________________
327 void AliPHOSPIDv0::PlotDispersionCuts()const
329 // produces a plot of the dispersion cut
330 TCanvas* lambdas = new TCanvas("lambdas","Cuts on the ellipse axis",200,10,700,500);
332 if(fIDOptions.Contains("ell",TString::kIgnoreCase ) ){
333 TF2 * ell = new TF2("Elliptic Cuts",fFormula->GetName(),0,3,0,3) ;
334 ell->SetMinimum(0.0000001) ;
335 ell->SetMaximum(0.001) ;
336 ell->SetLineStyle(1) ;
337 ell->SetLineWidth(2) ;
341 if( fIDOptions.Contains("dis",TString::kIgnoreCase ) ){
342 TF2 * dsp = new TF2("dispersion","(y<x)*(x*x+y*y < [0]*[0])",0,3,0,3) ;
343 dsp->SetParameter(0,fDispersion) ;
344 dsp->SetMinimum(0.0000001) ;
345 dsp->SetMaximum(0.001) ;
346 dsp->SetLineStyle(1) ;
347 dsp->SetLineColor(2) ;
348 dsp->SetLineWidth(2) ;
351 if(fIDOptions.Contains("ell",TString::kIgnoreCase ) )
359 //____________________________________________________________________________
360 TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * )const
362 // Calculates the momentum direction:
363 // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint
364 // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints
365 // However because of the poor position resolution of PPSD the direction is always taken as if we were
368 TVector3 dir(0,0,0) ;
370 TVector3 emcglobalpos ;
373 emc->GetGlobalPosition(emcglobalpos, dummy) ;
376 // The following commented code becomes valid once the PPSD provides
377 // a reasonable position resolution, at least as good as EMC !
378 // TVector3 ppsdlglobalpos ;
379 // TVector3 ppsduglobalpos ;
380 // if( fPpsdLowRecPoint ){ // certainly a photon that has concerted
381 // fPpsdLowRecPoint->GetGlobalPosition(ppsdlglobalpos, mdummy) ;
382 // dir = emcglobalpos - ppsdlglobalpos ;
383 // if( fPpsdUpRecPoint ){ // not looks like a charged
384 // fPpsdUpRecPoint->GetGlobalPosition(ppsduglobalpos, mdummy) ;
385 // dir = ( dir + emcglobalpos - ppsduglobalpos ) * 0.5 ;
388 // else { // looks like a neutral
389 // dir = emcglobalpos ;
393 dir.SetZ( -dir.Z() ) ; // why ?
396 // One can not access MC information in the reconstruction!!
397 // PLEASE FIT IT, EITHER BY TAKING 0,0,0 OR ACCESSING THE
398 // VERTEX DIAMOND FROM CDB GRP FOLDER.
399 //account correction to the position of IP
400 // Float_t xo,yo,zo ; //Coordinates of the origin
401 // gAlice->Generator()->GetOrigin(xo,yo,zo) ;
402 // TVector3 origin(xo,yo,zo);
403 // dir = dir - origin ;
407 //____________________________________________________________________________
408 void AliPHOSPIDv0::PrintRecParticles(Option_t * option)
410 // Print table of reconstructed particles
413 message = "Found %d RecParticles\n" ;
414 AliInfo(Form(message.Data(),
415 fRecParticles->GetEntriesFast() )) ;
417 if(strstr(option,"all")) { // printing found TS
418 AliInfo(" PARTICLE Index" ) ;
421 for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
422 AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;
424 Text_t particle[100];
425 switch(rp->GetType()) {
426 case AliPHOSFastRecParticle::kNEUTRALEMFAST:
427 strncpy( particle, "NEUTRAL EM FAST",100);
429 case AliPHOSFastRecParticle::kNEUTRALHAFAST:
430 strncpy(particle, "NEUTRAL HA FAST",100);
432 case AliPHOSFastRecParticle::kNEUTRALEMSLOW:
433 strncpy(particle, "NEUTRAL EM SLOW",100);
435 case AliPHOSFastRecParticle::kNEUTRALHASLOW:
436 strncpy(particle, "NEUTRAL HA SLOW",100);
438 case AliPHOSFastRecParticle::kCHARGEDEMFAST:
439 strncpy(particle, "CHARGED EM FAST",100);
441 case AliPHOSFastRecParticle::kCHARGEDHAFAST:
442 strncpy(particle, "CHARGED HA FAST",100);
444 case AliPHOSFastRecParticle::kCHARGEDEMSLOW:
445 strncpy(particle, "CHARGEDEMSLOW",100);
447 case AliPHOSFastRecParticle::kCHARGEDHASLOW:
448 strncpy(particle, "CHARGED HA SLOW",100);
452 // Int_t * primaries;
454 // primaries = rp->GetPrimaries(nprimaries);
456 AliInfo(Form(" %s %d",
457 particle, rp->GetIndexInList())) ;