PWGHF/hfe/AliHFEpidEMCAL.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 // Class for EMCAL PID
  17 // Implements the abstract base class AliHFEpidBase
  18 // IsInitialized() does the PID decision with energy and
  19 // momentum matching (e/p)
  20 //
  21 // Authors:
  22 //  Shingo Sakai
  23 //
  24 //
  25 //
  26 #include <TMath.h>
  27 #include "AliESDInputHandler.h"
  28 #include "AliESDpid.h"
  29
  30 #include "AliHFEdetPIDqa.h"
  31 #include "AliHFEpidEMCAL.h"
  32 #include "AliHFEpidQAmanager.h"
  33
  34 #include "AliHFEemcalPIDqa.h"
  35 //#include "AliVCluster.h"
  36 //#include "AliVCaloCells.h"
  37 //#include "AliVEvent.h"
  38 #include "AliLog.h"
  39 #include "AliPID.h"
  40 //#include "AliESDEvent.h"
  41 //#include "AliESDtrack.h"
  42 #include "AliHFEpidEMCAL.h"
  43
  44 ClassImp(AliHFEpidEMCAL)
  45
  46 //___________________________________________________________________
  47 AliHFEpidEMCAL::AliHFEpidEMCAL():
  48   AliHFEpidBase()
  49   , fPID(NULL)
  50   , feopMim(0.9)
  51   , feopMax(1.3)
  52 {
  53   //
  54   // Constructor
  55   //
  56 }
  57
  58 //___________________________________________________________________
  59 AliHFEpidEMCAL::AliHFEpidEMCAL(const Char_t *name):
  60   AliHFEpidBase(name)
  61   , fPID(NULL)
  62   , feopMim(0.9)
  63   , feopMax(1.3)
  64 {
  65   //
  66   // Constructor
  67   //
  68 }
  69
  70 //___________________________________________________________________
  71 AliHFEpidEMCAL::AliHFEpidEMCAL(const AliHFEpidEMCAL &c):
  72   AliHFEpidBase("")
  73   , fPID(NULL)
  74   , feopMim(0.9)
  75   , feopMax(1.3)
  76 {
  77   //
  78   // Copy operator
  79   //
  80
  81   c.Copy(*this);
  82 }
  83 //___________________________________________________________________
  84 AliHFEpidEMCAL &AliHFEpidEMCAL::operator=(const AliHFEpidEMCAL &ref){
  85   //
  86   // Assignment operator
  87   //
  88
  89   if(this != &ref){
  90     ref.Copy(*this);
  91   }
  92
  93   return *this;
  94 }
  95 //___________________________________________________________________
  96 AliHFEpidEMCAL::~AliHFEpidEMCAL(){
  97   //
  98   // Destructor
  99   //
 100   if(fPID) delete fPID;
 101 }
 102 //___________________________________________________________________
 103 void AliHFEpidEMCAL::Copy(TObject &ref) const {
 104   //
 105   // Performs the copying of the object
 106   //
 107   AliHFEpidEMCAL &target = dynamic_cast<AliHFEpidEMCAL &>(ref);
 108
 109   target.fPID = fPID;
 110   target.feopMax = feopMax;
 111   target.feopMim = feopMim;
 112
 113   AliHFEpidBase::Copy(ref);
 114 }
 115 //___________________________________________________________________
 116 Bool_t AliHFEpidEMCAL::InitializePID(Int_t /*run*/){
 117   //
 118   // InitializePID: EMCAL experts have to implement code here
 119   //
 120   return kTRUE;
 121 }
 122
 123 //___________________________________________________________________
 124 Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
 125 { //Function to a return a code indicating whether or not an electron candidate is selected
 126   //
 127   //
 128
 129   if(track==NULL)return 0;
 130
 131   if(!fkPIDResponse) return 0;
 132   AliDebug(2, "PID object available");
 133   // EMCal not fESDpid  (s.s Feb. 11)
 134
 135   const AliVTrack *trk = dynamic_cast<const AliVTrack *>(track->GetRecTrack());
 136   if (trk == NULL) return 0;
 137   //AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 138   if(!(trk->GetStatus() & AliESDtrack::kEMCALmatch)) return 0;
 139   AliDebug(2, "Track Has EMCAL PID");
 140
 141   if(pidqa)
 142   pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID);
 143   // not QA for EMCal, will be added (s.s Feb. 11)
 144
 145
 146      double feopMimCut = 0.0;
 147      double feopMaxCut = 0.0;
 148
 149      feopMimCut = feopMim;
 150      feopMaxCut = feopMax;
 151
 152      if(feopMax >900)  // not use fix cut
 153        {
 154         feopMimCut = CalEopCutMim(track->GetRecTrack(),0);
 155         feopMaxCut = CalEopCutMax(track->GetRecTrack(),0);
 156        }
 157      if(feopMax < -900)  // not use fix cut for MC
 158        {
 159         feopMimCut = CalEopCutMim(track->GetRecTrack(),1);
 160         feopMaxCut = CalEopCutMax(track->GetRecTrack(),1);
 161        }
 162
 163  //printf("eop cuts org; %g; %g \n",feopMim,feopMax);
 164  //printf("eop cuts ; %g; %g \n",feopMimCut,feopMaxCut);
 165
 166   Double_t eop = MomentumEnergyMatchV2(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
 167   AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop));
 168   AliDebug(2, Form("E/p cut ; %g ; %g \n", feopMim,feopMax));
 169   Int_t pdg = 0;
 170   //if(eop>feopMim && eop<feopMax){
 171   //if(eop>feopMim && eop<feopMax){
 172   if(eop>feopMimCut && eop<feopMaxCut){
 173     pdg = 11;
 174     //printf("eop cuts ; %g; %g ; %g\n",feopMimCut,feopMaxCut,eop);
 175     if(pidqa)
 176     pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID);
 177   }
 178   else
 179   {
 180    pdg = 211; // EMCal doesn't separate pi,k.p by e/p. return pion code as same as TRD
 181   }
 182
 183     AliDebug(1, Form("eID %g ; %d \n",eop,pdg));
 184     //printf("eID %g ; %d \n",eop,pdg);
 185
 186   return pdg;
 187 }
 188
 189
 190 //__________________________________________________________________________
 191 Double_t AliHFEpidEMCAL::CalEopCutMax(const AliVParticle *const track, Int_t flageop) const
 192 {
 193  // max nSig cuts for 10d
 194
 195   double maxCut = 0.0;
 196   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 197   if(esdtrack==NULL)return maxCut;
 198   double pt = esdtrack->Pt();
 199
 200   if(flageop<0.5)
 201     { //<--- new para for updated non-liniarity
 202      double meanP[3] = {1.08833e+00,1.17837e+00,7.47923e-02};
 203      double sigP[3] = {9.93208e-05,7.13074e-04,2.45145e-04};
 204      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt);
 205      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
 206      maxCut = mean+3.0*sig;
 207     }
 208   else
 209     {
 210      double meanP[3] = {9.94183e-01,1.23952e+00,3.55571e-01};
 211      double sigP[3] = {4.61207e-02,-3.39978e-02,4.26198e-01};
 212      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt);
 213      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
 214      maxCut = mean+3.0*sig;
 215     }
 216
 217   //printf("eop cuts org; %g; %g ; %g\n",feopMim,feopMax,pt);
 218
 219      return maxCut;
 220 }
 221
 222 Double_t AliHFEpidEMCAL::CalEopCutMim(const AliVParticle *const track, Int_t flageop) const
 223 {
 224  // mim nsig cuts for 10d
 225   double mimCut = 0.0;
 226   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 227   if(esdtrack==NULL)return mimCut;
 228   double pt = esdtrack->Pt();
 229
 230   if(flageop<0.5) // real
 231      {
 232       //printf("real");
 233      double meanP[3] = {1.08833e+00,1.17837e+00,7.47923e-02};
 234      double sigP[3] = {9.93208e-05,7.13074e-04,2.45145e-04};
 235      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt);
 236      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
 237      mimCut = mean-3.0*sig;
 238     }
 239    else // MC
 240     {
 241       //printf("MC");
 242      double meanP[3] = {9.94183e-01,1.23952e+00,3.55571e-01};
 243      double sigP[3] = {4.61207e-02,-3.39978e-02,4.26198e-01};
 244      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt);
 245      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
 246      mimCut = mean-3.0*sig;
 247     }
 248      return mimCut;
 249 }
 250
 251 //___________________________________________________________________________
 252 Double_t AliHFEpidEMCAL::MomentumEnergyMatchV2(const AliVParticle *const track) const
 253 { // Returns e/p & Rmatch
 254
 255   Double_t matchclsE = 9999.9;
 256   double feop = -9999.9;
 257
 258   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 259   if(esdtrack==NULL)return feop;
 260   const AliESDEvent *evt = esdtrack->GetESDEvent();
 261
 262    Int_t icl = esdtrack->GetEMCALcluster();
 263
 264    AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
 265    if(cluster==NULL) {return feop;}
 266    if(!cluster->IsEMCAL()) {return feop;}
 267
 268        matchclsE = cluster->E();
 269
 270    if(matchclsE<9999.0) feop = matchclsE/esdtrack->P();
 271
 272    return feop;
 273
 274 }
 275
 276
 277 /*
 278 //____________________________________________________________________________________
 279 Double_t AliHFEpidEMCAL::MomentumEnergyMatchV1(const AliVParticle *track) const
 280 { // Returns e/p if an electron is matched
 281
 282   Float_t  clsPos[3];
 283   Double_t trkPos[3];
 284   Double_t matchclsE = 9999.9;
 285
 286   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 287   AliESDEvent *evt = esdtrack->GetESDEvent();
 288   Double_t  magF = evt->GetMagneticField();
 289   Double_t magSign = 1.0;
 290   if(magF<0)magSign = -1.0;
 291   //printf("magF ; %g ; %g \n", magF,magSign);
 292
 293   if (!TGeoGlobalMagField::Instance()->GetField()) {
 294           printf("Loading field map...\n");
 295           //AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
 296           AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
 297           TGeoGlobalMagField::Instance()->SetField(field);
 298                     }
 299
 300   AliEMCALTrack *emctrack = new AliEMCALTrack(*esdtrack);
 301   Double_t fieldB[3];
 302   emctrack->GetBxByBz(fieldB);
 303   //printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
 304
 305   for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
 306    {
 307
 308    AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
 309    if(!cluster->IsEMCAL()) continue;
 310
 311    cluster->GetPosition(clsPos);
 312    if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  continue;
 313    emctrack->GetXYZ(trkPos);
 314
 315    TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
 316    TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
 317
 318    Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi();  // track cluster matching
 319    Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
 320
 321    double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
 322
 323    if(rmatch<0.02)
 324       {
 325        matchclsE = cluster->E();
 326       }
 327   }
 328    delete emctrack;
 329
 330    double feop = -9999.9;
 331    if(matchclsE<9999) feop = matchclsE/esdtrack->P();
 332
 333    //   if(feop!=-9999.9)printf("%f\n",feop) ;
 334    return feop;
 335
 336 }
 337 */
 338