PWG4/JCORRAN/AliJPhoton.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 notifce   *
  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: AliJPhoton.cxx,v 1.2 2008/01/21 11:56:39 djkim Exp $
  17
  18 ////////////////////////////////////////////////////
  19 //
  20 //  \file AliJPhoton.cxx
  21 //  \brief
  22 //  \author J. Rak, D.J.Kim, R.Diaz (University of Jyvaskyla)
  23 //  \email: djkim@jyu.fi
  24 //  \version $Revision: 1.1 $
  25 //  \date $Date: 2008/05/02 11:56:39 $
  26 //
  27 //  class for alice photons
  28 ////////////////////////////////////////////////////
  29
  30 #include <iostream>
  31 #include <fstream>
  32 #include <stdlib.h>
  33 #include <stdio.h>
  34
  35 #include <TFormula.h>
  36 #include "JConst.h"
  37 #include "AliPhJBaseTrack.h"
  38 #include "AliPhJPhoton.h"
  39 #include "AliJPhoton.h"
  40
  41
  42 //ClassImp(AliJPhoton)
  43
  44
  45 //______________________________________________________________________________
  46 AliJPhoton::AliJPhoton() :
  47   AliPhJPhoton(),
  48   fCaloType(kUndef),
  49   fDistToBadChannel(-999),
  50   fDispersion(-999),
  51   fM20(-999),
  52   fM02(-999),
  53   fEmcCpvDist(-999),
  54   fNCells(0),
  55   fSuperModuleId(-999),
  56   fCellsAbsId(0x0),
  57   fCellsAmpFraction(0x0),
  58   fMCid(-999)
  59 {
  60   // default constructor
  61   for(Int_t i=0;i<kUnknownAli+1;i++) fCaloPID[i]=-1;
  62
  63   SetPID((Double_t*)NULL);
  64 }
  65
  66 //______________________________________________________________________________
  67 AliJPhoton::AliJPhoton(const AliJPhoton& a):
  68   AliPhJPhoton(a),
  69   fCaloType(a.fCaloType),
  70   fDistToBadChannel(a.fDistToBadChannel),
  71   fDispersion(a.fDispersion),
  72   fM20(a.fM20),
  73   fM02(a.fM02),
  74   fEmcCpvDist(a.fEmcCpvDist),
  75   fNCells(a.fNCells),
  76   fSuperModuleId(a.fSuperModuleId),
  77   fCellsAbsId(a.fCellsAbsId),
  78   fCellsAmpFraction(a.fCellsAmpFraction),
  79   fMCid(a.fMCid)
  80 {
  81   //copy contructor
  82   for(Int_t i=0;i<kUnknownAli+1;i++) fCaloPID[i] = a.fCaloPID[i];
  83 }
  84
  85 //______________________________________________________________________________
  86 void  AliJPhoton::SetCellsAbsId(const UShort_t *array)
  87 {
  88     //  Set the array of cell absId numbers
  89     if (fNCells) {
  90         if(fCellsAbsId){ delete [] fCellsAbsId; fCellsAbsId = NULL; }
  91         fCellsAbsId = new  UShort_t[fNCells];
  92         for (Int_t i = 0; i < fNCells; i++) fCellsAbsId[i] = array[i];
  93     }
  94 }
  95
  96 //_______________________________________________________________________
  97 void  AliJPhoton::SetCellsAmplitudeFraction(const Double32_t *array)
  98 {
  99     //  Set the array of cell amplitude fraction
 100     if (fNCells) {
 101         if(fCellsAmpFraction){ delete [] fCellsAmpFraction; fCellsAmpFraction = NULL;}
 102         fCellsAmpFraction = new  Double32_t[fNCells];
 103         for (Int_t i = 0; i < fNCells; i++) fCellsAmpFraction[i] = array[i];
 104     }
 105 }
 106
 107 //______________________________________________________________________________
 108 void  AliJPhoton::SetPID(const Double_t *pid) {
 109    //set pid
 110     if(pid){
 111        for(Int_t i=0; i<kUnknownAli+1; ++i) fCaloPID[i]=pid[i];
 112        SetProbPhot(fCaloPID[kPhotonAli]);
 113     }else{
 114       for(Int_t i=0; i<kUnknownAli+1; fCaloPID[i++]=0.){}
 115       fCaloPID[kUnknownAli]=1.;
 116     }
 117 }
 118
 119
 120 //______________________________________________________________________________
 121 particleType AliJPhoton::GetParticleType() {
 122
 123     // return the most problable particle type
 124     // Note: following the AliCaloPID implementation
 125
 126     //Init default weights
 127     Float_t wPhoton = 0.75 ;
 128     Float_t wPi0 = 0.8 ;
 129     Float_t wElectron = 0.5 ;
 130     Float_t wCharged = 0.5 ;
 131     Float_t wNeutral = 0.5 ;
 132
 133     Bool_t usePHOSweightFormula = kTRUE;
 134     //Formula to set the PID weight threshold for photon or pizero
 135     // ALICE-INT-2005-016 (2007)
 136     TFormula* wPhotonPHOSFormula =
 137     new TFormula("photonWeight","0.75*(x<40)+ 0.68*(x>=100)+(x>=40 && x<100)*(0.98+x*(6e-3)-x*x*(2e-04)+x*x*x*(1.1e-06))");
 138     TFormula* wPi0PHOSFormula =
 139     new TFormula("pi0Weight","0.80*(x<65)+ 0.915*(x>=100)+(x>=65 && x-x*(1.95e-3)-x*x*(4.31e-05)+x*x*x*(3.61e-07))");
 140
 141     if(fCaloType == kPHOSCalo && usePHOSweightFormula){
 142        wPhoton  = wPhotonPHOSFormula->Eval(GetE()) ;
 143        wPi0 = wPi0PHOSFormula->Eval(GetE());
 144       // cout << "wPhotonF " << wPhoton << " wphoton " << fCaloPID[kPhotonAli] << " wPi0F " << wPi0 << " wpi0 " << fCaloPID[kPi0Ali] << endl;
 145     }
 146
 147     if(fCaloType == kEMCALCalo){
 148        wPhoton   =  0.8 ;
 149        wPi0      =  0.5 ;
 150        wElectron =  0.8 ;
 151        wCharged  =  0.5 ;
 152        wNeutral  =  0.5 ;
 153     }
 154
 155
 156    particleType pdg = kHadron;
 157
 158   //Select most probable ID
 159   if(fCaloType == kPHOSCalo){
 160     if(fCaloPID[kPhotonAli] > wPhoton) pdg = kPhoton ;
 161     else if(fCaloPID[kPi0Ali] > wPi0) pdg = kPizero ;
 162     //else if(fCaloPID[kElectronAli] > wElectron)  pdg = electron ;
 163     //else if(fCaloPID[kEleConAli] >  wElectron) pdg = electronCon ;
 164     //else if(chargedHadronWeight > wCharged) pdg = chargedHadron ;
 165     //else if(neutralHadronWeight > wNeutral) pdg = neutralHadron ;
 166     //else if(allChargedWeight >  allNeutralWeight)
 167     // pdg = chargedUnknown ;
 168     //else
 169     //  pdg = neutralUnknown ;
 170   }
 171   else{//EMCAL
 172     //Temporal solution, electrons and photons not differenciated
 173     if(fCaloPID[kPhotonAli] + fCaloPID[kElectronAli]  > wPhoton) pdg = kPhoton ;
 174     else if(fCaloPID[kPi0Ali] > wPi0) pdg = kPizero ;
 175     //else if(chargedHadronWeight + neutralHadronWeight > wCharged) pdg = chargedHadron ;
 176     //else if(neutralHadronWeight + chargedHadronWeight > wNeutral) pdg = neutralHadron ;
 177     //else pdg =  neutralUnknown ;
 178
 179   }
 180
 181   delete wPhotonPHOSFormula;
 182   delete wPi0PHOSFormula;
 183
 184   //if(fCaloType == kPHOSCalo) cout << "pdg " << pdg << endl;
 185
 186   return pdg ;
 187
 188 }
 189
 190 //_____________________________________________________________________________
 191 AliJPhoton&  AliJPhoton::operator=(const  AliJPhoton& photon){
 192   //overloaded =
 193   if(this != &photon) {
 194     AliPhJPhoton::operator=(photon);
 195     fCaloType = photon.fCaloType;
 196     for(Int_t i=0; i<kUnknownAli+1; i++){
 197       fCaloPID[i] = photon.fCaloPID[i];
 198     }
 199     fDistToBadChannel = photon.fDistToBadChannel;
 200     fDispersion    = photon.fDispersion;
 201     fM20           = photon.fM20;
 202     fM02           = photon.fM02;
 203     fEmcCpvDist    = photon.fEmcCpvDist;
 204     fNCells        = photon.fNCells;
 205     fSuperModuleId =  photon.fSuperModuleId;
 206     fCellsAbsId    = photon.fCellsAbsId;
 207     fCellsAmpFraction = photon.fCellsAmpFraction;
 208     fMCid          = photon.fMCid;
 209   }
 210   return *this;
 211 }
 212