Updates for the ROOT trunk
[u/mrichter/AliRoot.git] / PHOS / AliPHOSPIDv0.cxx
CommitLineData
b91d88dc 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
702ab87e 18/* History of cvs commits:
19 *
20 * $Log$
9a2cdbdf 21 * Revision 1.15 2007/03/06 06:57:05 kharlov
22 * DP:calculation of distance to CPV done in TSM
23 *
26aa7e4a 24 * Revision 1.14 2006/09/07 18:31:08 kharlov
25 * Effective c++ corrections (T.Pocheptsov)
26 *
3f7dbdb7 27 * Revision 1.13 2005/05/28 14:19:04 schutz
28 * Compilation warnings fixed by T.P.
29 *
702ab87e 30 */
31
b91d88dc 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 ---
b91d88dc 74#include "TTree.h"
b91d88dc 75#include "TF2.h"
76#include "TFormula.h"
77#include "TCanvas.h"
9a2cdbdf 78#include "TClonesArray.h"
7f78a025 79
b91d88dc 80#include "TBenchmark.h"
81// --- Standard library ---
82
b91d88dc 83// --- AliRoot header files ---
351dd634 84#include "AliLog.h"
b91d88dc 85#include "AliPHOSPIDv0.h"
e957fea8 86#include "AliPHOSEmcRecPoint.h"
b91d88dc 87#include "AliPHOSTrackSegment.h"
b91d88dc 88#include "AliPHOSRecParticle.h"
89#include "AliPHOSGeometry.h"
b91d88dc 90
91ClassImp( AliPHOSPIDv0)
92
93//____________________________________________________________________________
3f7dbdb7 94AliPHOSPIDv0::AliPHOSPIDv0():
3f7dbdb7 95 fIDOptions("dis time"),
3f7dbdb7 96 fClusterizer(0),
97 fTSMaker(0),
98 fFormula(0),
99 fDispersion(0.f),
100 fCpvEmcDistance(0.f),
9a2cdbdf 101 fTimeGate(2.e-9f)
b91d88dc 102{
103 // default ctor
b91d88dc 104}
105
106//____________________________________________________________________________
9a2cdbdf 107AliPHOSPIDv0::AliPHOSPIDv0(AliPHOSGeometry *geom) :
108 AliPHOSPID(geom),
3f7dbdb7 109 fIDOptions("dis time"),
3f7dbdb7 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),
9a2cdbdf 115 fTimeGate(2.e-9f)
b91d88dc 116{
9a2cdbdf 117 //ctor
3f7dbdb7 118}
b91d88dc 119
3f7dbdb7 120//____________________________________________________________________________
121AliPHOSPIDv0::AliPHOSPIDv0(const AliPHOSPIDv0 & rhs) :
122 AliPHOSPID(rhs),
3f7dbdb7 123 fIDOptions(rhs.fIDOptions),
3f7dbdb7 124 fClusterizer(rhs.fClusterizer),
125 fTSMaker(rhs.fTSMaker),
126 fFormula(rhs.fFormula),
127 fDispersion(rhs.fDispersion),
128 fCpvEmcDistance(rhs.fCpvEmcDistance),
9a2cdbdf 129 fTimeGate(rhs.fTimeGate)
3f7dbdb7 130{
131 //Copy ctor, the same as compiler-generated, possibly wrong if
132 //someone implements dtor correctly.
133}
134
135//____________________________________________________________________________
136AliPHOSPIDv0 & AliPHOSPIDv0::operator = (const AliPHOSPIDv0 & rhs)
137{
138 //Copy-assignment, emulates compiler generated, possibly wrong.
139 AliPHOSPID::operator = (rhs);
929af5d7 140 if (this != &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;
148 }
149 else {
150 AliFatal("Self assignment!");
151 }
3f7dbdb7 152 return *this;
b91d88dc 153}
154
155//____________________________________________________________________________
156AliPHOSPIDv0::~AliPHOSPIDv0()
3f7dbdb7 157{
158 //Empty dtor, fFormula leaks
b91d88dc 159}
160
26aa7e4a 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//}
b91d88dc 189
190//____________________________________________________________________________
9a2cdbdf 191void AliPHOSPIDv0::TrackSegments2RecParticles(Option_t * option)
b91d88dc 192{
193 //Steering method
9a2cdbdf 194
b91d88dc 195 if(strstr(option,"tim"))
196 gBenchmark->Start("PHOSPID");
197
198 if(strstr(option,"print")) {
88cb7938 199 Print() ;
b91d88dc 200 return ;
201 }
202
9a2cdbdf 203 AliInfo(Form("%d emc clusters, %d track segments",
204 fEMCRecPoints->GetEntriesFast(),
205 fTrackSegments->GetEntriesFast())) ;
fbf811ec 206
9a2cdbdf 207 MakeRecParticles() ;
b91d88dc 208
9a2cdbdf 209 if(strstr(option,"deb"))
210 PrintRecParticles(option) ;
b91d88dc 211
b91d88dc 212 if(strstr(option,"tim")){
213 gBenchmark->Stop("PHOSPID");
9a2cdbdf 214 AliInfo(Form("took %f seconds for PID",
215 gBenchmark->GetCpuTime("PHOSPID")));
b91d88dc 216 }
217
218}
b91d88dc 219
220//____________________________________________________________________________
7f78a025 221void AliPHOSPIDv0::MakeRecParticles()
222{
223 // Reconstructs the particles from the tracksegments
b91d88dc 224
9a2cdbdf 225 fRecParticles->Clear();
88cb7938 226
9a2cdbdf 227 TIter next(fTrackSegments) ;
b91d88dc 228 AliPHOSTrackSegment * ts ;
229 Int_t index = 0 ;
230 AliPHOSRecParticle * rp ;
231
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 ) ;
235
236 while ( (ts = (AliPHOSTrackSegment *)next()) ) {
237
9a2cdbdf 238 new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
239 rp = (AliPHOSRecParticle *)fRecParticles->At(index) ;
a278df55 240 rp->SetTrackSegment(index) ;
b91d88dc 241 rp->SetIndexInList(index) ;
242
243 AliPHOSEmcRecPoint * emc = 0 ;
244 if(ts->GetEmcIndex()>=0)
9a2cdbdf 245 emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
c78764fb 246 else
247 AliFatal(Form("ts->GetEmcIndex() return illegal index %d",ts->GetEmcIndex()));
b91d88dc 248
249 AliPHOSRecPoint * cpv = 0 ;
250 if(ts->GetCpvIndex()>=0)
9a2cdbdf 251 cpv = (AliPHOSRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ;
b91d88dc 252
253 //set momentum and energy first
254 Float_t e = emc->GetEnergy() ;
255 TVector3 dir = GetMomentumDirection(emc,cpv) ;
256 dir.SetMag(e) ;
257
258 rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
259 rp->SetCalcMass(0);
260
261 //now set type (reconstructed) of the particle
262 Int_t showerprofile = 0; // 0 narrow and 1 wide
263
264 if(ellips){
265 Float_t lambda[2] ;
266 emc->GetElipsAxis(lambda) ;
267 if(fFormula->Eval(lambda[0],lambda[1]) <= 0 )
268 showerprofile = 1 ; // not narrow
269 }
270
271 if(disp)
272 if(emc->GetDispersion() > fDispersion )
273 showerprofile = 1 ; // not narrow
274
275 Int_t slow = 0 ;
276 if(time)
277 if(emc->GetTime() > fTimeGate )
278 slow = 0 ;
279
280 // Looking at the CPV detector
281 Int_t cpvdetector= 0 ; //1 hit and 0 no hit
282 if(cpv)
26aa7e4a 283 if(ts->GetCpvDistance("R") < fCpvEmcDistance)
b91d88dc 284 cpvdetector = 1 ;
285
286 Int_t type = showerprofile + 2 * slow + 4 * cpvdetector ;
287 rp->SetType(type) ;
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);
294 index++ ;
295 }
296
297}
298
299//____________________________________________________________________________
702ab87e 300void AliPHOSPIDv0:: Print(const Option_t *) const
b91d88dc 301{
302 // Print the parameters used for the particle type identification
21cd0c07 303 TString message ;
88cb7938 304 message = "=============== AliPHOSPIDv0 ================\n" ;
21cd0c07 305 message += "Making PID\n" ;
21cd0c07 306 message += "with parameters:\n" ;
307 message += " Maximal EMC - CPV distance (cm) %f\n" ;
351dd634 308 AliInfo(Form( message.Data(),
351dd634 309 fCpvEmcDistance ));
21cd0c07 310
311 if(fIDOptions.Contains("dis",TString::kIgnoreCase ))
351dd634 312 AliInfo(Form(" dispersion cut %f", fDispersion )) ;
21cd0c07 313 if(fIDOptions.Contains("ell",TString::kIgnoreCase ))
351dd634 314 AliInfo(Form(" Eliptic cuts function: %s",
315 fFormula->GetTitle() )) ;
21cd0c07 316 if(fIDOptions.Contains("tim",TString::kIgnoreCase ))
351dd634 317 AliInfo(Form(" Time Gate used: %f", fTimeGate)) ;
b91d88dc 318}
319
320//____________________________________________________________________________
8e8eae84 321void AliPHOSPIDv0::SetShowerProfileCut(const char * formula)
b91d88dc 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 if(fFormula)
326 delete fFormula;
327 fFormula = new TFormula("Lambda Cut",formula) ;
328}
b91d88dc 329
330//____________________________________________________________________________
331void AliPHOSPIDv0::PlotDispersionCuts()const
332{
333 // produces a plot of the dispersion cut
334 TCanvas* lambdas = new TCanvas("lambdas","Cuts on the ellipse axis",200,10,700,500);
335
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) ;
342 ell->Draw() ;
343 }
344
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) ;
353 dsp->SetNpx(200) ;
354 dsp->SetNpy(200) ;
355 if(fIDOptions.Contains("ell",TString::kIgnoreCase ) )
356 dsp->Draw("same") ;
357 else
358 dsp->Draw() ;
359 }
360 lambdas->Update();
361}
362
363//____________________________________________________________________________
8f2a3661 364TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * )const
b91d88dc 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 TVector3 dir(0,0,0) ;
373
374 TVector3 emcglobalpos ;
375 TMatrix dummy ;
376
377 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 dir = emcglobalpos ;
397 dir.SetZ( -dir.Z() ) ; // why ?
398 dir.SetMag(1.) ;
399
9a2cdbdf 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.
b91d88dc 403 //account correction to the position of IP
9a2cdbdf 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 ;
b91d88dc 408
409 return dir ;
410}
411//____________________________________________________________________________
412void AliPHOSPIDv0::PrintRecParticles(Option_t * option)
413{
414 // Print table of reconstructed particles
415
21cd0c07 416 TString message ;
9a2cdbdf 417 message = "Found %d RecParticles\n" ;
351dd634 418 AliInfo(Form(message.Data(),
9a2cdbdf 419 fRecParticles->GetEntriesFast() )) ;
21cd0c07 420
b91d88dc 421 if(strstr(option,"all")) { // printing found TS
351dd634 422 AliInfo(" PARTICLE Index" ) ;
88cb7938 423
b91d88dc 424 Int_t index ;
9a2cdbdf 425 for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
426 AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;
b91d88dc 427
a7307087 428 Text_t particle[100];
b91d88dc 429 switch(rp->GetType()) {
430 case AliPHOSFastRecParticle::kNEUTRALEMFAST:
14895588 431 strncpy( particle, "NEUTRAL EM FAST",100);
b91d88dc 432 break;
433 case AliPHOSFastRecParticle::kNEUTRALHAFAST:
14895588 434 strncpy(particle, "NEUTRAL HA FAST",100);
b91d88dc 435 break;
436 case AliPHOSFastRecParticle::kNEUTRALEMSLOW:
14895588 437 strncpy(particle, "NEUTRAL EM SLOW",100);
b91d88dc 438 break ;
439 case AliPHOSFastRecParticle::kNEUTRALHASLOW:
14895588 440 strncpy(particle, "NEUTRAL HA SLOW",100);
b91d88dc 441 break ;
442 case AliPHOSFastRecParticle::kCHARGEDEMFAST:
14895588 443 strncpy(particle, "CHARGED EM FAST",100);
b91d88dc 444 break ;
445 case AliPHOSFastRecParticle::kCHARGEDHAFAST:
14895588 446 strncpy(particle, "CHARGED HA FAST",100);
b91d88dc 447 break ;
448 case AliPHOSFastRecParticle::kCHARGEDEMSLOW:
14895588 449 strncpy(particle, "CHARGEDEMSLOW",100);
b91d88dc 450 break ;
451 case AliPHOSFastRecParticle::kCHARGEDHASLOW:
14895588 452 strncpy(particle, "CHARGED HA SLOW",100);
b91d88dc 453 break ;
454 }
455
456 // Int_t * primaries;
457 // Int_t nprimaries;
458 // primaries = rp->GetPrimaries(nprimaries);
459
351dd634 460 AliInfo(Form(" %s %d",
461 particle, rp->GetIndexInList())) ;
b91d88dc 462 }
21cd0c07 463 }
b91d88dc 464}