L1phase shift corrected
[u/mrichter/AliRoot.git] / PHOS / AliPHOSPIDv0.cxx
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
18 /* History of cvs commits:
19  *
20  * $Log$
21  * Revision 1.15  2007/03/06 06:57:05  kharlov
22  * DP:calculation of distance to CPV done in TSM
23  *
24  * Revision 1.14  2006/09/07 18:31:08  kharlov
25  * Effective c++ corrections (T.Pocheptsov)
26  *
27  * Revision 1.13  2005/05/28 14:19:04  schutz
28  * Compilation warnings fixed by T.P.
29  *
30  */
31
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 ---
74 #include "TTree.h"
75 #include "TF2.h"
76 #include "TFormula.h"
77 #include "TCanvas.h"
78 #include "TClonesArray.h"
79
80 #include "TBenchmark.h"
81 // --- Standard library ---
82
83 // --- AliRoot header files ---
84 #include "AliLog.h"
85 #include "AliPHOSPIDv0.h"
86 #include "AliPHOSEmcRecPoint.h"
87 #include "AliPHOSTrackSegment.h"
88 #include "AliPHOSRecParticle.h"
89 #include "AliPHOSGeometry.h"
90
91 ClassImp( AliPHOSPIDv0) 
92
93 //____________________________________________________________________________
94 AliPHOSPIDv0::AliPHOSPIDv0():
95   fIDOptions("dis time"),
96   fClusterizer(0),
97   fTSMaker(0),
98   fFormula(0),
99   fDispersion(0.f),
100   fCpvEmcDistance(0.f),
101   fTimeGate(2.e-9f)
102
103   // default ctor
104 }
105
106 //____________________________________________________________________________
107 AliPHOSPIDv0::AliPHOSPIDv0(AliPHOSGeometry *geom) : 
108   AliPHOSPID(geom),
109   fIDOptions("dis time"),
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),
115   fTimeGate(2.e-9f)
116
117   //ctor
118 }
119
120 //____________________________________________________________________________
121 AliPHOSPIDv0::AliPHOSPIDv0(const AliPHOSPIDv0 & rhs) :
122   AliPHOSPID(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)
130 {
131   //Copy ctor, the same as compiler-generated, possibly wrong if
132   //someone implements dtor correctly.
133 }
134   
135 //____________________________________________________________________________
136 AliPHOSPIDv0 & AliPHOSPIDv0::operator = (const AliPHOSPIDv0 & rhs)
137 {
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;
147
148   return *this;
149 }
150
151 //____________________________________________________________________________
152 AliPHOSPIDv0::~AliPHOSPIDv0()
153 {
154   //Empty dtor, fFormula leaks 
155 }
156
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 //}
185
186 //____________________________________________________________________________
187 void  AliPHOSPIDv0::TrackSegments2RecParticles(Option_t * option) 
188 {
189   //Steering method
190
191   if(strstr(option,"tim"))
192     gBenchmark->Start("PHOSPID");
193   
194   if(strstr(option,"print")) {
195     Print() ; 
196     return ; 
197   }
198
199   AliInfo(Form("%d emc clusters, %d track segments", 
200                fEMCRecPoints->GetEntriesFast(), 
201                fTrackSegments->GetEntriesFast())) ;
202
203   MakeRecParticles() ;
204     
205   if(strstr(option,"deb"))
206     PrintRecParticles(option) ;
207
208   if(strstr(option,"tim")){
209     gBenchmark->Stop("PHOSPID");
210     AliInfo(Form("took %f seconds for PID", 
211                  gBenchmark->GetCpuTime("PHOSPID"))); 
212   }
213   
214 }
215
216 //____________________________________________________________________________
217 void  AliPHOSPIDv0::MakeRecParticles()
218 {
219   // Reconstructs the particles from the tracksegments
220
221   fRecParticles->Clear();
222   
223   TIter next(fTrackSegments) ; 
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     
234     new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
235     rp = (AliPHOSRecParticle *)fRecParticles->At(index) ; 
236     rp->SetTrackSegment(index) ;
237     rp->SetIndexInList(index) ;
238     
239     AliPHOSEmcRecPoint * emc = 0 ;
240     if(ts->GetEmcIndex()>=0)
241       emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
242     else
243       AliFatal(Form("ts->GetEmcIndex() return illegal index %d",ts->GetEmcIndex()));
244     
245     AliPHOSRecPoint    * cpv = 0 ;
246     if(ts->GetCpvIndex()>=0)
247       cpv = (AliPHOSRecPoint *)   fCPVRecPoints->At(ts->GetCpvIndex()) ;
248     
249     //set momentum and energy first
250     Float_t    e = emc->GetEnergy() ;
251     TVector3 dir = GetMomentumDirection(emc,cpv) ; 
252     dir.SetMag(e) ;
253
254     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
255     rp->SetCalcMass(0);
256
257     //now set type (reconstructed) of the particle    
258     Int_t showerprofile = 0;  // 0 narrow and 1 wide
259     
260     if(ellips){
261       Float_t lambda[2] ;
262       emc->GetElipsAxis(lambda) ;
263       if(fFormula->Eval(lambda[0],lambda[1]) <= 0 )
264         showerprofile = 1 ;  // not narrow
265     }
266     
267     if(disp)
268       if(emc->GetDispersion() > fDispersion )
269         showerprofile = 1 ;  // not narrow
270     
271     Int_t slow = 0 ;
272     if(time)
273       if(emc->GetTime() > fTimeGate )
274         slow = 0 ; 
275         
276     // Looking at the CPV detector
277     Int_t cpvdetector= 0 ;  //1 hit and 0 no hit     
278     if(cpv)
279       if(ts->GetCpvDistance("R") < fCpvEmcDistance) 
280         cpvdetector = 1 ;  
281     
282     Int_t type = showerprofile + 2 * slow  + 4 * cpvdetector ;
283     rp->SetType(type) ; 
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);
290     index++ ; 
291   }
292   
293 }
294
295 //____________________________________________________________________________
296 void  AliPHOSPIDv0:: Print(const Option_t *) const
297 {
298   // Print the parameters used for the particle type identification
299   TString message ; 
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(),  
305        fCpvEmcDistance ));
306
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)) ;
314 }
315
316 //____________________________________________________________________________
317 void  AliPHOSPIDv0::SetShowerProfileCut(const char * formula)
318 {
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.
321   if(fFormula) 
322     delete fFormula; 
323   fFormula = new TFormula("Lambda Cut",formula) ;
324 }
325
326 //____________________________________________________________________________
327 void  AliPHOSPIDv0::PlotDispersionCuts()const
328 {
329   // produces a plot of the dispersion cut
330   TCanvas*  lambdas = new TCanvas("lambdas","Cuts on the ellipse axis",200,10,700,500);
331  
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) ;
338     ell->Draw() ;
339   }
340   
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) ;
349     dsp->SetNpx(200) ;
350     dsp->SetNpy(200) ;
351     if(fIDOptions.Contains("ell",TString::kIgnoreCase ) )
352       dsp->Draw("same") ;
353     else
354       dsp->Draw() ;
355   }
356   lambdas->Update();
357 }
358
359 //____________________________________________________________________________
360 TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * )const 
361
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 
366   //  in case 1.
367
368   TVector3 dir(0,0,0) ; 
369   
370   TVector3 emcglobalpos ;
371   TMatrix  dummy ;
372   
373   emc->GetGlobalPosition(emcglobalpos, dummy) ;
374   
375  
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 ; 
386   //      }
387   //   }
388   //   else { // looks like a neutral
389   //    dir = emcglobalpos ;  
390   //  }
391   
392   dir = emcglobalpos ;  
393   dir.SetZ( -dir.Z() ) ;   // why ?  
394   dir.SetMag(1.) ;
395
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 ;
404
405   return dir ;  
406 }
407 //____________________________________________________________________________
408 void AliPHOSPIDv0::PrintRecParticles(Option_t * option)
409 {
410   // Print table of reconstructed particles
411
412   TString message ; 
413   message = "Found %d RecParticles\n" ; 
414   AliInfo(Form(message.Data(), 
415                fRecParticles->GetEntriesFast() )) ;   
416
417   if(strstr(option,"all")) {  // printing found TS
418     AliInfo("  PARTICLE   Index"  ) ; 
419    
420     Int_t index ;
421     for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
422       AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;       
423       
424       Text_t particle[100];
425       switch(rp->GetType()) {
426       case  AliPHOSFastRecParticle::kNEUTRALEMFAST:
427         strncpy( particle, "NEUTRAL EM FAST",100);
428         break;
429       case  AliPHOSFastRecParticle::kNEUTRALHAFAST:
430         strncpy(particle, "NEUTRAL HA FAST",100);
431         break;
432       case  AliPHOSFastRecParticle::kNEUTRALEMSLOW:
433         strncpy(particle, "NEUTRAL EM SLOW",100);
434         break ;
435       case  AliPHOSFastRecParticle::kNEUTRALHASLOW: 
436         strncpy(particle, "NEUTRAL HA SLOW",100);
437         break ;
438       case  AliPHOSFastRecParticle::kCHARGEDEMFAST:
439         strncpy(particle, "CHARGED EM FAST",100);
440         break ;
441       case  AliPHOSFastRecParticle::kCHARGEDHAFAST:
442         strncpy(particle, "CHARGED HA FAST",100);
443         break ; 
444       case  AliPHOSFastRecParticle::kCHARGEDEMSLOW:
445         strncpy(particle, "CHARGEDEMSLOW",100);
446         break ;
447       case  AliPHOSFastRecParticle::kCHARGEDHASLOW:
448         strncpy(particle, "CHARGED HA SLOW",100);
449         break ; 
450       }
451       
452       //    Int_t * primaries; 
453       //    Int_t nprimaries;
454       //    primaries = rp->GetPrimaries(nprimaries);
455       
456       AliInfo(Form("          %s     %d",  
457                    particle, rp->GetIndexInList())) ;
458     }
459   }  
460 }