PWG4/AliGammaMCDataReader.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.1.2.1  2007/07/26 10:32:09  schutz
  22  * new analysis classes in the the new analysis framework
  23  *
  24  *
  25  */
  26
  27 //_________________________________________________________________________
  28 // Class for reading data (Kinematics and ESDs) in order to do prompt gamma correlations
  29 //  Class created from old AliPHOSGammaJet
  30 //  (see AliRoot versions previous Release 4-09)
  31 //
  32 //*-- Author: Gustavo Conesa (LNF-INFN)
  33 //////////////////////////////////////////////////////////////////////////////
  34
  35
  36 // --- ROOT system ---
  37 #include "Riostream.h"
  38 #include <TParticle.h>
  39
  40 //---- ANALYSIS system ----
  41 #include "AliGammaMCDataReader.h"
  42 #include "AliLog.h"
  43 #include "AliESDEvent.h"
  44 #include "AliESDVertex.h"
  45 #include "AliESDCaloCluster.h"
  46 #include "AliStack.h"
  47
  48 ClassImp(AliGammaMCDataReader)
  49
  50 //____________________________________________________________________________
  51   AliGammaMCDataReader::AliGammaMCDataReader() :
  52     AliGammaReader(),
  53     fEMCALPID(0),fPHOSPID(0),
  54     fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.), fPHOSPhotonWeight(0.),  fPHOSPi0Weight(0.)
  55 {
  56   //Default Ctor
  57
  58   //Initialize parameters
  59   fDataType=kMCData;
  60   InitParameters();
  61
  62 }
  63
  64 //____________________________________________________________________________
  65 AliGammaMCDataReader::AliGammaMCDataReader(const AliGammaMCDataReader & g) :
  66   AliGammaReader(g),
  67   fEMCALPID(g.fEMCALPID),
  68   fPHOSPID(g.fPHOSPID),
  69   fEMCALPhotonWeight(g.fEMCALPhotonWeight),
  70   fEMCALPi0Weight(g.fEMCALPi0Weight),
  71   fPHOSPhotonWeight(g.fPHOSPhotonWeight),
  72   fPHOSPi0Weight(g.fPHOSPi0Weight)
  73 {
  74   // cpy ctor
  75 }
  76
  77 //_________________________________________________________________________
  78 AliGammaMCDataReader & AliGammaMCDataReader::operator = (const AliGammaMCDataReader & source)
  79 {
  80   // assignment operator
  81
  82   if(&source == this) return *this;
  83
  84   fEMCALPID = source.fEMCALPID ;
  85   fPHOSPID = source.fPHOSPID ;
  86   fEMCALPhotonWeight = source. fEMCALPhotonWeight ;
  87   fEMCALPi0Weight = source.fEMCALPi0Weight ;
  88   fPHOSPhotonWeight = source.fPHOSPhotonWeight ;
  89   fPHOSPi0Weight = source.fPHOSPi0Weight ;
  90
  91   return *this;
  92
  93 }
  94
  95 //____________________________________________________________________________
  96 void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
  97                                               TClonesArray * plCTS,
  98                                               TClonesArray * plEMCAL,
  99                                               TClonesArray * plPHOS, TClonesArray *){
 100
 101   //Create a list of particles from the ESD. These particles have been measured
 102   //by the Central Tracking system (TPC+ITS), PHOS and EMCAL
 103
 104   AliESDEvent* esd = (AliESDEvent*) data;
 105   AliStack* stack = (AliStack*) kine;
 106
 107   Int_t npar  = 0 ;
 108   Float_t *pid = new Float_t[AliPID::kSPECIESN];
 109   AliDebug(3,"Fill particle lists");
 110
 111   //Get vertex for momentum calculation
 112   Double_t v[3] ; //vertex ;
 113   esd->GetVertex()->GetXYZ(v) ;
 114
 115   //########### PHOS ##############
 116
 117   Int_t begphos = esd->GetFirstPHOSCluster();
 118   Int_t endphos = esd->GetFirstPHOSCluster() +
 119     esd->GetNumberOfPHOSClusters() ;
 120   Int_t indexPH = plPHOS->GetEntries() ;
 121   AliDebug(3,Form("First PHOS particle %d, last particle %d", begphos,endphos));
 122
 123   for (npar =  begphos; npar <  endphos; npar++) {//////////////PHOS track loop
 124     AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
 125     if(clus->GetTrackMatched()==-1 ){
 126       //Create a TParticle to fill the particle list
 127       TLorentzVector momentum ;
 128       clus->GetMomentum(momentum, v);
 129       pid=clus->GetPid();
 130       Int_t pdg = 22;
 131
 132       if(fPHOSPID){
 133         if( pid[AliPID::kPhoton] > fPHOSPhotonWeight) pdg=22;
 134         if( pid[AliPID::kPi0] > fPHOSPi0Weight) pdg=111;
 135         else pdg = 0 ;
 136       }
 137
 138       TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
 139                                            momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
 140
 141       AliDebug(4,Form("PHOS clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
 142
 143       //Look  if parent is prompt photon
 144       Int_t label = clus->GetLabel();
 145       if(label>=0){
 146         TParticle * pmother = stack->Particle(label);
 147         Int_t imother = pmother->GetFirstMother();
 148         if(imother == 6 || imother == 7)
 149           pmother->SetFirstMother(22);
 150       }
 151
 152       new((*plPHOS)[indexPH++])   TParticle(*particle) ;
 153
 154     }//not charged
 155   }//cluster loop
 156
 157   //########### CTS (TPC+ITS) #####################
 158   Int_t begtpc   = 0 ;
 159   Int_t endtpc   = esd->GetNumberOfTracks() ;
 160   Int_t indexCh  = plCTS->GetEntries() ;
 161   AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
 162
 163   for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
 164     AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
 165
 166     //We want tracks fitted in the detectors:
 167     ULong_t status=AliESDtrack::kTPCrefit;
 168     status|=AliESDtrack::kITSrefit;
 169
 170     //We want tracks whose PID bit is set:
 171     //     ULong_t status =AliESDtrack::kITSpid;
 172     //     status|=AliESDtrack::kTPCpid;
 173
 174     if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
 175       // Do something with the tracks which were successfully
 176       // re-fitted
 177       Double_t en = 0; //track ->GetTPCsignal() ;
 178       Double_t mom[3];
 179       track->GetPxPyPz(mom) ;
 180       Double_t px = mom[0];
 181       Double_t py = mom[1];
 182       Double_t pz = mom[2]; //Check with TPC people if this is correct.
 183       Int_t pdg = 11; //Give any charged PDG code, in this case electron.
 184       //I just want to tag the particle as charged
 185        TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
 186                                             px, py, pz, en, v[0], v[1], v[2], 0);
 187
 188       //TParticle * particle = new TParticle() ;
 189       //particle->SetMomentum(px,py,pz,en) ;
 190
 191       new((*plCTS)[indexCh++])       TParticle(*particle) ;
 192     }
 193   }
 194
 195   //################ EMCAL ##############
 196
 197   Int_t indexEM  = plEMCAL->GetEntries() ;
 198   Int_t begem = esd->GetFirstEMCALCluster();
 199   Int_t endem = esd->GetFirstEMCALCluster() +
 200     esd->GetNumberOfEMCALClusters() ;
 201
 202   AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
 203
 204   for (npar =  begem; npar <  endem; npar++) {//////////////EMCAL track loop
 205     AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
 206     Int_t clustertype= clus->GetClusterType();
 207     if(clustertype == AliESDCaloCluster::kClusterv1 && clus->GetTrackMatched()==-1 ){
 208
 209       TLorentzVector momentum ;
 210       clus->GetMomentum(momentum, v);
 211       pid=clus->GetPid();
 212       Int_t pdg = 22;
 213
 214       if(fEMCALPID){
 215         if( pid[AliPID::kPhoton] > fEMCALPhotonWeight) pdg=22;
 216         else if( pid[AliPID::kPi0] > fEMCALPi0Weight) pdg=111;
 217         else pdg = 0;
 218       }
 219
 220       TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
 221                                            momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
 222       AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
 223
 224       //Look  if parent is prompt photon
 225       Int_t label = clus->GetLabel();
 226       if(label>=0){
 227         TParticle * pmother = stack->Particle(label);
 228         Int_t imother = pmother->GetFirstMother();
 229         if(imother == 6 || imother == 7)
 230           pmother->SetFirstMother(22);
 231       }
 232
 233       new((*plEMCAL)[indexEM++])   TParticle(*particle) ;
 234
 235     }//not charged, not pseudocluster
 236   }//Cluster loop
 237
 238   AliDebug(3,"Particle lists filled");
 239
 240 }
 241
 242   //____________________________________________________________________________
 243 void AliGammaMCDataReader::InitParameters()
 244 {
 245
 246   //Initialize the parameters of the analysis.
 247
 248   //Fill particle lists when PID is ok
 249   fEMCALPID = kFALSE;
 250   fPHOSPID = kFALSE;
 251   fEMCALPhotonWeight = 0.5 ;
 252   fEMCALPi0Weight = 0.5 ;
 253   fPHOSPhotonWeight = 0.8 ;
 254   fPHOSPi0Weight = 0.5 ;
 255
 256 }
 257
 258
 259 void AliGammaMCDataReader::Print(const Option_t * opt) const
 260 {
 261
 262   //Print some relevant parameters set for the analysis
 263   if(! opt)
 264     return;
 265
 266   Info("Print", "%s %s", GetName(), GetTitle() ) ;
 267   printf("PHOS PID on?               =     %d\n",  fPHOSPID) ;
 268   printf("EMCAL PID  on?         =     %d\n",  fEMCALPID) ;
 269   printf("PHOS PID weight , photon  %f, pi0 %f\n\n",  fPHOSPhotonWeight,  fPHOSPi0Weight) ;
 270   printf("EMCAL PID weight, photon %f, pi0 %f\n",   fEMCALPhotonWeight,  fEMCALPi0Weight) ;
 271
 272 }