git://git.uio.no / u/mrichter/AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Unused DAs removed.
[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);
3f7dbdb7 140 fIDOptions = rhs.fIDOptions;
3f7dbdb7 141 fClusterizer = rhs.fClusterizer;
142 fTSMaker = rhs.fTSMaker;
143 fFormula = rhs.fFormula;
144 fDispersion = rhs.fDispersion;
145 fCpvEmcDistance = rhs.fCpvEmcDistance;
146 fTimeGate = rhs.fTimeGate;
3f7dbdb7 147
148 return *this;
b91d88dc 149}
150
151//____________________________________________________________________________
152AliPHOSPIDv0::~AliPHOSPIDv0()
3f7dbdb7 153{
154 //Empty dtor, fFormula leaks
b91d88dc 155}
156
26aa7e4a 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//}
b91d88dc 185
186//____________________________________________________________________________
9a2cdbdf 187void AliPHOSPIDv0::TrackSegments2RecParticles(Option_t * option)
b91d88dc 188{
189 //Steering method
9a2cdbdf 190
b91d88dc 191 if(strstr(option,"tim"))
192 gBenchmark->Start("PHOSPID");
193
194 if(strstr(option,"print")) {
88cb7938 195 Print() ;
b91d88dc 196 return ;
197 }
198
9a2cdbdf 199 AliInfo(Form("%d emc clusters, %d track segments",
200 fEMCRecPoints->GetEntriesFast(),
201 fTrackSegments->GetEntriesFast())) ;
fbf811ec 202
9a2cdbdf 203 MakeRecParticles() ;
b91d88dc 204
9a2cdbdf 205 if(strstr(option,"deb"))
206 PrintRecParticles(option) ;
b91d88dc 207
b91d88dc 208 if(strstr(option,"tim")){
209 gBenchmark->Stop("PHOSPID");
9a2cdbdf 210 AliInfo(Form("took %f seconds for PID",
211 gBenchmark->GetCpuTime("PHOSPID")));
b91d88dc 212 }
213
214}
b91d88dc 215
216//____________________________________________________________________________
7f78a025 217void AliPHOSPIDv0::MakeRecParticles()
218{
219 // Reconstructs the particles from the tracksegments
b91d88dc 220
9a2cdbdf 221 fRecParticles->Clear();
88cb7938 222
9a2cdbdf 223 TIter next(fTrackSegments) ;
b91d88dc 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
9a2cdbdf 234 new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
235 rp = (AliPHOSRecParticle *)fRecParticles->At(index) ;
a278df55 236 rp->SetTrackSegment(index) ;
b91d88dc 237 rp->SetIndexInList(index) ;
238
239 AliPHOSEmcRecPoint * emc = 0 ;
240 if(ts->GetEmcIndex()>=0)
9a2cdbdf 241 emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
b91d88dc 242
243 AliPHOSRecPoint * cpv = 0 ;
244 if(ts->GetCpvIndex()>=0)
9a2cdbdf 245 cpv = (AliPHOSRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ;
b91d88dc 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)
26aa7e4a 277 if(ts->GetCpvDistance("R") < fCpvEmcDistance)
b91d88dc 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//____________________________________________________________________________
702ab87e 294void AliPHOSPIDv0:: Print(const Option_t *) const
b91d88dc 295{
296 // Print the parameters used for the particle type identification
21cd0c07 297 TString message ;
88cb7938 298 message = "=============== AliPHOSPIDv0 ================\n" ;
21cd0c07 299 message += "Making PID\n" ;
21cd0c07 300 message += "with parameters:\n" ;
301 message += " Maximal EMC - CPV distance (cm) %f\n" ;
351dd634 302 AliInfo(Form( message.Data(),
351dd634 303 fCpvEmcDistance ));
21cd0c07 304
305 if(fIDOptions.Contains("dis",TString::kIgnoreCase ))
351dd634 306 AliInfo(Form(" dispersion cut %f", fDispersion )) ;
21cd0c07 307 if(fIDOptions.Contains("ell",TString::kIgnoreCase ))
351dd634 308 AliInfo(Form(" Eliptic cuts function: %s",
309 fFormula->GetTitle() )) ;
21cd0c07 310 if(fIDOptions.Contains("tim",TString::kIgnoreCase ))
351dd634 311 AliInfo(Form(" Time Gate used: %f", fTimeGate)) ;
b91d88dc 312}
313
314//____________________________________________________________________________
8e8eae84 315void AliPHOSPIDv0::SetShowerProfileCut(const char * formula)
b91d88dc 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}
b91d88dc 323
324//____________________________________________________________________________
325void 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//____________________________________________________________________________
8f2a3661 358TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * )const
b91d88dc 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
9a2cdbdf 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.
b91d88dc 397 //account correction to the position of IP
9a2cdbdf 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 ;
b91d88dc 402
403 return dir ;
404}
405//____________________________________________________________________________
406void AliPHOSPIDv0::PrintRecParticles(Option_t * option)
407{
408 // Print table of reconstructed particles
409
21cd0c07 410 TString message ;
9a2cdbdf 411 message = "Found %d RecParticles\n" ;
351dd634 412 AliInfo(Form(message.Data(),
9a2cdbdf 413 fRecParticles->GetEntriesFast() )) ;
21cd0c07 414
b91d88dc 415 if(strstr(option,"all")) { // printing found TS
351dd634 416 AliInfo(" PARTICLE Index" ) ;
88cb7938 417
b91d88dc 418 Int_t index ;
9a2cdbdf 419 for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
420 AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;
b91d88dc 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
351dd634 454 AliInfo(Form(" %s %d",
455 particle, rp->GetIndexInList())) ;
b91d88dc 456 }
21cd0c07 457 }
b91d88dc 458}