Unused DAs removed.
[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     
243     AliPHOSRecPoint    * cpv = 0 ;
244     if(ts->GetCpvIndex()>=0)
245       cpv = (AliPHOSRecPoint *)   fCPVRecPoints->At(ts->GetCpvIndex()) ;
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)
277       if(ts->GetCpvDistance("R") < fCpvEmcDistance) 
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 //____________________________________________________________________________
294 void  AliPHOSPIDv0:: Print(const Option_t *) const
295 {
296   // Print the parameters used for the particle type identification
297   TString message ; 
298   message  = "=============== AliPHOSPIDv0 ================\n" ;
299   message += "Making PID\n" ;
300   message += "with parameters:\n"  ;
301   message += "    Maximal EMC - CPV  distance (cm) %f\n" ;
302   AliInfo(Form( message.Data(),  
303        fCpvEmcDistance ));
304
305   if(fIDOptions.Contains("dis",TString::kIgnoreCase ))
306     AliInfo(Form("                    dispersion cut %f",  fDispersion )) ;
307   if(fIDOptions.Contains("ell",TString::kIgnoreCase ))
308     AliInfo(Form("             Eliptic cuts function: %s",  
309                  fFormula->GetTitle() )) ;
310   if(fIDOptions.Contains("tim",TString::kIgnoreCase ))
311     AliInfo(Form("             Time Gate used: %f",  fTimeGate)) ;
312 }
313
314 //____________________________________________________________________________
315 void  AliPHOSPIDv0::SetShowerProfileCut(const char * formula)
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 }
323
324 //____________________________________________________________________________
325 void  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 //____________________________________________________________________________
358 TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * )const 
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
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.
397   //account correction to the position of IP
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 ;
402
403   return dir ;  
404 }
405 //____________________________________________________________________________
406 void AliPHOSPIDv0::PrintRecParticles(Option_t * option)
407 {
408   // Print table of reconstructed particles
409
410   TString message ; 
411   message = "Found %d RecParticles\n" ; 
412   AliInfo(Form(message.Data(), 
413                fRecParticles->GetEntriesFast() )) ;   
414
415   if(strstr(option,"all")) {  // printing found TS
416     AliInfo("  PARTICLE   Index"  ) ; 
417    
418     Int_t index ;
419     for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
420       AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;       
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       
454       AliInfo(Form("          %s     %d",  
455                    particle, rp->GetIndexInList())) ;
456     }
457   }  
458 }