Charged jets (pPb): Improved trackcut analysis
[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   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   }
152   return *this;
153 }
154
155 //____________________________________________________________________________
156 AliPHOSPIDv0::~AliPHOSPIDv0()
157 {
158   //Empty dtor, fFormula leaks 
159 }
160
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 //}
189
190 //____________________________________________________________________________
191 void  AliPHOSPIDv0::TrackSegments2RecParticles(Option_t * option) 
192 {
193   //Steering method
194
195   if(strstr(option,"tim"))
196     gBenchmark->Start("PHOSPID");
197   
198   if(strstr(option,"print")) {
199     Print() ; 
200     return ; 
201   }
202
203   AliInfo(Form("%d emc clusters, %d track segments", 
204                fEMCRecPoints->GetEntriesFast(), 
205                fTrackSegments->GetEntriesFast())) ;
206
207   MakeRecParticles() ;
208     
209   if(strstr(option,"deb"))
210     PrintRecParticles(option) ;
211
212   if(strstr(option,"tim")){
213     gBenchmark->Stop("PHOSPID");
214     AliInfo(Form("took %f seconds for PID", 
215                  gBenchmark->GetCpuTime("PHOSPID"))); 
216   }
217   
218 }
219
220 //____________________________________________________________________________
221 void  AliPHOSPIDv0::MakeRecParticles()
222 {
223   // Reconstructs the particles from the tracksegments
224
225   fRecParticles->Clear();
226   
227   TIter next(fTrackSegments) ; 
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     
238     new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
239     rp = (AliPHOSRecParticle *)fRecParticles->At(index) ; 
240     rp->SetTrackSegment(index) ;
241     rp->SetIndexInList(index) ;
242     
243     AliPHOSEmcRecPoint * emc = 0 ;
244     if(ts->GetEmcIndex()>=0)
245       emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
246     else
247       AliFatal(Form("ts->GetEmcIndex() return illegal index %d",ts->GetEmcIndex()));
248     
249     AliPHOSRecPoint    * cpv = 0 ;
250     if(ts->GetCpvIndex()>=0)
251       cpv = (AliPHOSRecPoint *)   fCPVRecPoints->At(ts->GetCpvIndex()) ;
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)
283       if(ts->GetCpvDistance("R") < fCpvEmcDistance) 
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 //____________________________________________________________________________
300 void  AliPHOSPIDv0:: Print(const Option_t *) const
301 {
302   // Print the parameters used for the particle type identification
303   TString message ; 
304   message  = "=============== AliPHOSPIDv0 ================\n" ;
305   message += "Making PID\n" ;
306   message += "with parameters:\n"  ;
307   message += "    Maximal EMC - CPV  distance (cm) %f\n" ;
308   AliInfo(Form( message.Data(),  
309        fCpvEmcDistance ));
310
311   if(fIDOptions.Contains("dis",TString::kIgnoreCase ))
312     AliInfo(Form("                    dispersion cut %f",  fDispersion )) ;
313   if(fIDOptions.Contains("ell",TString::kIgnoreCase ))
314     AliInfo(Form("             Eliptic cuts function: %s",  
315                  fFormula->GetTitle() )) ;
316   if(fIDOptions.Contains("tim",TString::kIgnoreCase ))
317     AliInfo(Form("             Time Gate used: %f",  fTimeGate)) ;
318 }
319
320 //____________________________________________________________________________
321 void  AliPHOSPIDv0::SetShowerProfileCut(const char * formula)
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 }
329
330 //____________________________________________________________________________
331 void  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 //____________________________________________________________________________
364 TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * )const 
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
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.
403   //account correction to the position of IP
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 ;
408
409   return dir ;  
410 }
411 //____________________________________________________________________________
412 void AliPHOSPIDv0::PrintRecParticles(Option_t * option)
413 {
414   // Print table of reconstructed particles
415
416   TString message ; 
417   message = "Found %d RecParticles\n" ; 
418   AliInfo(Form(message.Data(), 
419                fRecParticles->GetEntriesFast() )) ;   
420
421   if(strstr(option,"all")) {  // printing found TS
422     AliInfo("  PARTICLE   Index"  ) ; 
423    
424     Int_t index ;
425     for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
426       AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;       
427       
428       Text_t particle[100];
429       switch(rp->GetType()) {
430       case  AliPHOSFastRecParticle::kNEUTRALEMFAST:
431         strncpy( particle, "NEUTRAL EM FAST",100);
432         break;
433       case  AliPHOSFastRecParticle::kNEUTRALHAFAST:
434         strncpy(particle, "NEUTRAL HA FAST",100);
435         break;
436       case  AliPHOSFastRecParticle::kNEUTRALEMSLOW:
437         strncpy(particle, "NEUTRAL EM SLOW",100);
438         break ;
439       case  AliPHOSFastRecParticle::kNEUTRALHASLOW: 
440         strncpy(particle, "NEUTRAL HA SLOW",100);
441         break ;
442       case  AliPHOSFastRecParticle::kCHARGEDEMFAST:
443         strncpy(particle, "CHARGED EM FAST",100);
444         break ;
445       case  AliPHOSFastRecParticle::kCHARGEDHAFAST:
446         strncpy(particle, "CHARGED HA FAST",100);
447         break ; 
448       case  AliPHOSFastRecParticle::kCHARGEDEMSLOW:
449         strncpy(particle, "CHARGEDEMSLOW",100);
450         break ;
451       case  AliPHOSFastRecParticle::kCHARGEDHASLOW:
452         strncpy(particle, "CHARGED HA SLOW",100);
453         break ; 
454       }
455       
456       //    Int_t * primaries; 
457       //    Int_t nprimaries;
458       //    primaries = rp->GetPrimaries(nprimaries);
459       
460       AliInfo(Form("          %s     %d",  
461                    particle, rp->GetIndexInList())) ;
462     }
463   }  
464 }