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 }