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