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 }