PWG3/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 "AliAODpidUtil.h"
  29 #include "AliAODTrack.h"
  30 #include "AliESDpid.h"
  31
  32 #include "AliHFEdetPIDqa.h"
  33 //#include "AliHFEpidEMCAL.h"
  34 #include "AliHFEpidQAmanager.h"
  35
  36 //#include "AliEMCALRecoUtils.h"
  37 //#include "AliEMCALGeometry.h"
  38 #include "AliVCluster.h"
  39 #include "AliVCaloCells.h"
  40 #include "AliVEvent.h"
  41 #include "AliLog.h"
  42 #include "AliEMCALPIDUtils.h"
  43 #include "AliPID.h"
  44 #include "AliESDEvent.h"
  45 #include "AliESDtrack.h"
  46 //#include "AliEMCALTrack.h"
  47 //#include "AliEMCALTracker.h"
  48 #include "AliHFEpidEMCAL.h"
  49
  50 ClassImp(AliHFEpidEMCAL)
  51
  52 //___________________________________________________________________
  53 AliHFEpidEMCAL::AliHFEpidEMCAL():
  54   AliHFEpidBase()
  55   , fPID(NULL)
  56   , feopMim(0.8)
  57   , feopMax(1.4)
  58 {
  59   //
  60   // Constructor
  61   //
  62 }
  63
  64 //___________________________________________________________________
  65 AliHFEpidEMCAL::AliHFEpidEMCAL(const Char_t *name):
  66   AliHFEpidBase(name)
  67   , fPID(NULL)
  68   , feopMim(0.8)
  69   , feopMax(1.4)
  70 {
  71   //
  72   // Constructor
  73   //
  74 }
  75
  76 //___________________________________________________________________
  77 AliHFEpidEMCAL::AliHFEpidEMCAL(const AliHFEpidEMCAL &c):
  78   AliHFEpidBase("")
  79   , fPID(NULL)
  80   , feopMim(0.8)
  81   , feopMax(1.4)
  82 {
  83   //
  84   // Copy operator
  85   //
  86
  87   c.Copy(*this);
  88 }
  89 //___________________________________________________________________
  90 AliHFEpidEMCAL &AliHFEpidEMCAL::operator=(const AliHFEpidEMCAL &ref){
  91   //
  92   // Assignment operator
  93   //
  94
  95   if(this != &ref){
  96     ref.Copy(*this);
  97   }
  98
  99   return *this;
 100 }
 101 //___________________________________________________________________
 102 AliHFEpidEMCAL::~AliHFEpidEMCAL(){
 103   //
 104   // Destructor
 105   //
 106   if(fPID) delete fPID;
 107 }
 108 //___________________________________________________________________
 109 void AliHFEpidEMCAL::Copy(TObject &ref) const {
 110   //
 111   // Performs the copying of the object
 112   //
 113   AliHFEpidEMCAL &target = dynamic_cast<AliHFEpidEMCAL &>(ref);
 114
 115   target.fPID = fPID;
 116   target.feopMax = feopMax;
 117   target.feopMim = feopMim;
 118
 119   AliHFEpidBase::Copy(ref);
 120 }
 121 //___________________________________________________________________
 122 Bool_t AliHFEpidEMCAL::InitializePID(){
 123   //
 124   // InitializePID: EMCAL experts have to implement code here
 125   //
 126   return kTRUE;
 127 }
 128
 129 //___________________________________________________________________
 130 Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager */*pidqa*/) const
 131 { //Function to a return a code indicating whether or not an electron candidate is selected
 132   //
 133   //
 134   if((!fESDpid && track->IsESDanalysis()) || (!fAODpid && track->IsAODanalysis())) return 0;
 135   AliDebug(2, "PID object available");
 136   // EMCal not fESDpid  (s.s Feb. 11)
 137
 138   /*
 139   AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 140   if(!((dynamic_cast<const AliVTrack *>(track->GetRecTrack()))->GetStatus() & AliESDtrack::kEMCALpid)) return 0;
 141   AliDebug(2, "Track Has EMCAL PID");
 142   */
 143   //if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID);
 144   // not QA for EMCal, will be added (s.s Feb. 11)
 145
 146   Double_t eop = 0.;//MomentumEnergyMatch(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
 147   AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop));
 148   Int_t pdg = 0;
 149   if(eop>feopMim && eop<feopMax){
 150     pdg = 11;
 151     //if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID);
 152   }
 153   else
 154   {
 155    pdg = 211; // EMCal doesn't separate pi,k.p by e/p. return pion code as same as TRD
 156   }
 157
 158     printf("eID %g ; %d \n",eop,pdg);
 159
 160   return pdg;
 161 }
 162 /*
 163 Double_t AliHFEpidEMCAL::MomentumEnergyMatch(const AliVParticle *track) const
 164 { // Returns e/p if an electron is matched
 165
 166   Float_t  clsPos[3];
 167   Double_t trkPos[3];
 168   Double_t matchclsE = 9999.9;
 169
 170   const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 171   AliESDEvent *evt = esdtrack->GetESDEvent();
 172   Double_t  magF = evt->GetMagneticField();
 173   Double_t magSign = 1.0;
 174   if(magF<0)magSign = -1.0;
 175   //printf("magF ; %g ; %g \n", magF,magSign);
 176
 177   if (!TGeoGlobalMagField::Instance()->GetField()) {
 178           printf("Loading field map...\n");
 179           //AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
 180           AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
 181           TGeoGlobalMagField::Instance()->SetField(field);
 182                     }
 183
 184   AliEMCALTrack *emctrack = new AliEMCALTrack(*esdtrack);
 185   Double_t fieldB[3];
 186   emctrack->GetBxByBz(fieldB);
 187   //printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
 188
 189   for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
 190    {
 191
 192    AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
 193    if(!cluster->IsEMCAL()) continue;
 194
 195    cluster->GetPosition(clsPos);
 196    if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  continue;
 197    emctrack->GetXYZ(trkPos);
 198
 199    TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
 200    TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
 201
 202    Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi();  // track cluster matching
 203    Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
 204
 205    double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
 206
 207    if(rmatch<0.02)
 208       {
 209        matchclsE = cluster->E();
 210       }
 211   }
 212    delete emctrack;
 213
 214    double feop = -9999.9;
 215    if(matchclsE<9999) feop = matchclsE/esdtrack->P();
 216
 217    //   if(feop!=-9999.9)printf("%f\n",feop) ;
 218    return feop;
 219
 220 }
 221 */
 222