Implement the new logging system (AliLog)
[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 //_________________________________________________________________________
19 // Implementation version v0 of the PHOS particle identifier 
20 // Particle identification based on the 
21 //     - CPV information, 
22 //     - Preshower information (in MIXT or GPS2 geometries)
23 //     - shower width.
24 //
25 // CPV or Preshower clusters should be closer in PHOS plane than fCpvEmcDistance (in cm).
26 // This parameter can be set by method SetCpvtoEmcDistanceCut(Float_t cut)  
27 //
28 // One can set desirable ID method by the function SetIdentificationMethod(option).
29 // Presently the following options can be used together or separately :
30 //     - "disp": use dispersion cut on shower width 
31 //               (width can be set by method SetDispersionCut(Float_t cut)
32 //     - "ell" : use cut on the axis of the ellipse, drawn around shower 
33 //       (this cut can be changed by SetShowerProfileCut(char* formula), 
34 //        where formula - any function of two variables f(lambda[0],lambda[1]).
35 //        Shower is considered as EM if f() > 0 )
36 // One can visualize current cuts calling method PlotDispersionCuts().    
37 //
38 // use case:
39 //  root [0] AliPHOSPIDv0 * p1 = new AliPHOSPIDv0("galice.root")
40 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
41 //  root [1] p1->SetIdentificationMethod("disp ellipse")
42 //  root [2] p1->ExecuteTask()
43 //  root [3] AliPHOSPIDv0 * p2 = new AliPHOSPIDv0("galice1.root","ts1")
44 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
45 //                // reading headers from file galice1.root and TrackSegments 
46 //                // with title "ts1"
47 //  root [4] p2->SetRecParticlesBranch("rp1")
48 //                // set file name for the branch RecParticles
49 //  root [5] p2->ExecuteTask("deb all time")
50 //                // available options
51 //                // "deb" - prints # of reconstructed particles
52 //                // "deb all" -  prints # and list of RecParticles
53 //                // "time" - prints benchmarking results
54 //                  
55 //*-- Author: Yves Schutz (SUBATECH)  & Gines Martinez (SUBATECH) & 
56 //            Dmitri Peressounko (SUBATECH & Kurchatov Institute)
57 //            Completely redesined by Dmitri Peressounko, March 2001
58
59 // --- ROOT system ---
60 #include "TTree.h"
61 #include "TF2.h"
62 #include "TFormula.h"
63 #include "TCanvas.h"
64
65 #include "TBenchmark.h"
66 // --- Standard library ---
67
68 // --- AliRoot header files ---
69 #include "AliLog.h"
70 #include "AliRun.h"
71 #include "AliGenerator.h"
72 #include "AliPHOSPIDv0.h"
73 #include "AliPHOSEmcRecPoint.h"
74 #include "AliPHOSTrackSegment.h"
75 #include "AliPHOSRecParticle.h"
76 #include "AliPHOSGeometry.h"
77 #include "AliPHOSLoader.h"
78
79 ClassImp( AliPHOSPIDv0) 
80
81 //____________________________________________________________________________
82 AliPHOSPIDv0::AliPHOSPIDv0():AliPHOSPID()
83
84   // default ctor
85   fFormula           = 0 ;
86   fDispersion        = 0. ; 
87   fCpvEmcDistance    = 0 ; 
88   fTimeGate          = 2.e-9 ;
89   fEventFolderName    = "" ; 
90   fTrackSegmentsTitle= "" ; 
91   fRecPointsTitle    = "" ; 
92   fRecParticlesTitle = "" ; 
93   fIDOptions         = "dis time" ; 
94   fRecParticlesInRun = 0 ;
95   fClusterizer = 0;
96   fTSMaker = 0;
97 }
98
99 //____________________________________________________________________________
100 AliPHOSPIDv0::AliPHOSPIDv0(const char * evFolderName,const char * name) : AliPHOSPID(evFolderName, name)
101
102   //ctor with the indication on where to look for the track segments
103
104   fFormula        = new TFormula("LambdaCuts","(x>1)*(x<2.5)*(y>0)*(y<x)") ;   
105   fDispersion     = 2.0 ; 
106   fCpvEmcDistance = 3.0 ;
107   fTimeGate          = 2.e-9 ;
108  
109   fEventFolderName     = GetTitle() ; 
110   fTrackSegmentsTitle = GetName();
111   fRecPointsTitle     = GetName();
112   fRecParticlesTitle  = GetName();
113   fIDOptions          = "dis time" ;
114     
115   fRecParticlesInRun = 0 ; 
116
117   Init() ;
118
119 }
120
121 //____________________________________________________________________________
122 AliPHOSPIDv0::~AliPHOSPIDv0()
123
124 }
125
126 //____________________________________________________________________________
127 Float_t  AliPHOSPIDv0::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t *  Axis)const
128 {
129   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
130  
131   const AliPHOSGeometry * geom = AliPHOSLoader::GetPHOSGeometry() ; 
132   TVector3 vecEmc ;
133   TVector3 vecCpv ;
134   
135   emc->GetLocalPosition(vecEmc) ;
136   cpv->GetLocalPosition(vecCpv) ; 
137   if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ 
138     
139     // Correct to difference in CPV and EMC position due to different distance to center.
140     // we assume, that particle moves from center
141     Float_t dCPV = geom->GetIPtoOuterCoverDistance();
142     Float_t dEMC = geom->GetIPtoCrystalSurface() ;
143     dEMC         = dEMC / dCPV ;
144     vecCpv = dEMC * vecCpv  - vecEmc ; 
145     if (Axis == "X") return vecCpv.X();
146     if (Axis == "Y") return vecCpv.Y();
147     if (Axis == "Z") return vecCpv.Z();
148     if (Axis == "R") return vecCpv.Mag();
149   } 
150  
151   return 100000000 ;
152 }
153
154 //____________________________________________________________________________
155 void  AliPHOSPIDv0::Exec(Option_t * option) 
156 {
157   //Steering method
158   
159   if( strcmp(GetName(), "")== 0 ) 
160     Init() ;
161   
162   if(strstr(option,"tim"))
163     gBenchmark->Start("PHOSPID");
164   
165   if(strstr(option,"print")) {
166     Print() ; 
167     return ; 
168   }
169
170   AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
171   if(runget == 0x0) 
172    {
173      AliError(Form("Can not find run getter in event folder \"%s\"",
174                    GetTitle()));
175      return;
176    }
177   
178   AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
179   if ( gime == 0 ) 
180    {
181      AliError("Could not obtain the Loader object !"); 
182      return ;
183    } 
184
185   if(gime->BranchExists("RecParticles") )
186     return ;
187
188   
189   Int_t nevents = runget->GetNumberOfEvents() ;       //(Int_t) gAlice->TreeE()->GetEntries() ;
190
191   Int_t ievent ;
192   
193   for(ievent = 0; ievent < nevents; ievent++){
194     runget->GetEvent(ievent);
195     AliInfo(Form("event %d %d %d", 
196                  ievent, gime->EmcRecPoints(), 
197                  gime->TrackSegments())) ;
198     MakeRecParticles() ;
199     
200     WriteRecParticles();
201     
202     if(strstr(option,"deb"))
203       PrintRecParticles(option) ;
204
205     //increment the total number of rec particles per run 
206     fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; 
207
208   }
209   
210   if(strstr(option,"tim")){
211     gBenchmark->Stop("PHOSPID");
212     AliInfo(Form("took %f seconds for PID %f seconds per event", 
213                  gBenchmark->GetCpuTime("PHOSPID"), 
214                  gBenchmark->GetCpuTime("PHOSPID")/nevents)) ; 
215   }
216   
217 }
218 //____________________________________________________________________________
219 void AliPHOSPIDv0::Init()
220 {
221   // Make all memory allocations that are not possible in default constructor
222   // Add the PID task to the list of PHOS tasks
223   
224   AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
225   if(runget == 0x0) 
226    {
227      AliError(Form("Can not find run getter in event folder \"%s\"",
228                    GetTitle()));
229      return;
230    }
231   
232   AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
233   if ( gime == 0 ) 
234    {
235      AliError("Could not obtain the Loader object !"); 
236      return ;
237    } 
238    
239   gime->PostPID(this);
240   gime->LoadRecParticles("UPDATE");
241   
242 }
243
244 //____________________________________________________________________________
245 void  AliPHOSPIDv0::MakeRecParticles()
246 {
247   // Reconstructs the particles from the tracksegments
248
249   TString taskName(GetName()) ; 
250   taskName.Remove(taskName.Index(Version())-1) ;
251
252   AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
253   if(runget == 0x0) 
254    {
255      AliError(Form("Can not find run getter in event folder \"%s\"",
256                    GetTitle()));
257      return;
258    }
259   
260   AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
261   if ( gime == 0 ) 
262    {
263      AliError("Could not obtain the Loader object !"); 
264      return ;
265    } 
266
267   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
268   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
269   TClonesArray * trackSegments = gime->TrackSegments() ; 
270   TClonesArray * recParticles  = gime->RecParticles() ; 
271
272   recParticles->Clear();
273   
274   TIter next(trackSegments) ; 
275   AliPHOSTrackSegment * ts ; 
276   Int_t index = 0 ; 
277   AliPHOSRecParticle * rp ; 
278   
279   Bool_t ellips = fIDOptions.Contains("ell",TString::kIgnoreCase ) ;
280   Bool_t disp   = fIDOptions.Contains("dis",TString::kIgnoreCase ) ;
281   Bool_t time   = fIDOptions.Contains("tim",TString::kIgnoreCase ) ;
282   
283   while ( (ts = (AliPHOSTrackSegment *)next()) ) {
284     
285     new( (*recParticles)[index] ) AliPHOSRecParticle() ;
286     rp = (AliPHOSRecParticle *)recParticles->At(index) ; 
287     rp->SetTrackSegment(index) ;
288     rp->SetIndexInList(index) ;
289     
290     AliPHOSEmcRecPoint * emc = 0 ;
291     if(ts->GetEmcIndex()>=0)
292       emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
293     
294     AliPHOSRecPoint    * cpv = 0 ;
295     if(ts->GetCpvIndex()>=0)
296       cpv = (AliPHOSRecPoint *)   cpvRecPoints->At(ts->GetCpvIndex()) ;
297     
298     //set momentum and energy first
299     Float_t    e = emc->GetEnergy() ;
300     TVector3 dir = GetMomentumDirection(emc,cpv) ; 
301     dir.SetMag(e) ;
302
303     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
304     rp->SetCalcMass(0);
305
306     //now set type (reconstructed) of the particle    
307     Int_t showerprofile = 0;  // 0 narrow and 1 wide
308     
309     if(ellips){
310       Float_t lambda[2] ;
311       emc->GetElipsAxis(lambda) ;
312       if(fFormula->Eval(lambda[0],lambda[1]) <= 0 )
313         showerprofile = 1 ;  // not narrow
314     }
315     
316     if(disp)
317       if(emc->GetDispersion() > fDispersion )
318         showerprofile = 1 ;  // not narrow
319     
320     Int_t slow = 0 ;
321     if(time)
322       if(emc->GetTime() > fTimeGate )
323         slow = 0 ; 
324         
325     // Looking at the CPV detector
326     Int_t cpvdetector= 0 ;  //1 hit and 0 no hit     
327     if(cpv)
328       if(GetDistance(emc, cpv,  "R") < fCpvEmcDistance) 
329         cpvdetector = 1 ;  
330     
331     Int_t type = showerprofile + 2 * slow  + 4 * cpvdetector ;
332     rp->SetType(type) ; 
333     rp->SetProductionVertex(0,0,0,0);
334     rp->SetFirstMother(-1);
335     rp->SetLastMother(-1);
336     rp->SetFirstDaughter(-1);
337     rp->SetLastDaughter(-1);
338     rp->SetPolarisation(0,0,0);
339     index++ ; 
340   }
341   
342 }
343
344 //____________________________________________________________________________
345 void  AliPHOSPIDv0:: Print() const
346 {
347   // Print the parameters used for the particle type identification
348   TString message ; 
349   message  = "=============== AliPHOSPIDv0 ================\n" ;
350   message += "Making PID\n" ;
351   message += "    Headers file:               %s\n" ; 
352   message += "    RecPoints branch title:     %s\n" ;
353   message += "    TrackSegments Branch title: %s\n" ; 
354   message += "    RecParticles Branch title   %s\n" ;  
355   message += "with parameters:\n"  ;
356   message += "    Maximal EMC - CPV  distance (cm) %f\n" ;
357   AliInfo(Form( message.Data(),  
358        GetTitle(), 
359        fRecPointsTitle.Data(), 
360        fTrackSegmentsTitle.Data(), 
361        fRecParticlesTitle.Data(), 
362        fCpvEmcDistance ));
363
364   if(fIDOptions.Contains("dis",TString::kIgnoreCase ))
365     AliInfo(Form("                    dispersion cut %f",  fDispersion )) ;
366   if(fIDOptions.Contains("ell",TString::kIgnoreCase ))
367     AliInfo(Form("             Eliptic cuts function: %s",  
368                  fFormula->GetTitle() )) ;
369   if(fIDOptions.Contains("tim",TString::kIgnoreCase ))
370     AliInfo(Form("             Time Gate used: %f",  fTimeGate)) ;
371 }
372
373 //____________________________________________________________________________
374 void  AliPHOSPIDv0::SetShowerProfileCut(const char * formula)
375 {
376   //set shape of the cut on the axis of ellipce, drown around shouer
377   //shower considered "narrow" if Formula(lambda[0],lambda[1]) > 0.
378   if(fFormula) 
379     delete fFormula; 
380   fFormula = new TFormula("Lambda Cut",formula) ;
381 }
382 //____________________________________________________________________________
383 void  AliPHOSPIDv0::WriteRecParticles()
384 {
385   // Saves the reconstructed particles too a file
386  
387   AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
388   if(runget == 0x0) 
389    {
390      AliError(Form("Can not find run getter in event folder \"%s\"",
391                    GetTitle()));
392      return;
393    }
394   
395   AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
396   if ( gime == 0 ) 
397    {
398      AliError("Could not obtain the Loader object !"); 
399      return ;
400    } 
401
402   TClonesArray * recParticles = gime->RecParticles() ; 
403   recParticles->Expand(recParticles->GetEntriesFast() ) ;
404
405   TTree * treeR = gime->TreeR();
406   
407   if(!treeR){
408     gime->MakeTree("R");
409     treeR = gime->TreeR() ;
410   }
411   
412   //First rp
413   Int_t bufferSize = 32000 ;    
414   TBranch * rpBranch = treeR->Branch("PHOSRP",&recParticles,bufferSize);
415   rpBranch->SetTitle(fRecParticlesTitle);
416   
417   //second, pid
418   Int_t splitlevel = 0 ; 
419   AliPHOSPIDv0 * pid = this ;
420   TBranch * pidBranch = treeR->Branch("AliPHOSPID","AliPHOSPIDv0",&pid,bufferSize,splitlevel);
421   pidBranch->SetTitle(fRecParticlesTitle.Data());
422   
423   rpBranch->Fill() ;
424   pidBranch->Fill() ;
425
426   gime->WriteRecParticles("OVERWRITE");
427   gime->WritePID("OVERWRITE");
428
429 }
430
431 //____________________________________________________________________________
432 void  AliPHOSPIDv0::PlotDispersionCuts()const
433 {
434   // produces a plot of the dispersion cut
435   TCanvas*  lambdas = new TCanvas("lambdas","Cuts on the ellipse axis",200,10,700,500);
436  
437   if(fIDOptions.Contains("ell",TString::kIgnoreCase ) ){
438     TF2 * ell = new TF2("Elliptic Cuts",fFormula->GetName(),0,3,0,3) ;
439     ell->SetMinimum(0.0000001) ;
440     ell->SetMaximum(0.001) ;
441     ell->SetLineStyle(1) ;
442     ell->SetLineWidth(2) ;
443     ell->Draw() ;
444   }
445   
446   if( fIDOptions.Contains("dis",TString::kIgnoreCase ) ){
447     TF2 * dsp = new TF2("dispersion","(y<x)*(x*x+y*y < [0]*[0])",0,3,0,3) ;
448     dsp->SetParameter(0,fDispersion) ;
449     dsp->SetMinimum(0.0000001) ;
450     dsp->SetMaximum(0.001) ;
451     dsp->SetLineStyle(1) ;
452     dsp->SetLineColor(2) ;
453     dsp->SetLineWidth(2) ;
454     dsp->SetNpx(200) ;
455     dsp->SetNpy(200) ;
456     if(fIDOptions.Contains("ell",TString::kIgnoreCase ) )
457       dsp->Draw("same") ;
458     else
459       dsp->Draw() ;
460   }
461   lambdas->Update();
462 }
463
464 //____________________________________________________________________________
465 TVector3 AliPHOSPIDv0::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * )const 
466
467   // Calculates the momentum direction:
468   //   1. if only a EMC RecPoint, direction is given by IP and this RecPoint
469   //   2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints 
470   //  However because of the poor position resolution of PPSD the direction is always taken as if we were 
471   //  in case 1.
472
473   TVector3 dir(0,0,0) ; 
474   
475   TVector3 emcglobalpos ;
476   TMatrix  dummy ;
477   
478   emc->GetGlobalPosition(emcglobalpos, dummy) ;
479   
480  
481   // The following commented code becomes valid once the PPSD provides 
482   // a reasonable position resolution, at least as good as EMC ! 
483   //   TVector3 ppsdlglobalpos ;
484   //   TVector3 ppsduglobalpos ;
485   //   if( fPpsdLowRecPoint ){ // certainly a photon that has concerted
486   //     fPpsdLowRecPoint->GetGlobalPosition(ppsdlglobalpos, mdummy) ; 
487   //     dir = emcglobalpos -  ppsdlglobalpos ; 
488   //     if( fPpsdUpRecPoint ){ // not looks like a charged       
489   //        fPpsdUpRecPoint->GetGlobalPosition(ppsduglobalpos, mdummy) ; 
490   //        dir = ( dir +  emcglobalpos -  ppsduglobalpos ) * 0.5 ; 
491   //      }
492   //   }
493   //   else { // looks like a neutral
494   //    dir = emcglobalpos ;  
495   //  }
496   
497   dir = emcglobalpos ;  
498   dir.SetZ( -dir.Z() ) ;   // why ?  
499   dir.SetMag(1.) ;
500
501   //account correction to the position of IP
502   Float_t xo,yo,zo ; //Coordinates of the origin
503   gAlice->Generator()->GetOrigin(xo,yo,zo) ;
504   TVector3 origin(xo,yo,zo);
505   dir = dir - origin ;
506
507   return dir ;  
508 }
509 //____________________________________________________________________________
510 void AliPHOSPIDv0::PrintRecParticles(Option_t * option)
511 {
512   // Print table of reconstructed particles
513
514   AliRunLoader* runget = AliRunLoader::GetRunLoader(GetTitle());
515   if(runget == 0x0) 
516    {
517      AliError(Form("Can not find run getter in event folder \"%s\"",
518                    GetTitle()));
519      return;
520    }
521   
522   AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(runget->GetLoader("PHOSLoader"));
523   if ( gime == 0 ) 
524    {
525      AliError("Could not obtain the Loader object !"); 
526      return ;
527    } 
528
529   TString taskName(GetName()) ; 
530   taskName.Remove(taskName.Index(Version())-1) ;
531   TClonesArray * recParticles = gime->RecParticles() ; 
532   
533   TString message ; 
534   message  = "event %d\n" ; 
535   message += "       found %d RecParticles\n" ; 
536   AliInfo(Form(message.Data(), 
537                gAlice->GetEvNumber(), recParticles->GetEntriesFast() )) ;   
538
539   if(strstr(option,"all")) {  // printing found TS
540     AliInfo("  PARTICLE   Index"  ) ; 
541    
542     Int_t index ;
543     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
544       AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ;       
545       
546       Text_t particle[11];
547       switch(rp->GetType()) {
548       case  AliPHOSFastRecParticle::kNEUTRALEMFAST:
549         strcpy( particle, "NEUTRAL EM FAST");
550         break;
551       case  AliPHOSFastRecParticle::kNEUTRALHAFAST:
552         strcpy(particle, "NEUTRAL HA FAST");
553         break;
554       case  AliPHOSFastRecParticle::kNEUTRALEMSLOW:
555         strcpy(particle, "NEUTRAL EM SLOW");
556         break ;
557       case  AliPHOSFastRecParticle::kNEUTRALHASLOW: 
558         strcpy(particle, "NEUTRAL HA SLOW");
559         break ;
560       case  AliPHOSFastRecParticle::kCHARGEDEMFAST:
561         strcpy(particle, "CHARGED EM FAST") ;
562         break ;
563       case  AliPHOSFastRecParticle::kCHARGEDHAFAST:
564         strcpy(particle, "CHARGED HA FAST") ;
565         break ; 
566       case  AliPHOSFastRecParticle::kCHARGEDEMSLOW:
567         strcpy(particle, "CHARGEDEMSLOW") ;
568         break ;
569       case  AliPHOSFastRecParticle::kCHARGEDHASLOW:
570         strcpy(particle, "CHARGED HA SLOW") ;
571         break ; 
572       }
573       
574       //    Int_t * primaries; 
575       //    Int_t nprimaries;
576       //    primaries = rp->GetPrimaries(nprimaries);
577       
578       AliInfo(Form("          %s     %d",  
579                    particle, rp->GetIndexInList())) ;
580     }
581   }  
582 }
583
584
585