PWG4/AliGammaDataReader.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16 /* $Id$ */
  17
  18 /* History of cvs commits:
  19  *
  20  * $Log$
  21  * Revision 1.2  2007/08/17 12:40:04  schutz
  22  * New analysis classes by Gustavo Conesa
  23  *
  24  * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
  25  * new analysis classes in the the new analysis framework
  26  *
  27  *
  28  */
  29
  30 //_________________________________________________________________________
  31 // Class for reading data (ESDs) in order to do prompt gamma correlations
  32 //  Class created from old AliPHOSGammaJet
  33 //  (see AliRoot versions previous Release 4-09)
  34 //
  35 //*-- Author: Gustavo Conesa (LNF-INFN)
  36 //////////////////////////////////////////////////////////////////////////////
  37
  38
  39 // --- ROOT system ---
  40 #include "Riostream.h"
  41 #include <TParticle.h>
  42
  43 //---- ANALYSIS system ----
  44 #include "AliGammaDataReader.h"
  45 #include "AliLog.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDVertex.h"
  48 #include "AliESDCaloCluster.h"
  49
  50 ClassImp(AliGammaDataReader)
  51
  52 //____________________________________________________________________________
  53 AliGammaDataReader::AliGammaDataReader() :
  54   AliGammaReader()
  55 {
  56   //Default Ctor
  57
  58   //Initialize parameters
  59   fDataType=kData;
  60
  61 }
  62
  63 //____________________________________________________________________________
  64 AliGammaDataReader::AliGammaDataReader(const AliGammaDataReader & g) :
  65   AliGammaReader(g)
  66 {
  67   // cpy ctor
  68 }
  69
  70 //_________________________________________________________________________
  71 AliGammaDataReader & AliGammaDataReader::operator = (const AliGammaDataReader & source)
  72 {
  73   // assignment operator
  74
  75   if(&source == this) return *this;
  76
  77   return *this;
  78
  79 }
  80
  81 //____________________________________________________________________________
  82 void AliGammaDataReader::CreateParticleList(TObject * data, TObject *,
  83                                             TClonesArray * plCTS,
  84                                             TClonesArray * plEMCAL,
  85                                             TClonesArray * plPHOS,
  86                                             TClonesArray * ,
  87                                             TClonesArray * ,
  88                                             TClonesArray * ){
  89
  90   //Create a list of particles from the ESD. These particles have been measured
  91   //by the Central Tracking system (TPC+ITS), PHOS and EMCAL
  92   //Also create particle list with mothers.
  93
  94   AliESDEvent* esd = (AliESDEvent*) data;
  95
  96   Int_t npar  = 0 ;
  97   Double_t *pid = new Double_t[AliPID::kSPECIESN];
  98    AliDebug(3,"Fill particle lists");
  99
 100   //Get vertex for momentum calculation
 101   Double_t v[3] ; //vertex ;
 102   esd->GetVertex()->GetXYZ(v) ;
 103
 104   //########### PHOS ##############
 105
 106   Int_t begphos = 0;
 107   Int_t endphos = esd->GetNumberOfCaloClusters() ;
 108   Int_t indexPH = plPHOS->GetEntries() ;
 109   Int_t begem = 0;
 110   Int_t endem = endphos;
 111
 112   AliDebug(3,Form("First Calorimeter particle %d, last particle %d", begphos,endphos));
 113
 114   for (npar =  begphos; npar <  endphos; npar++) {//////////////PHOS track loop
 115     AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
 116
 117     //Check if it is PHOS cluster
 118     if(!clus->IsEMCAL())
 119       begem=npar+1;
 120     else
 121       continue;
 122
 123     AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
 124     if(clus->GetTrackMatched()==-1){
 125       //Create a TParticle to fill the particle list
 126       TLorentzVector momentum ;
 127       clus->GetMomentum(momentum, v);
 128       Double_t phi = momentum.Phi();
 129       if(phi<0) phi+=TMath::TwoPi() ;
 130       if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fPHOSEtaCut &&
 131          phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] ) {
 132
 133         pid=clus->GetPid();
 134         Int_t pdg = 22;
 135
 136         if(IsPHOSPIDOn()){
 137           AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
 138                           momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
 139                           pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
 140
 141           Float_t wPhoton =  fPHOSPhotonWeight;
 142           Float_t wPi0 =  fPHOSPi0Weight;
 143
 144           if(fPHOSWeightFormula){
 145             wPhoton = fPHOSPhotonWeightFormula->Eval(momentum.E()) ;
 146             wPi0 =    fPHOSPi0WeightFormula->Eval(momentum.E());
 147           }
 148
 149           if(pid[AliPID::kPhoton] > wPhoton)
 150             pdg = kPhoton ;
 151           else if(pid[AliPID::kPi0] > wPi0)
 152             pdg = kPi0 ;
 153           else if(pid[AliPID::kElectron] > fPHOSElectronWeight)
 154             pdg = kElectron ;
 155           else if(pid[AliPID::kEleCon] > fPHOSElectronWeight)
 156             pdg = kEleCon ;
 157           else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fPHOSChargeWeight)
 158             pdg = kChargedHadron ;
 159           else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fPHOSNeutralWeight)
 160             pdg = kNeutralHadron ;
 161
 162           else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >
 163                   pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron])
 164             pdg = kChargedUnknown  ;
 165           else
 166             pdg = kNeutralUnknown ;
 167           //neutral cluster, unidentifed.
 168         }
 169
 170         if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown ){//keep only neutral particles in the array
 171           TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
 172                                                momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
 173
 174           AliDebug(4,Form("PHOS added: pdg %d, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
 175
 176           new((*plPHOS)[indexPH++])   TParticle(*particle) ;
 177         }
 178         else AliDebug(4,Form("PHOS charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f\n",
 179                              pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
 180
 181       }//pt, eta, phi cut
 182       else      AliDebug(4,"Particle not added");
 183     }//not charged
 184   }//cluster loop
 185
 186   //########### CTS (TPC+ITS) #####################
 187   Int_t begtpc   = 0 ;
 188   Int_t endtpc   = esd->GetNumberOfTracks() ;
 189   Int_t indexCh  = plCTS->GetEntries() ;
 190   AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
 191
 192   for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
 193     AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
 194
 195     //We want tracks fitted in the detectors:
 196     ULong_t status=AliESDtrack::kTPCrefit;
 197     status|=AliESDtrack::kITSrefit;
 198
 199     //We want tracks whose PID bit is set:
 200     //     ULong_t status =AliESDtrack::kITSpid;
 201     //     status|=AliESDtrack::kTPCpid;
 202
 203     if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
 204       // Do something with the tracks which were successfully
 205       // re-fitted
 206       Double_t en = 0; //track ->GetTPCsignal() ;
 207       Double_t mom[3];
 208       track->GetPxPyPz(mom) ;
 209       Double_t px = mom[0];
 210       Double_t py = mom[1];
 211       Double_t pz = mom[2]; //Check with TPC people if this is correct.
 212       Int_t pdg = 11; //Give any charged PDG code, in this case electron.
 213       //I just want to tag the particle as charged
 214        TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
 215                                             px, py, pz, en, v[0], v[1], v[2], 0);
 216
 217       //TParticle * particle = new TParticle() ;
 218       //particle->SetMomentum(px,py,pz,en) ;
 219        if(particle->Pt() > fChargedPtCut && TMath::Abs(particle->Eta())<fCTSEtaCut)
 220          new((*plCTS)[indexCh++])       TParticle(*particle) ;
 221     }
 222   }
 223
 224   //################ EMCAL ##############
 225
 226   Int_t indexEM  = plEMCAL->GetEntries() ;
 227
 228   for (npar =  begem; npar <  endem; npar++) {//////////////EMCAL track loop
 229     AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
 230     Int_t clustertype= clus->GetClusterType();
 231     AliDebug(4,Form("EMCAL clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
 232
 233     if(clustertype == AliESDCaloCluster::kEMCALClusterv1 && clus->GetTrackMatched()==-1 ){
 234
 235       TLorentzVector momentum ;
 236       clus->GetMomentum(momentum, v);
 237       Double_t phi = momentum.Phi();
 238       if(phi<0) phi+=TMath::TwoPi() ;
 239       if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fEMCALEtaCut &&
 240         phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] ) {
 241
 242         pid=clus->GetPid();
 243         Int_t pdg = 22;
 244
 245         if(IsEMCALPIDOn()){
 246           AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
 247                           momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
 248                           pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
 249
 250           if(pid[AliPID::kPhoton] > fEMCALPhotonWeight)
 251             pdg = kPhoton ;
 252           else if(pid[AliPID::kPi0] > fEMCALPi0Weight)
 253             pdg = kPi0 ;
 254           else if(pid[AliPID::kElectron] > fEMCALElectronWeight)
 255             pdg = kElectron ;
 256           else if(pid[AliPID::kEleCon] > fEMCALElectronWeight)
 257             pdg = kEleCon ;
 258           else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fEMCALChargeWeight)
 259             pdg = kChargedHadron ;
 260           else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fEMCALNeutralWeight)
 261             pdg = kNeutralHadron ;
 262           else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >
 263                   pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron])
 264             pdg = kChargedUnknown ;
 265           else
 266             pdg = kNeutralUnknown ;
 267         }
 268
 269         if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown){//keep only neutral particles in the array
 270
 271           TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
 272                                                momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
 273           AliDebug(4,Form("EMCAL cluster added: pdg %f, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
 274
 275           new((*plEMCAL)[indexEM++])   TParticle(*particle) ;
 276         }
 277         else AliDebug(4,Form("EMCAL charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f",
 278                              pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
 279
 280       }//pt, phi, eta cut
 281       else      AliDebug(4,"Particle not added");
 282     }//not charged, not pseudocluster
 283   }//Cluster loop
 284
 285   AliDebug(3,Form("Particle lists filled, PHOS clusters %d, EMCAL clusters %d, CTS tracks %d",
 286                   indexPH,indexEM,indexCh));
 287
 288 }
 289