PWGGA/PHOSTasks/ClusterSelection/AliPHOSClusterSelection.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 #include "AliAnalysisManager.h"
  19 #include "AliInputEventHandler.h"
  20 #include "AliMultiInputEventHandler.h"
  21 #include "AliVEvent.h"
  22 #include "AliESDEvent.h"
  23 #include "AliAODEvent.h"
  24 #include "AliESDCaloCluster.h"
  25 #include "AliVParticle.h"
  26 #include "AliVTrack.h"
  27 #include "AliAODCluster.h"
  28 #include "AliVCaloCells.h"
  29 #include "AliOADBContainer.h"
  30 #include "AliPHOSGeometry.h"
  31 #include "AliAnalysisManager.h"
  32
  33
  34
  35 #include "AliPHOSClusterSelection.h"
  36
  37 AliPHOSClusterSelection::AliPHOSClusterSelection()
  38   : fMinChargedParticleTrackDistance(-1.),
  39     fNotUnfolded(false),
  40     fMaxDispR2(-1.),
  41     fIsCore(false),
  42     fMaxTOF(-1.)
  43 {
  44   // Defaults to the most lenient selection allowable
  45         return;
  46 }
  47
  48 AliPHOSClusterSelection::~AliPHOSClusterSelection()
  49 {
  50 }
  51
  52 Bool_t AliPHOSClusterSelection::IsSelected(AliVCluster* cluster) const
  53 {
  54   return IsSelectedCPV(cluster)
  55     && IsSelectedUnfolded(cluster)
  56     && IsSelectedDisp(cluster)
  57     && IsSelectedTOF(cluster);
  58 }
  59
  60 Bool_t AliPHOSClusterSelection::IsSelectedCPV(AliVCluster* cluster) const
  61 {
  62   if( 0 > fMinChargedParticleTrackDistance )//No selection on CPV
  63     return true;
  64
  65   else{
  66     Bool_t cpvBit = true; //Changed to false if there is a matching track in the requested radius.
  67
  68     Double_t dx=cluster->GetTrackDx();//Position for Track matching
  69     Double_t dz=cluster->GetTrackDz();
  70
  71
  72     AliESDEvent* EventESD = dynamic_cast<AliESDEvent*> (AliPHOSClusterSelection::GetCurrentEvent());//dynamic cast to test for ESD or AOD
  73     AliAODEvent* EventAOD = dynamic_cast<AliAODEvent*> (AliPHOSClusterSelection::GetCurrentEvent());
  74     Double_t mf = 0.; //
  75     if(EventAOD) mf = EventAOD->GetMagneticField(); //Positive for ++ and negative for --
  76     else if(EventESD) mf = EventESD->GetMagneticField(); //Positive for ++ and negative for --
  77
  78
  79     if(EventESD){//if ESD
  80       AliESDCaloCluster * ESDcluster = static_cast<AliESDCaloCluster*> (cluster);//Know its ESD so static cast is fine?
  81       TArrayI * itracks = ESDcluster-> GetTracksMatched() ;
  82       if(itracks->GetSize()>0){
  83         Int_t iTr = itracks->At(0);
  84         if(iTr>=0 && iTr<EventESD->GetNumberOfTracks()){
  85           AliVParticle* track = EventESD->GetTrack(iTr);
  86           Double_t pt = track->Pt() ;
  87           Short_t charge = track->Charge() ;
  88           Double_t r = AliPHOSClusterSelection::TestCPV(dx, dz, pt, charge, mf);
  89           cpvBit=(r>fMinChargedParticleTrackDistance) ;
  90         }
  91       }
  92     }
  93
  94
  95     if(EventAOD){//if AOD
  96       AliAODCluster * AODcluster = static_cast<AliAODCluster*> (cluster);
  97       int nTracksMatched = AODcluster->GetNTracksMatched();
  98       if(nTracksMatched > 0) {
  99         AliVTrack* track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0));
 100         if ( track ) {
 101           Double_t pt = track->Pt();
 102           Short_t charge = track->Charge();
 103           Double_t r = AliPHOSClusterSelection::TestCPV(dx, dz, pt, charge, mf);
 104           cpvBit=(r>fMinChargedParticleTrackDistance) ;
 105         }
 106       }
 107     }
 108     return cpvBit;
 109   }
 110 }
 111
 112
 113 Bool_t AliPHOSClusterSelection::IsSelectedUnfolded(AliVCluster* cluster) const
 114 {
 115   if(!fNotUnfolded)
 116     return true;
 117   else{
 118     Bool_t NotUnfolded = cluster->GetNExMax()<2;//True if it was not unfolded
 119     return !NotUnfolded;
 120   }
 121 }
 122
 123 Bool_t AliPHOSClusterSelection::IsSelectedDisp(AliVCluster* cluster) const
 124 {
 125   if(0 > fMaxDispR2)
 126     return true;
 127   else{
 128     Double_t m02 = 0.,m20 = 0.;
 129     if(!fIsCore){//No core calculation
 130       m02 = cluster->GetM02();
 131       m20 = cluster->GetM20();
 132       return AliPHOSClusterSelection::TestLambda(cluster->E(),m20,m02) ;
 133     }
 134     else{//DispCore
 135       AliVCaloCells* cells = static_cast<AliVCaloCells*> (AliPHOSClusterSelection::GetCurrentEvent()->GetPHOSCells());//Need the cells
 136       AliPHOSClusterSelection::EvalCoreLambdas(cluster, cells, m02, m20);
 137       return AliPHOSClusterSelection::TestLambda(cluster->E(),m20,m02);
 138     }
 139   }
 140
 141 }
 142
 143 Bool_t AliPHOSClusterSelection::IsSelectedTOF(AliVCluster* cluster) const
 144 {
 145   if(0 > fMaxTOF)
 146     return true;
 147   else{
 148     // Time of Flight (TOF)
 149     Double_t tof = cluster->GetTOF();//Time of Flight for the cluster
 150     return  TMath::Abs(tof) < fMaxTOF;//True if the cut is passed
 151   }
 152 }
 153
 154 AliPHOSClusterSelection* AliPHOSClusterSelection::SetMinChargedParticleTrackDistance(Float_t distance)
 155 {
 156   // 'distance' set the minimal allowable distance between the cluster
 157   // and the nearest extrapolated track.
 158   // If 'distance' is negative, then all clusters are sellected, the selection
 159   // being "not applied" or "disabled".
 160
 161   fMinChargedParticleTrackDistance = distance;
 162   return this;
 163 }
 164
 165 AliPHOSClusterSelection* AliPHOSClusterSelection::SetNotUnfolded(Bool_t notUnfolded)
 166 {
 167   //if notUnfolded true, it rejects Unfolded Clusters
 168   fNotUnfolded = notUnfolded;
 169   return this;
 170 }
 171
 172 AliPHOSClusterSelection* AliPHOSClusterSelection::SetMaxDispR2(Float_t maxR2)
 173 {
 174   // 'maxR2' sets the maximum allowed dispersion.
 175   // If 'maxR2' is negative, then all clusters are selected, the selection
 176   // being "not applied" or "disabled".
 177
 178   fMaxDispR2 = maxR2;
 179   return this;
 180 }
 181
 182 AliPHOSClusterSelection* AliPHOSClusterSelection::SetIsCore(Bool_t isCore)
 183 {
 184   // 'isCore' sets wether core version of Disp is used. -1 gives no core.
 185
 186   fIsCore = isCore;
 187   return this;
 188 }
 189
 190 AliPHOSClusterSelection* AliPHOSClusterSelection::SetMaxTOF(Float_t maxTOF)
 191 {
 192   // 'maxTOF' sets the maximum allowed time of flight for the cluster.
 193   // If 'maxTOF' is negative, all clusters are selected and the selection is "disabled".
 194
 195   fMaxTOF = maxTOF;
 196   return this;
 197 }
 198
 199 TString AliPHOSClusterSelection::ToString() const
 200 {
 201   // returns a string an quasi-unique string for whatever selection
 202   // parameters the instance contains. The uniqueness of the string
 203   // is limited by the precision given in the formatting of the string.
 204   // Take care that the precision is sufficient for your needs.
 205
 206   TString string ="v1";
 207   string+=Form("_%.2f", fMinChargedParticleTrackDistance);
 208   string+=Form("_%i", fNotUnfolded);
 209   string+=Form("_%.2f", fMaxDispR2);
 210   string+=Form("_%i", fIsCore);
 211   string+=Form("_%.2f", fMaxTOF);
 212   return string;
 213 }
 214
 215
 216 Float_t AliPHOSClusterSelection::GetMinChargedParticleTrackDistance(const TString& string)
 217 {
 218   TObjArray* array = string.Tokenize("_");
 219   TObjString* objString = (TObjString*) array->At(kMinChargedParticleTrackDistance);
 220   Float_t flt = objString->String().Atof();
 221   delete array;
 222   return flt;
 223 }
 224
 225 Bool_t AliPHOSClusterSelection::GetUnfolded(const TString& string)
 226 {
 227   TObjArray* array = string.Tokenize("_");
 228   TObjString* objString = (TObjString*) array->At(kNotUnfolded);
 229   Float_t flt = objString->String().Atoi();
 230   delete array;
 231   return flt;
 232 }
 233
 234 Float_t AliPHOSClusterSelection::GetMaxDispR2(const TString& string)
 235 {
 236   TObjArray* array = string.Tokenize("_");
 237   TObjString* objString = (TObjString*) array->At(kMaxDispR2);
 238   Float_t flt = objString->String().Atof();
 239   delete array;
 240   return flt;
 241 }
 242
 243 Bool_t AliPHOSClusterSelection::GetIsCore(const TString& string)
 244 {
 245   TObjArray* array = string.Tokenize("_");
 246   TObjString* objString = (TObjString*) array->At(kIsCore);
 247   Float_t flt = objString->String().Atoi();
 248   delete array;
 249   return flt;
 250 }
 251
 252 Float_t AliPHOSClusterSelection::GetMaxTOF(const TString& string)
 253 {
 254   TObjArray* array = string.Tokenize("_");
 255   TObjString* objString = (TObjString*) array->At(kMaxTOF);
 256   Float_t flt = objString->String().Atof();
 257   delete array;
 258   return flt;
 259 }
 260
 261
 262 Double_t AliPHOSClusterSelection::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge, Double_t mf) const {
 263   //Parameterization of LHC10h period
 264   //_true if neutral_
 265   //Copied from Pi0Flow task
 266
 267   Double_t meanX=0;
 268   Double_t meanZ=0.;
 269   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
 270                          6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
 271   Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
 272
 273   if(mf<0.){ //field --
 274     meanZ = -0.468318 ;
 275     if(charge>0)
 276       meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
 277     else
 278       meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
 279   }
 280   else{ //Field ++
 281     meanZ= -0.468318;
 282     if(charge>0)
 283       meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
 284     else
 285       meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
 286   }
 287
 288   Double_t rz=(dz-meanZ)/sz ;
 289   Double_t rx=(dx-meanX)/sx ;
 290   return TMath::Sqrt(rx*rx+rz*rz) ;
 291 }
 292
 293 //____________________________________________________________________________
 294 void  AliPHOSClusterSelection::EvalCoreLambdas(AliVCluster  * clu, AliVCaloCells * cells,Double_t &m02, Double_t &m20) const
 295 {
 296   //calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum
 297   //Copied from pi0flowtask
 298   Int_t fRunNumber;
 299   const Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
 300   // Calculates the center of gravity in the local PHOS-module coordinates
 301   Float_t wtot = 0;
 302   const Int_t mulDigit=clu->GetNCells() ;
 303   Double_t xc[mulDigit] ;
 304   Double_t zc[mulDigit] ;
 305   Double_t wi[mulDigit] ;
 306   Double_t x = 0 ;
 307   Double_t z = 0 ;
 308   const Double_t logWeight=4.5 ;
 309   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
 310     Int_t relid[4] ;
 311     Float_t xi=0. ;
 312     Float_t zi=0. ;
 313     Int_t absId = clu->GetCellAbsId(iDigit) ;
 314
 315
 316     AliESDEvent* EventESD = dynamic_cast<AliESDEvent*> (AliPHOSClusterSelection::GetCurrentEvent());//dynamic cast to test for ESD or AOD
 317
 318     AliAODEvent* EventAOD = dynamic_cast<AliAODEvent*> (AliPHOSClusterSelection::GetCurrentEvent());
 319     if(EventESD)
 320       fRunNumber = EventESD->GetRunNumber() ;
 321     if(EventAOD)
 322       fRunNumber = EventESD->GetRunNumber() ;
 323
 324     AliOADBContainer geomContainer("phosGeo");//Initialize Geometry
 325     geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
 326     TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
 327     AliPHOSGeometry * fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP") ;
 328     for(Int_t mod=0; mod<5; mod++) {
 329       if(!matrixes->At(mod)) {
 330         AliInfo(Form("No PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
 331         continue;
 332       }
 333
 334
 335       fPHOSGeo->AbsToRelNumbering(absId, relid) ;
 336       fPHOSGeo->RelPosInModule(relid, xi, zi);
 337       xc[iDigit]=xi ;
 338       zc[iDigit]=zi ;
 339       Double_t ei = elist[iDigit]*cells->GetCellAmplitude(absId) ;
 340       wi[iDigit]=0. ;
 341       if (clu->E()>0 && ei>0) {
 342         wi[iDigit] = TMath::Max( 0., logWeight + TMath::Log( ei / clu->E() ) ) ;
 343         Double_t w=wi[iDigit];
 344         x    += xc[iDigit] * w ;
 345         z    += zc[iDigit] * w ;
 346         wtot += w ;
 347       }
 348     }
 349     if (wtot>0) {
 350       x /= wtot ;
 351       z /= wtot ;
 352     }
 353
 354     wtot = 0. ;
 355     Double_t dxx  = 0.;
 356     Double_t dzz  = 0.;
 357     Double_t dxz  = 0.;
 358     Double_t xCut = 0. ;
 359     Double_t zCut = 0. ;
 360     for(Int_t iDigit1=0; iDigit1<mulDigit; iDigit1++) {
 361       Double_t w=wi[iDigit1];
 362       if (w>0.) {
 363         Double_t xi1= xc[iDigit1] ;
 364         Double_t zi1= zc[iDigit1] ;
 365         if((xi1-x)*(xi1-x)+(zi1-z)*(zi1-z) < 4.5*4.5){
 366           xCut += w * xi1 ;
 367           zCut += w * zi1 ;
 368           dxx  += w * xi1 * xi1 ;
 369           dzz  += w * zi1 * zi1 ;
 370           dxz  += w * xi1 * zi1 ;
 371           wtot += w ;
 372         }
 373       }
 374
 375     }
 376     if (wtot>0) {
 377       xCut/= wtot ;
 378       zCut/= wtot ;
 379       dxx /= wtot ;
 380       dzz /= wtot ;
 381       dxz /= wtot ;
 382       dxx -= xCut * xCut ;
 383       dzz -= zCut * zCut ;
 384       dxz -= xCut * zCut ;
 385
 386       m02 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
 387       m20 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
 388     }
 389     else {
 390       m20=m02=0.;
 391     }
 392
 393   }
 394 }
 395
 396 //_____________________________________________________________________________
 397 Bool_t AliPHOSClusterSelection::TestLambda(Double_t pt,Double_t l1,Double_t l2) const
 398 {
 399   //Evaluates if lambdas correspond to photon cluster
 400   //Tuned using pp data
 401   //copied from Pi0FlowTask
 402   Double_t l2Mean, l1Mean, l2Sigma, l1Sigma, c, R2;
 403   if(! fIsCore ){
 404     l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
 405     l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
 406     l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
 407     l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
 408     c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
 409     R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
 410       0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
 411       0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
 412   }
 413   else{
 414     //For core radius R=4.5
 415     l1Mean  = 1.150200 + 0.097886/(1.+1.486645*pt+0.000038*pt*pt) ;
 416     l2Mean = 1.574706 + 0.997966*exp(-0.895075*pt)-0.010666*pt ;
 417     l1Sigma = 0.100255 + 0.337177*exp(-0.517684*pt)+0.001170*pt ;
 418     l2Sigma = 0.232580 + 0.573401*exp(-0.735903*pt)-0.002325*pt ;
 419     c = -0.110983 -0.017353/(1.-1.836995*pt+0.934517*pt*pt) ;
 420     R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
 421       0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
 422       0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
 423   }
 424   return (R2<fMaxDispR2*fMaxDispR2);//TODO: Check if this should be 2. order
 425 }
 426
 427
 428 AliVEvent* AliPHOSClusterSelection::GetCurrentEvent() const
 429 {
 430   // Hackish way of getting the current event.
 431   // Its probably not appropriate to call this function outside of
 432   // AliAnalysisTaskSE::UserExec
 433
 434   AliAnalysisManager* analysisManager = dynamic_cast<AliAnalysisManager*>(AliAnalysisManager::GetAnalysisManager());
 435   AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(analysisManager->GetInputEventHandler());
 436   AliMultiInputEventHandler *multiInputHandler = dynamic_cast<AliMultiInputEventHandler *>(inputHandler);
 437   if (multiInputHandler)
 438     inputHandler = dynamic_cast<AliInputEventHandler *>(multiInputHandler->GetFirstInputEventHandler());
 439
 440   AliVEvent* inputEvent = dynamic_cast<AliVEvent*>(inputHandler->GetEvent());
 441   if( ! inputEvent )
 442     AliError("Was not able to retrieve event!");
 443
 444   return inputEvent;
 445 }