]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PHOS/AliPHOSPIDv1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Correcter eror by A. Dainese
[u/mrichter/AliRoot.git] / PHOS / AliPHOSPIDv1.cxx
CommitLineData
6ad0bfa0 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
b2a60966 16/* $Id$ */
17
6ad0bfa0 18//_________________________________________________________________________
b2a60966 19// Implementation version v1 of the PHOS particle identifier
7acf6008 20// Particle identification based on the
148b2bba 21// - RCPV: distance from CPV recpoint to EMCA recpoint.
22// - TOF
23// - PCA: Principal Components Analysis..
24// The identified particle has an identification number corresponding
25// to a 9 bits number:
26// -Bit 0 to 2: bit set if RCPV > fCpvEmcDistance (each bit corresponds
27// to a different efficiency-purity point of the photon identification)
28// -Bit 3 to 5: bit set if TOF < fTimeGate (each bit corresponds
29// to a different efficiency-purity point of the photon identification)
30// -Bit 6 to 9: bit set if Principal Components are
31// inside an ellipse defined by fX_center, fY_center, fA, fB, fAngle
32// (each bit corresponds to a different efficiency-purity point of the
33// photon identification)
a4e98857 34//
7acf6008 35//
a4e98857 36// use case:
37// root [0] AliPHOSPIDv1 * p1 = new AliPHOSPIDv1("galice.root")
38// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
39// root [1] p1->SetIdentificationMethod("disp ellipse")
40// root [2] p1->ExecuteTask()
148b2bba 41// root [3] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","v1")
a4e98857 42// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
148b2bba 43// // reading headers from file galice1.root and create RecParticles with title v1
44 // TrackSegments and RecPoints with title "v1" are used
7acf6008 45// // set file name for the branch RecParticles
148b2bba 46// root [4] p2->ExecuteTask("deb all time")
7acf6008 47// // available options
48// // "deb" - prints # of reconstructed particles
49// // "deb all" - prints # and list of RecParticles
50// // "time" - prints benchmarking results
51//
148b2bba 52// root [5] AliPHOSPIDv1 * p3 = new AliPHOSPIDv1("galice1.root","v1","v0")
53// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
54// // reading headers from file galice1.root and create RecParticles with title v1
55 // RecPoints and TrackSegments with title "v0" are used
56// root [6] p3->ExecuteTask()
7acf6008 57//*-- Author: Yves Schutz (SUBATECH) & Gines Martinez (SUBATECH) &
148b2bba 58// Gustavo Conesa April 2002
6ad0bfa0 59
60// --- ROOT system ---
7acf6008 61#include "TROOT.h"
62#include "TTree.h"
63#include "TFile.h"
64#include "TF2.h"
65#include "TFormula.h"
66#include "TCanvas.h"
67#include "TFolder.h"
68#include "TSystem.h"
69#include "TBenchmark.h"
148b2bba 70#include "TMatrixD.h"
71#include "TPrincipal.h"
72#include "TSystem.h"
73
6ad0bfa0 74// --- Standard library ---
75
de9ec31b 76#include <iostream.h>
148b2bba 77#include <fstream.h>
7acf6008 78#include <iomanip.h>
6ad0bfa0 79
80// --- AliRoot header files ---
81
7acf6008 82#include "AliRun.h"
83#include "AliGenerator.h"
84#include "AliPHOS.h"
26d4b141 85#include "AliPHOSPIDv1.h"
7acf6008 86#include "AliPHOSClusterizerv1.h"
6ad0bfa0 87#include "AliPHOSTrackSegment.h"
7acf6008 88#include "AliPHOSTrackSegmentMakerv1.h"
6ad0bfa0 89#include "AliPHOSRecParticle.h"
7b7c1533 90#include "AliPHOSGeometry.h"
91#include "AliPHOSGetter.h"
6ad0bfa0 92
26d4b141 93ClassImp( AliPHOSPIDv1)
6ad0bfa0 94
1cb7c1ee 95//____________________________________________________________________________
96AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID()
97{
a4e98857 98 // default ctor
148b2bba 99
7b7c1533 100 fHeaderFileName = "" ;
101 fTrackSegmentsTitle= "" ;
102 fRecPointsTitle = "" ;
103 fRecParticlesTitle = "" ;
2685bf00 104 fRecParticlesInRun = 0 ;
148b2bba 105 fClusterizer = 0 ;
106 fTSMaker = 0 ;
107 fFrom = "" ;
108
7acf6008 109}
110
111//____________________________________________________________________________
148b2bba 112AliPHOSPIDv1::AliPHOSPIDv1(const char * headerFile,const char * name, const char * from) : AliPHOSPID(headerFile, name)
7acf6008 113{
a4e98857 114 //ctor with the indication on where to look for the track segments
7b7c1533 115
116 fHeaderFileName = GetTitle() ;
117 fTrackSegmentsTitle = GetName() ;
118 fRecPointsTitle = GetName() ;
148b2bba 119 fRecParticlesTitle = GetName() ;
7b7c1533 120 TString tempo(GetName()) ;
bf8f1fbd 121 tempo.Append(":") ;
7b7c1533 122 tempo.Append(Version()) ;
bf8f1fbd 123 SetName(tempo) ;
2b60655b 124 fRecParticlesInRun = 0 ;
148b2bba 125 if ( from == 0 )
126 fFrom = name ;
127 else
128 fFrom = from ;
2bd5457f 129 Init() ;
7acf6008 130
131}
7b7c1533 132
7acf6008 133//____________________________________________________________________________
134AliPHOSPIDv1::~AliPHOSPIDv1()
135{
148b2bba 136 delete [] fX ; // Principal input
137 delete [] fP ; // Principal components
138 delete fParameters ; // Matrix of Parameters
7acf6008 139}
2bd5457f 140
148b2bba 141//____________________________________________________________________________
142void AliPHOSPIDv1::Init()
143{
144 // Make all memory allocations that are not possible in default constructor
145 // Add the PID task to the list of PHOS tasks
146
147 if ( strcmp(GetTitle(), "") == 0 )
148 SetTitle("galice.root") ;
149
150 TString taskName(GetName()) ;
151 taskName.Remove(taskName.Index(Version())-1) ;
152
153 // PCA : To do the Principal Components Analysis it is necessary
154 // the Principal file, which is opened here
155 fX = new double[7]; // Data for the PCA
156 fP = new double[7]; // Eigenvalues of the PCA
157 fFileName = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ;
158 TFile f( fFileName.Data(), "read" ) ;
159 fPrincipal = dynamic_cast<TPrincipal*> (f.Get("principal")) ;
160 f.Close() ;
161
162 // Initialization of the Parameters matrix. In the File Parameters.dat
163 // are all the parameters. These are introduced in a matrix of 21x3 elements.
164 // All the parameters defined in this file are, in order of row (there are
165 // 3 rows per parameter): CpvtoEmcDistanceCut, TimeGate (and the ellipse
166 // parameters), X_center, Y_center, a, b, angle. Each row of a given parameter
167 // depends on the cluster energy range (0.3-1,1-2, >2 )
168 // Each column designes the parameters for a point in the Efficiency-Purity
169 // of the photon identification P1(0.959,0.625), P2(0.919,0.835), P3(0.833,0.901).
170
171 fParameters = new TMatrixD(21,3) ;
172
173 fFileNamePar = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat");
174 ifstream paramFile(fFileNamePar, ios::in) ;
175
176 Int_t i,j ;
7b7c1533 177
148b2bba 178 for(i = 0; i< 21; i++){
179 for(j = 0; j< 3; j++){
180 paramFile >> (*fParameters)(i,j) ;
181 }
182 }
183 paramFile.close();
184
185 AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), fFrom.Data()) ;
186
187 gime->SetRecParticlesTitle(taskName) ;
188 if ( gime == 0 ) {
189 cerr << "ERROR: AliPHOSPIDv1::Init -> Could not obtain the Getter object !" << endl ;
190 return ;
191 }
192
193 gime->PostPID(this) ;
194 // create a folder on the white board //YSAlice/WhiteBoard/RecParticles/PHOS/recparticlesName
195 gime->PostRecParticles(taskName.Data() ) ;
196
197}
69183710 198//____________________________________________________________________________
148b2bba 199Double_t AliPHOSPIDv1::GetCpvtoEmcDistanceCut(const Float_t Cluster_En, const TString Eff_Pur)const
200{
201 // Get CpvtoEmcDistanceCut parameter depending on the cluster energy and
202 // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
203 // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
204 // EFFICIENCY by PURITY)
205
206 Int_t cluster = -1 ;
207 Int_t eff_pur = -1 ;
208
209 // Check the cluster energy range
210 if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
211 if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
212 if( Cluster_En > 2.0) cluster = 2 ;
213
214 if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
215 if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
216 if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
217
218 if(cluster ==-1){
219 cout<<"Invalid Cluster Energy option"<<endl;
220 }
221 else if(eff_pur ==-1){
222 cout<<"Invalid Efficiency-Purity option"<<endl;
223 }
224 else{
225 return (*fParameters)(cluster,eff_pur) ;
226 }
227}
228//____________________________________________________________________________
229
230Double_t AliPHOSPIDv1::GetTimeGate(const Float_t Cluster_En, const TString Eff_Pur) const
231{
232 // Get TimeGate parameter depending on the cluster energy and
233 // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
234 // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
235 // EFFICIENCY by PURITY)
236 Int_t cluster = -1 ;
237 Int_t eff_pur = -1 ;
238
239 if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
240 if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
241 if( Cluster_En > 2.0) cluster = 2 ;
242
243 if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
244 if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
245 if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
246 if(cluster ==-1){
247 cout<<"Invalid Cluster Energy option"<<endl;
248 }
249 else if(eff_pur ==-1){
250 cout<<"Invalid Efficiency-Purity option"<<endl;
251 }
252 else{
253 return (*fParameters)(cluster+3,eff_pur) ;
254 }
255}
256//_____________________________________________________________________________
7acf6008 257Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * Axis)const
69183710 258{
259 // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
148b2bba 260
7b7c1533 261 const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
69183710 262 TVector3 vecEmc ;
7acf6008 263 TVector3 vecCpv ;
148b2bba 264 if(cpv){
265 emc->GetLocalPosition(vecEmc) ;
266 cpv->GetLocalPosition(vecCpv) ;
267 if(emc->GetPHOSMod() == cpv->GetPHOSMod()){
268 // Correct to difference in CPV and EMC position due to different distance to center.
269 // we assume, that particle moves from center
270 Float_t dCPV = geom->GetIPtoOuterCoverDistance();
271 Float_t dEMC = geom->GetIPtoCrystalSurface() ;
272 dEMC = dEMC / dCPV ;
273 vecCpv = dEMC * vecCpv - vecEmc ;
274 if (Axis == "X") return vecCpv.X();
275 if (Axis == "Y") return vecCpv.Y();
276 if (Axis == "Z") return vecCpv.Z();
277 if (Axis == "R") return vecCpv.Mag();
7acf6008 278 }
148b2bba 279
280 return 100000000 ;
281 }
7acf6008 282 return 100000000 ;
69183710 283}
284
6ad0bfa0 285//____________________________________________________________________________
148b2bba 286Int_t AliPHOSPIDv1::GetPrincipalSign(Double_t* P, Int_t cluster, Int_t eff_pur )const
287{
288 //This method gives if the PCA of the particle are inside a defined ellipse
289
290 // Get the parameters that define the ellipse stored in the
291 // fParameters matrix.
292 Double_t X_center = (*fParameters)(cluster+6,eff_pur) ;
293 Double_t Y_center = (*fParameters)(cluster+9,eff_pur) ;
294 Double_t A = (*fParameters)(cluster+12,eff_pur) ;
295 Double_t B = (*fParameters)(cluster+15,eff_pur) ;
296 Double_t Angle = (*fParameters)(cluster+18,eff_pur) ;
297
298 Int_t prinsign;
299 Double_t Dx = 0. ;
300 Double_t Delta = 0. ;
301 Double_t Y = 0. ;
302 Double_t Y_1 = 0. ;
303 Double_t Y_2 = 0. ;
304 Double_t Pi = TMath::Pi() ;
305 Double_t Cos_Theta = TMath::Cos(Pi*Angle/180.) ;
306 Double_t Sin_Theta = TMath::Sin(Pi*Angle/180.) ;
307
308 Dx = P[0] - X_center ;
309 Delta = 4.*A*A*A*B* (A*A*Cos_Theta*Cos_Theta
310 + B*B*Sin_Theta*Sin_Theta - Dx*Dx) ;
311 if (Delta < 0.)
312 {prinsign=0;}
313
314 else if (Delta == 0.)
315 {
316 Y = Cos_Theta*Sin_Theta*(A*A - B*B)*Dx /
317 (A*A*Cos_Theta*Cos_Theta + B*B*Sin_Theta*Sin_Theta) ;
318 Y += Y_center ;
319 if(P[1]==Y )
320 {prinsign=1;}
321 else
322 {prinsign=0;}
323 }
324 else
325 {
326 Y_1 = (Cos_Theta*Sin_Theta*(A*A - B*B) *Dx +
327 TMath::Sqrt(Delta)/2.)/(A*A*Cos_Theta*Cos_Theta +
328 B*B*Sin_Theta*Sin_Theta) ;
329 Y_2 = (Cos_Theta*Sin_Theta*(A*A - B*B) *Dx -
330 TMath::Sqrt(Delta)/2.)/(A*A*Cos_Theta*Cos_Theta
331 + B*B*Sin_Theta*Sin_Theta) ;
332 Y_1 += Y_center ;
333 Y_2 += Y_center ;
334 if ((P[1]<=Y_1) && (P[1]>=Y_2))
335 {prinsign=1;}
336 else
337 {prinsign=0;}
338 }
339 return prinsign;
340}
341
342//____________________________________________________________________________
343 void AliPHOSPIDv1::SetEllipseParameters(Float_t Cluster_En, TString Eff_Pur, Float_t x, Float_t y,Float_t a, Float_t b,Float_t angle)
344{
345
346 // Set all ellipse parameters depending on the cluster energy and
347 // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
348 // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
349 // EFFICIENCY by PURITY)
350
351 Int_t cluster = -1 ;
352 Int_t eff_pur = -1 ;
353
354 if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
355 if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
356 if( Cluster_En > 2.0) cluster= 2 ;
357
358 if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
359 if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
360 if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
361
362 if(cluster ==-1){
363 cout<<"Invalid Cluster Energy option"<<endl;
364 }
365 else if(eff_pur ==-1){
366 cout<<"Invalid Efficiency-Purity option"<<endl;
367 }
368 else{
369 (*fParameters)(cluster+6,eff_pur) = x ;
370 (*fParameters)(cluster+9,eff_pur) = y ;
371 (*fParameters)(cluster+12,eff_pur) = a ;
372 (*fParameters)(cluster+15,eff_pur) = b ;
373 (*fParameters)(cluster+18,eff_pur) = angle ;
374 }
375
376}
377//__________________________________________________________________________
378void AliPHOSPIDv1::SetEllipseXCenter(Float_t Cluster_En, TString Eff_Pur, Float_t x)
379{
380 // Set the ellipse parameter x_center depending on the custer energy and
381 // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
382 // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
383 // EFFICIENCY by PURITY)
384
385 Int_t cluster = -1 ;
386 Int_t eff_pur = -1 ;
387
388 if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
389 if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
390 if( Cluster_En > 2.0) cluster = 2 ;
391
392 if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
393 if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
394 if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
395
396 if(cluster ==-1){
397 cout<<"Invalid Cluster Energy option"<<endl;
398 }
399 else if(eff_pur ==-1){
400 cout<<"Invalid Efficiency-Purity option"<<endl;
401 }
402 else{
403 (*fParameters)(cluster+6,eff_pur) = x ;
404 }
405}
406//_________________________________________________________________________
407void AliPHOSPIDv1::SetEllipseYCenter(Float_t Cluster_En, TString Eff_Pur, Float_t y)
408{
409
410 // Set the ellipse parameter y_center depending on the cluster energy and
411 // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
412 // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
413 // EFFICIENCY by PURITY)
414
415 Int_t cluster = -1 ;
416 Int_t eff_pur = -1 ;
417
418 if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
419 if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
420 if( Cluster_En > 2.0) cluster = 2 ;
421
422 if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
423 if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
424 if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
425
426 if(cluster ==-1){
427 cout<<"Invalid Cluster Energy option"<<endl;
428 }
429 else if(eff_pur ==-1){
430 cout<<"Invalid Efficiency-Purity option"<<endl ;
431 }
432 else{
433 (*fParameters)(cluster+9,eff_pur) = y ;
434 }
435}
436//_________________________________________________________________________
437void AliPHOSPIDv1::SetEllipseAParameter(Float_t Cluster_En, TString Eff_Pur, Float_t a)
438{
439 // Set the ellipse parameter a depending on the cluster energy and
440 // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
441 // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
442 // EFFICIENCY by PURITY)
443
444 Int_t cluster = -1 ;
445 Int_t eff_pur = -1 ;
446
447 if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
448 if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
449 if( Cluster_En > 2.0) cluster = 2 ;
450
451 if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
452 if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
453 if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
454
455 if(cluster ==-1){
456 cout<<"Invalid Cluster Energy option"<<endl;
457 }
458 else if(eff_pur ==-1){
459 cout<<"Invalid Efficiency-Purity option"<<endl;
460 }
461 else{
462 (*fParameters)(cluster+12,eff_pur) = a ;
463 }
464}
465//________________________________________________________________________
466void AliPHOSPIDv1::SetEllipseBParameter(Float_t Cluster_En, TString Eff_Pur, Float_t b)
467{
468 // Set the ellipse parameter b depending on the cluster energy and
469 // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
470 // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
471 // EFFICIENCY by PURITY)
472
473 Int_t cluster = -1 ;
474 Int_t eff_pur = -1 ;
475
476 if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
477 if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
478 if( Cluster_En > 2.0) cluster = 2 ;
479
480 if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
481 if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
482 if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
483
484 if(cluster ==-1){
485 cout<<"Invalid Cluster Energy option"<<endl;
486 }
487 else if(eff_pur ==-1){
488 cout<<"Invalid Efficiency-Purity option"<<endl;
489 }
490 else{
491 (*fParameters)(cluster+15,eff_pur) = b ;
492 }
493}
494//________________________________________________________________________
495void AliPHOSPIDv1::SetEllipseAngle(Float_t Cluster_En, TString Eff_Pur, Float_t angle)
496{
497
498 // Set the ellipse parameter angle depending on the cluster energy and
499 // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
500 // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
501 // EFFICIENCY by PURITY)
502
503 Int_t cluster = -1 ;
504 Int_t eff_pur = -1 ;
505
506 if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
507 if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
508 if( Cluster_En > 2.0) cluster = 2 ;
509
510 if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
511 if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
512 if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
513
514 if(cluster ==-1){
515 cout<<"Invalid Cluster Energy option"<<endl;
516 }
517 else if(eff_pur ==-1){
518 cout<<"Invalid Efficiency-Purity option"<<endl;
519 }
520 else{
521 (*fParameters)(cluster+18,eff_pur) = angle ;
522 }
523}
524//_____________________________________________________________________________
525void AliPHOSPIDv1::SetCpvtoEmcDistanceCut(Float_t Cluster_En, TString Eff_Pur, Float_t cut)
526{
527
528 // Set the parameter Cpvto EmcDistanceCut depending on the cluster energy and
529 // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
530 // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
531 // EFFICIENCY by PURITY)
532
533 Int_t cluster = -1 ;
534 Int_t eff_pur = -1 ;
535
536 if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
537 if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
538 if( Cluster_En > 2.0) cluster = 2 ;
539
540 if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
541 if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
542 if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
543
544 if(cluster ==-1){
545 cout<<"Invalid Cluster Energy option"<<endl;
546 }
547 else if(eff_pur ==-1){
548 cout<<"Invalid Efficiency-Purity option"<<endl;
549 }
550 else{
551 (*fParameters)(cluster,eff_pur) = cut ;
552 }
553}
554//_____________________________________________________________________________
555void AliPHOSPIDv1::SetTimeGate(Float_t Cluster_En, TString Eff_Pur, Float_t gate)
556{
557
558 // Set the parameter TimeGate depending on the cluster energy and
559 // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
560 // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
561 // EFFICIENCY by PURITY)
562
563 Int_t cluster = -1 ;
564 Int_t eff_pur = -1 ;
565
566 if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
567 if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
568 if( Cluster_En > 2.0) cluster = 2 ;
569
570 if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
571 if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
572 if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
573
574 if(cluster ==-1){
575 cout<<"Invalid Cluster Energy option"<<endl;
576 }
577 else if(eff_pur ==-1){
578 cout<<"Invalid Efficiency-Purity option"<<endl;
579 }
580 else{
581 (*fParameters)(cluster+3,eff_pur) = gate ;
582 }
583}
584//____________________________________________________________________________
585
7acf6008 586void AliPHOSPIDv1::Exec(Option_t * option)
6ad0bfa0 587{
f035f6ce 588 //Steering method
9688c1dd 589
bf8f1fbd 590 if( strcmp(GetName(), "")== 0 )
7acf6008 591 Init() ;
bf8f1fbd 592
593 if(strstr(option,"tim"))
7acf6008 594 gBenchmark->Start("PHOSPID");
bf8f1fbd 595
596 if(strstr(option,"print")) {
7b7c1533 597 Print("") ;
598 return ;
bf8f1fbd 599 }
9688c1dd 600
148b2bba 601 //cout << gDirectory->GetName() << endl ;
602
bf8f1fbd 603 gAlice->GetEvent(0) ;
148b2bba 604
7b7c1533 605 //check, if the branch with name of this" already exits?
606 TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
607 TIter next(lob) ;
608 TBranch * branch = 0 ;
609 Bool_t phospidfound = kFALSE, pidfound = kFALSE ;
610
611 TString taskName(GetName()) ;
bf8f1fbd 612 taskName.Remove(taskName.Index(Version())-1) ;
9688c1dd 613
7b7c1533 614 while ( (branch = (TBranch*)next()) && (!phospidfound || !pidfound) ) {
615 if ( (strcmp(branch->GetName(), "PHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
616 phospidfound = kTRUE ;
617
618 else if ( (strcmp(branch->GetName(), "AliPHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
619 pidfound = kTRUE ;
620 }
9688c1dd 621
7b7c1533 622 if ( phospidfound || pidfound ) {
623 cerr << "WARNING: AliPHOSPIDv1::Exec -> RecParticles and/or PIDtMaker branch with name "
624 << taskName.Data() << " already exits" << endl ;
625 return ;
626 }
bf8f1fbd 627
7b7c1533 628 Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
629 Int_t ievent ;
148b2bba 630 AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
7b7c1533 631 for(ievent = 0; ievent < nevents; ievent++){
148b2bba 632 gime->Event(ievent,"RA") ;
633
7acf6008 634 MakeRecParticles() ;
9688c1dd 635
7b7c1533 636 WriteRecParticles(ievent);
9688c1dd 637
7acf6008 638 if(strstr(option,"deb"))
639 PrintRecParticles(option) ;
94de8339 640
641 //increment the total number of rec particles per run
148b2bba 642 fRecParticlesInRun += gime->RecParticles(taskName)->GetEntriesFast() ;
94de8339 643
7acf6008 644 }
9688c1dd 645
7acf6008 646 if(strstr(option,"tim")){
647 gBenchmark->Stop("PHOSPID");
648 cout << "AliPHOSPID:" << endl ;
649 cout << " took " << gBenchmark->GetCpuTime("PHOSPID") << " seconds for PID "
7b7c1533 650 << gBenchmark->GetCpuTime("PHOSPID")/nevents << " seconds per event " << endl ;
7acf6008 651 cout << endl ;
652 }
9688c1dd 653
7b7c1533 654}
655
7acf6008 656//____________________________________________________________________________
657void AliPHOSPIDv1::MakeRecParticles(){
658
b2a60966 659 // Makes a RecParticle out of a TrackSegment
148b2bba 660
bf8f1fbd 661 TString taskName(GetName()) ;
662 taskName.Remove(taskName.Index(Version())-1) ;
148b2bba 663
7b7c1533 664 AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
148b2bba 665 TObjArray * emcRecPoints = gime->EmcRecPoints(fFrom) ;
666 TObjArray * cpvRecPoints = gime->CpvRecPoints(fFrom) ;
667 TClonesArray * trackSegments = gime->TrackSegments(fFrom) ;
668 if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) {
669 cerr << "ERROR: AliPHOSPIDv1::MakeRecParticles -> RecPoints or TrackSegments with name "
670 << fFrom << " not found ! " << endl ;
671 abort() ;
672 }
bf8f1fbd 673 TClonesArray * recParticles = gime->RecParticles(taskName) ;
01a599c9 674 recParticles->Clear();
148b2bba 675
676
7b7c1533 677 TIter next(trackSegments) ;
7acf6008 678 AliPHOSTrackSegment * ts ;
6ad0bfa0 679 Int_t index = 0 ;
09fc14a0 680 AliPHOSRecParticle * rp ;
7acf6008 681
7acf6008 682 while ( (ts = (AliPHOSTrackSegment *)next()) ) {
fad3e5b9 683
7b7c1533 684 new( (*recParticles)[index] ) AliPHOSRecParticle() ;
685 rp = (AliPHOSRecParticle *)recParticles->At(index) ;
7acf6008 686 rp->SetTraskSegment(index) ;
9688c1dd 687 rp->SetIndexInList(index) ;
148b2bba 688
7acf6008 689 AliPHOSEmcRecPoint * emc = 0 ;
690 if(ts->GetEmcIndex()>=0)
7b7c1533 691 emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
fad3e5b9 692
7acf6008 693 AliPHOSRecPoint * cpv = 0 ;
694 if(ts->GetCpvIndex()>=0)
7b7c1533 695 cpv = (AliPHOSRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ;
fad3e5b9 696
7acf6008 697 //set momentum and energy first
698 Float_t e = emc->GetEnergy() ;
9688c1dd 699 TVector3 dir = GetMomentumDirection(emc,cpv) ;
7acf6008 700 dir.SetMag(e) ;
701
702 rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
703 rp->SetCalcMass(0);
fad3e5b9 704
148b2bba 705 // Now set type (reconstructed) of the particle
706
707 // Choose the cluster energy range
708 Int_t cluster = 0 ; // Ellipse and rcpv cut in function of the cluster energy
709 if((e > 0.3)&&(e <= 1.0)) cluster = 0 ;
710 if((e > 1.0)&&(e <= 2.0)) cluster = 1 ;
711 if( e > 2.0) cluster = 2 ;
712
713 // Loop of Efficiency-Purity (the 3 points of purity or efficiency are taken
714 // into account to set the particle identification)
715 for(Int_t eff_pur = 0; eff_pur < 3 ; eff_pur++){
716
717 // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, 1st,
718 // 2nd or 3rd bit (depending on the efficiency-purity point )is set to 1 .
719 if(GetDistance(emc, cpv, "R") > (*fParameters)(cluster,eff_pur) )
720 rp->SetPIDBit(eff_pur) ;
721
722 // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th
723 // bit (depending on the efficiency-purity point )is set to 1
724 if(emc->GetTime()< (*fParameters)(cluster+3,eff_pur))
725 rp->SetPIDBit(eff_pur+3) ;
726
727 // Looking PCA. Define and calculate the data (X), introduce in the function
728 // X2P that gives the components (P).
729 Float_t fSpher = 0. ;
730 Float_t fEmaxdtotal = 0. ;
731 Float_t lambda[2] ;
732
733 emc->GetElipsAxis(lambda) ;
734 if((lambda[0]+lambda[1])!=0) fSpher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
fad3e5b9 735
148b2bba 736 fEmaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy();
7acf6008 737
148b2bba 738 fX[0] = lambda[0] ;
739 fX[1] = lambda[1] ;
740 fX[2] = emc->GetDispersion() ;
741 fX[3] = fSpher ;
742 fX[4] = emc->GetMultiplicity() ;
743 fX[5] = fEmaxdtotal ;
744 fX[6] = emc->GetCoreEnergy() ;
7acf6008 745
148b2bba 746 fPrincipal->X2P(fX,fP);
747
748 //If we are inside the ellipse, 7th, 8th or 9th
749 // bit (depending on the efficiency-purity point )is set to 1
750 if(GetPrincipalSign(fP,cluster,eff_pur) == 1)
751 rp->SetPIDBit(eff_pur+6) ;
752
753 }
e747b8da 754 rp->SetProductionVertex(0,0,0,0);
755 rp->SetFirstMother(-1);
756 rp->SetLastMother(-1);
757 rp->SetFirstDaughter(-1);
758 rp->SetLastDaughter(-1);
759 rp->SetPolarisation(0,0,0);
6ad0bfa0 760 index++ ;
761 }
7acf6008 762
6ad0bfa0 763}
764
09fc14a0 765//____________________________________________________________________________
7acf6008 766void AliPHOSPIDv1:: Print(Option_t * option) const
09fc14a0 767{
b2a60966 768 // Print the parameters used for the particle type identification
f035f6ce 769 cout << "=============== AliPHOSPID1 ================" << endl ;
770 cout << "Making PID "<< endl ;
771 cout << " Headers file: " << fHeaderFileName.Data() << endl ;
772 cout << " RecPoints branch title: " << fRecPointsTitle.Data() << endl ;
7b7c1533 773 cout << " TrackSegments Branch title: " << fTrackSegmentsTitle.Data() << endl ;
774 cout << " RecParticles Branch title " << fRecParticlesTitle.Data() << endl;
f035f6ce 775 cout << "with parameters: " << endl ;
148b2bba 776 // cout << " Maximal EMC - CPV distance (cm) " << fCpvEmcDistance << endl ;
777 //cout << " Time Gate used: " << fTimeGate << endl ;
778 //cout << " Principal Ellipse Parameters " << endl ;
779 //cout << " Ellipse center (x,y) (" << fX_center<<","<<fY_center<<")"<< endl;
780 //cout << " Ellipse focus (a,b) (" << fA<<","<<fB<<")"<< endl;
781 //cout << " Ellipse angle " << fAngle<< endl;
f035f6ce 782 cout << "============================================" << endl ;
09fc14a0 783}
784
7acf6008 785//____________________________________________________________________________
7b7c1533 786void AliPHOSPIDv1::WriteRecParticles(Int_t event)
09fc14a0 787{
7b7c1533 788
789 AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ;
bf8f1fbd 790 TString taskName(GetName()) ;
791 taskName.Remove(taskName.Index(Version())-1) ;
792 TClonesArray * recParticles = gime->RecParticles(taskName) ;
793 recParticles->Expand(recParticles->GetEntriesFast() ) ;
9688c1dd 794
7b7c1533 795 //Make branch in TreeR for RecParticles
7acf6008 796 char * filename = 0;
797 if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){ //generating file name
798 filename = new char[strlen(gAlice->GetBaseFile())+20] ;
799 sprintf(filename,"%s/PHOS.Reco.root",gAlice->GetBaseFile()) ;
800 }
bf8f1fbd 801
7acf6008 802 TDirectory *cwd = gDirectory;
803
804 //First rp
805 Int_t bufferSize = 32000 ;
7b7c1533 806 TBranch * rpBranch = gAlice->TreeR()->Branch("PHOSRP",&recParticles,bufferSize);
bf8f1fbd 807 rpBranch->SetTitle(fRecParticlesTitle);
7acf6008 808 if (filename) {
809 rpBranch->SetFile(filename);
810 TIter next( rpBranch->GetListOfBranches());
761e34c0 811 TBranch * sb ;
812 while ((sb=(TBranch*)next())) {
813 sb->SetFile(filename);
7acf6008 814 }
815 cwd->cd();
816 }
9688c1dd 817
7acf6008 818 //second, pid
819 Int_t splitlevel = 0 ;
820 AliPHOSPIDv1 * pid = this ;
7b7c1533 821 TBranch * pidBranch = gAlice->TreeR()->Branch("AliPHOSPID","AliPHOSPIDv1",&pid,bufferSize,splitlevel);
9688c1dd 822 pidBranch->SetTitle(fRecParticlesTitle.Data());
7acf6008 823 if (filename) {
824 pidBranch->SetFile(filename);
825 TIter next( pidBranch->GetListOfBranches());
761e34c0 826 TBranch * sb ;
827 while ((sb=(TBranch*)next())) {
828 sb->SetFile(filename);
7acf6008 829 }
830 cwd->cd();
831 }
832
761e34c0 833 rpBranch->Fill() ;
834 pidBranch->Fill() ;
9688c1dd 835
7acf6008 836 gAlice->TreeR()->Write(0,kOverwrite) ;
148b2bba 837 //pidBranch->Write(0,kOverwrite) ;
7b7c1533 838
148b2bba 839 delete [] filename ;
7acf6008 840}
69183710 841
7acf6008 842//____________________________________________________________________________
9688c1dd 843TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv)const
7acf6008 844{
845 // Calculates the momentum direction:
846 // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint
9688c1dd 847 // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints
848 // However because of the poor position resolution of PPSD the direction is always taken as if we were
7acf6008 849 // in case 1.
850
851 TVector3 dir(0,0,0) ;
852
853 TVector3 emcglobalpos ;
854 TMatrix dummy ;
855
856 emc->GetGlobalPosition(emcglobalpos, dummy) ;
857
148b2bba 858
7acf6008 859 dir = emcglobalpos ;
860 dir.SetZ( -dir.Z() ) ; // why ?
861 dir.SetMag(1.) ;
862
863 //account correction to the position of IP
864 Float_t xo,yo,zo ; //Coordinates of the origin
865 gAlice->Generator()->GetOrigin(xo,yo,zo) ;
866 TVector3 origin(xo,yo,zo);
867 dir = dir - origin ;
868
869 return dir ;
870}
871//____________________________________________________________________________
a4e98857 872void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
873{
dd5c4038 874 // Print table of reconstructed particles
875
7b7c1533 876 AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ;
bf8f1fbd 877
878 TString taskName(GetName()) ;
879 taskName.Remove(taskName.Index(Version())-1) ;
880 TClonesArray * recParticles = gime->RecParticles(taskName) ;
9688c1dd 881
01a599c9 882 cout << "AliPHOSPIDv1: event "<<gAlice->GetEvNumber() << endl ;
7b7c1533 883 cout << " found " << recParticles->GetEntriesFast() << " RecParticles " << endl ;
9688c1dd 884
7acf6008 885 if(strstr(option,"all")) { // printing found TS
886
887 cout << " PARTICLE "
888 << " Index " << endl ;
7acf6008 889
890 Int_t index ;
7b7c1533 891 for (index = 0 ; index < recParticles->GetEntries() ; index++) {
148b2bba 892 AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ;
893 //cout<<" Type " <<rp->GetType()<<endl;
894 //Text_t particle[11];
895 // switch(rp->GetType()) {
896// case AliPHOSFastRecParticle::kCHARGEDHASLOW:
897// strcpy(particle, "CHARGED HA SLOW") ;
898// break ;
899// case AliPHOSFastRecParticle::kNEUTRALHASLOW:
900// strcpy(particle, "NEUTRAL HA SLOW");
901// break ;
902// case AliPHOSFastRecParticle::kCHARGEDHAFAST:
903// strcpy(particle, "CHARGED HA FAST") ;
904// break ;
905// case AliPHOSFastRecParticle::kNEUTRALHAFAST:
906// strcpy(particle, "NEUTRAL HA FAST");
907// break;
908// case AliPHOSFastRecParticle::kCHARGEDEMSLOW:
909// strcpy(particle, "CHARGED EM SLOW") ;
910// break ;
911// case AliPHOSFastRecParticle::kNEUTRALEMSLOW:
912// strcpy(particle, "NEUTRAL EM SLOW");
913// break ;
914// case AliPHOSFastRecParticle::kCHARGEDEMFAST:
915// strcpy(particle, "CHARGED EM FAST") ;
916// break ;
917// case AliPHOSFastRecParticle::kNEUTRALEMFAST:
918// strcpy( particle, "NEUTRAL EM FAST");
919// break;
920// }
921
922 cout << setw(10) << rp->GetType() << " "
923 << setw(5) << rp->GetIndexInList() << " " <<endl;
924
7acf6008 925 }
926 cout << "-------------------------------------------" << endl ;
927 }
928
69183710 929}
930
7acf6008 931
932