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.5  2007/12/14 10:54:15  gustavo
  22  * Protection against cluster labels larger than kinematics entries
  23  *
  24  * Revision 1.4  2007/10/29 13:48:42  gustavo
  25  * Corrected coding violations
  26  *
  27  * Revision 1.2  2007/08/17 12:40:04  schutz
  28  * New analysis classes by Gustavo Conesa
  29  *
  30  * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
  31  * new analysis classes in the the new analysis framework
  32  *
  33  *
  34  */
  35
  36 //_________________________________________________________________________
  37 // Class for reading data (Kinematics and ESDs) in order to do prompt gamma correlations
  38 //  Class created from old AliPHOSGammaJet
  39 //  (see AliRoot versions previous Release 4-09)
  40 //
  41 //*-- Author: Gustavo Conesa (LNF-INFN)
  42 //////////////////////////////////////////////////////////////////////////////
  43
  44
  45 // --- ROOT system ---
  46 #include <TFormula.h>
  47 #include <TParticle.h>
  48
  49 //---- ANALYSIS system ----
  50 #include "AliGammaMCDataReader.h"
  51 #include "AliESDEvent.h"
  52 #include "AliESDVertex.h"
  53 #include "AliESDCaloCluster.h"
  54 #include "AliStack.h"
  55 #include "AliLog.h"
  56
  57 ClassImp(AliGammaMCDataReader)
  58
  59 //____________________________________________________________________________
  60   AliGammaMCDataReader::AliGammaMCDataReader() :
  61     AliGammaReader()
  62 {
  63   //Default Ctor
  64
  65   //Initialize parameters
  66   fDataType=kMCData;
  67
  68 }
  69
  70 //____________________________________________________________________________
  71 AliGammaMCDataReader::AliGammaMCDataReader(const AliGammaMCDataReader & g) :
  72   AliGammaReader(g)
  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
  85   return *this;
  86
  87 }
  88
  89 //____________________________________________________________________________
  90 void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
  91                                               TClonesArray * plCTS,
  92                                               TClonesArray * plEMCAL,
  93                                               TClonesArray * plPHOS,
  94                                               TClonesArray * plPrimCTS,
  95                                               TClonesArray * plPrimEMCAL,
  96                                               TClonesArray * plPrimPHOS){
  97
  98   //Create a list of particles from the ESD. These particles have been measured
  99   //by the Central Tracking system (TPC+ITS+...), PHOS and EMCAL
 100   //Also create particle list with mothers.
 101
 102   AliESDEvent* esd = (AliESDEvent*) data;
 103   AliStack* stack = (AliStack*) kine;
 104
 105   Int_t npar  = 0 ;
 106   Double_t *pid = new Double_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   //########### CALORIMETERS ##############
 114   Int_t nCaloCluster = esd->GetNumberOfCaloClusters() ;
 115   Int_t indexPH = plPHOS->GetEntries() ;
 116   Int_t indexEM = plEMCAL->GetEntries() ;
 117
 118   for (npar =  0; npar <  nCaloCluster; npar++) {//////////////CaloCluster loop
 119     AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve cluster from esd
 120     Int_t type = clus->GetClusterType();
 121
 122     //########### PHOS ##############
 123     if(fSwitchOnPHOS && type ==  AliESDCaloCluster::kPHOSCluster){
 124       AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
 125
 126       if(clus->GetTrackMatched()==-1){
 127         TLorentzVector momentum ;
 128         clus->GetMomentum(momentum, v);
 129         Double_t phi = momentum.Phi();
 130         if(phi<0) phi+=TMath::TwoPi() ;
 131         if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fPHOSEtaCut &&
 132            phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] ) {
 133
 134           pid=clus->GetPid();
 135           Int_t pdg = 22;
 136
 137           if(IsPHOSPIDOn()){
 138             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 ",
 139                             momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
 140                             pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
 141
 142             Float_t wPhoton =  fPHOSPhotonWeight;
 143             Float_t wPi0 =  fPHOSPi0Weight;
 144
 145             if(fPHOSWeightFormula){
 146               wPhoton = fPHOSPhotonWeightFormula->Eval(momentum.E()) ;
 147               wPi0 =    fPHOSPi0WeightFormula->Eval(momentum.E());
 148             }
 149
 150             if(pid[AliPID::kPhoton] > wPhoton)
 151               pdg = kPhoton ;
 152             else if(pid[AliPID::kPi0] > wPi0)
 153               pdg = kPi0 ;
 154             else if(pid[AliPID::kElectron] > fPHOSElectronWeight)
 155               pdg = kElectron ;
 156             else if(pid[AliPID::kEleCon] > fPHOSElectronWeight)
 157               pdg = kEleCon ;
 158             else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fPHOSChargeWeight)
 159               pdg = kChargedHadron ;
 160             else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fPHOSNeutralWeight)
 161               pdg = kNeutralHadron ;
 162
 163             else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >
 164                     pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron])
 165               pdg = kChargedUnknown  ;
 166             else
 167               pdg = kNeutralUnknown ;
 168             //neutral cluster, unidentifed.
 169           }
 170
 171           if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown ){//keep only neutral particles in the array
 172
 173             TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
 174             new((*plPHOS)[indexPH])   TParticle(*particle) ;
 175             AliDebug(4,Form("PHOS added: pdg %d, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
 176
 177             //###############
 178             //Check kinematics
 179             //###############
 180             TParticle * pmother = new TParticle();
 181             Int_t label = clus->GetLabel();
 182             if(label < stack->GetNtrack())
 183               pmother = GetMotherParticle(label,stack, "PHOS",momentum);
 184             else
 185               AliInfo(Form("PHOS: Bad label %d, too large, NPrimaries %d",label,stack->GetNtrack()));
 186             new((*plPrimPHOS)[indexPH])   TParticle(*pmother) ;
 187
 188             indexPH++;
 189           }
 190           else AliDebug(4,Form("PHOS charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f\n",
 191                                pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
 192
 193         }//pt, eta, phi cut
 194         else    AliDebug(4,"Particle not added");
 195       }//track-match?
 196     }//PHOS cluster
 197
 198     //################ EMCAL ##############
 199     else if(fSwitchOnEMCAL &&  type ==  AliESDCaloCluster::kEMCALClusterv1){
 200       AliDebug(4,Form("EMCAL clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
 201
 202       if(clus->GetTrackMatched()==-1 ){
 203         TLorentzVector momentum ;
 204         clus->GetMomentum(momentum, v);
 205         Double_t phi = momentum.Phi();
 206         if(phi<0) phi+=TMath::TwoPi() ;
 207         if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fEMCALEtaCut &&
 208            phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] ) {
 209
 210           pid=clus->GetPid();
 211           Int_t pdg = 22;
 212
 213           if(IsEMCALPIDOn()){
 214             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 ",
 215                             momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
 216                             pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
 217
 218             if(pid[AliPID::kPhoton] > fEMCALPhotonWeight)
 219               pdg = kPhoton ;
 220             else if(pid[AliPID::kPi0] > fEMCALPi0Weight)
 221               pdg = kPi0 ;
 222             else if(pid[AliPID::kElectron] > fEMCALElectronWeight)
 223               pdg = kElectron ;
 224             else if(pid[AliPID::kEleCon] > fEMCALElectronWeight)
 225               pdg = kEleCon ;
 226             else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fEMCALChargeWeight)
 227               pdg = kChargedHadron ;
 228             else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fEMCALNeutralWeight)
 229               pdg = kNeutralHadron ;
 230             else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >
 231                     pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron])
 232               pdg = kChargedUnknown ;
 233             else
 234               pdg = kNeutralUnknown ;
 235           }
 236
 237           if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown){//keep only neutral particles in the array
 238
 239             TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
 240             new((*plEMCAL)[indexEM])   TParticle(*particle) ;
 241             AliDebug(4,Form("EMCAL cluster added: pdg %f, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
 242
 243             //###############
 244             //Check kinematics
 245             //###############
 246             TParticle * pmother = new TParticle();
 247             Int_t label = clus->GetLabel();
 248             if(label < stack->GetNtrack())
 249               pmother = GetMotherParticle(label,stack, "EMCAL",momentum);
 250             else
 251               AliInfo(Form("EMCAL: Bad label %d, too large, NPrimaries %d",label,stack->GetNtrack()));
 252             new((*plPrimEMCAL)[indexEM])   TParticle(*pmother) ;
 253
 254             indexEM++;
 255           }
 256           else AliDebug(4,Form("EMCAL charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f",
 257                                pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
 258
 259         }//pt, phi, eta cut
 260         else    AliDebug(4,"Particle not added");
 261       }//track-matched
 262     }//EMCAL cluster
 263
 264   }//cluster loop
 265
 266
 267
 268   //########### CTS (TPC+ITS) #####################
 269   Int_t nTracks   = esd->GetNumberOfTracks() ;
 270   Int_t indexCh  = plCTS->GetEntries() ;
 271
 272   if(fSwitchOnCTS){
 273     AliDebug(3,Form("Number of tracks %d",nTracks));
 274
 275     for (npar =  0; npar <  nTracks; npar++) {////////////// track loop
 276       AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
 277
 278       //We want tracks fitted in the detectors:
 279       ULong_t status=AliESDtrack::kTPCrefit;
 280       status|=AliESDtrack::kITSrefit;
 281
 282       //We want tracks whose PID bit is set:
 283       //     ULong_t status =AliESDtrack::kITSpid;
 284       //     status|=AliESDtrack::kTPCpid;
 285
 286       if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
 287         // Do something with the tracks which were successfully
 288         // re-fitted
 289         Double_t en = 0; //track ->GetTPCsignal() ;
 290         Double_t mom[3];
 291         track->GetPxPyPz(mom) ;
 292         Double_t px = mom[0];
 293         Double_t py = mom[1];
 294         Double_t pz = mom[2]; //Check with TPC people if this is correct.
 295         Int_t pdg = 11; //Give any charged PDG code, in this case electron.
 296         //I just want to tag the particle as charged
 297
 298         //Check kinematics
 299         Int_t label = TMath::Abs(track->GetLabel());
 300         TParticle * pmother = new TParticle();
 301         pmother = stack->Particle(label);
 302
 303         TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
 304                                              px, py, pz, en, v[0], v[1], v[2], 0);
 305
 306         if(particle->Pt() > fChargedPtCut && TMath::Abs(particle->Eta())<fCTSEtaCut){
 307           new((*plCTS)[indexCh])       TParticle(*particle) ;
 308           new((*plPrimCTS)[indexCh])       TParticle(*pmother) ;
 309           indexCh++;
 310
 311         }//kinematic selection
 312       }//select track from refit
 313     }//track loop
 314   }//CTS
 315
 316   AliDebug(3,Form("Particle lists filled, tracks  %d , clusters: EMCAL %d, PHOS %d", indexCh,indexEM,indexPH));
 317
 318 }
 319
 320 TParticle * AliGammaMCDataReader:: GetMotherParticle(Int_t label, AliStack *stack, TString calo,  TLorentzVector momentum)
 321 {
 322   //Gets the primary particle and do some checks:
 323   //Check if primary is inside calorimeter and look the mother outsie
 324   //Check if mother is a decay photon, in which case check if decay was overlapped
 325
 326   Float_t minangle = 0;
 327   Float_t ipdist = 0;
 328   TParticle * pmother = new TParticle();
 329
 330   if(calo == "PHOS"){
 331     ipdist= fPHOSIPDistance;
 332     minangle = fPHOSMinAngle;
 333   }
 334   else if (calo == "EMCAL"){
 335     ipdist = fEMCALIPDistance;
 336     minangle = fEMCALMinAngle;
 337   }
 338
 339
 340   if(label>=0){
 341     pmother = stack->Particle(label);
 342     Int_t mostatus = pmother->GetStatusCode();
 343     Int_t mopdg    = TMath::Abs(pmother->GetPdgCode());
 344
 345     //Check if mother is a conversion inside the calorimeter
 346     //In such case, return the mother before the calorimeter.
 347     //First approximation.
 348     Float_t vy = TMath::Abs(pmother->Vy()) ;
 349
 350     if( mostatus == 0 && vy >= ipdist){
 351
 352       //cout<<"Calo: "<<calo<<" vy "<<vy<<" E clus "<<momentum.E()<<" Emother "<<pmother->Energy()<<" "
 353       //  <<pmother->GetName()<<endl;
 354
 355       while( vy >= ipdist){//inside calorimeter
 356         AliDebug(4,Form("Conversion inside calorimeter, mother vertex %0.2f, IP distance %0.2f", vy, ipdist));
 357         pmother =  stack->Particle(pmother->GetMother(0));
 358         vy = TMath::Abs(pmother->Vy()) ;
 359         //cout<<" label "<<label<<" Mother: "<<pmother->GetName()<<" E "<<pmother->Energy()<<" Status "<<pmother->GetStatusCode()<<"  and vertex "<<vy<<endl;
 360         mostatus = pmother->GetStatusCode();
 361         mopdg    = TMath::Abs(pmother->GetPdgCode());
 362       }//while vertex is inside calorimeter
 363       //cout<<"Calo: "<<calo<<" final vy "<<vy<<" E clus "<<momentum.E()<<" Emother "<<pmother->Energy()<<" "
 364       //  <<pmother->GetName()<<endl;
 365     }//check status and vertex
 366
 367     AliDebug(4,Form("%s, ESD E %2.2f, PID %d,  mother: E %2.2f, label %d, status %d,  vertex %3.2f, mother 2 %d, grandmother %d \n",
 368                     calo.Data(),momentum.E(),pmother->Energy(), label, pmother->GetPdgCode(),
 369                     pmother->GetStatusCode(), vy, pmother->GetMother(0), stack->GetPrimary(label)));
 370
 371     //Check Decay photons
 372     if(mopdg == 22){
 373
 374       //his mother was a pi0?
 375       TParticle * pmotherpi0 =  stack->Particle(pmother->GetMother(0));
 376       if( pmotherpi0->GetPdgCode() == 111){
 377
 378         AliDebug(4,Form(" %s: E cluster %2.2f, E gamma %2.2f, Mother Pi0, E %0.2f, status %d, daughters %d\n",
 379                         calo.Data(), momentum.E(),pmother->Energy(),pmotherpi0->Energy(),
 380                         pmotherpi0->GetStatusCode(), pmotherpi0->GetNDaughters()));
 381
 382         if(pmotherpi0->GetNDaughters() == 1) mostatus = 2 ; //signal that this photon can only be decay, not overlapped.
 383         else if(pmotherpi0->GetNDaughters() == 2){
 384
 385           TParticle * pd1 = stack->Particle(pmotherpi0->GetDaughter(0));
 386           TParticle * pd2 = stack->Particle(pmotherpi0->GetDaughter(1));
 387           //if(pmotherpi0->Energy()> 10 ){
 388 //        cout<<"Two "<<calo<<" daugthers, pi0 label "<<pmother->GetMother(0)<<" E :"<<pmotherpi0->Energy()<<endl;
 389 //        cout<<" 1) label "<<pmotherpi0->GetDaughter(0)<<" pdg "<<pd1->GetPdgCode()<<" E  "<<pd1->Energy()
 390 //            <<" phi "<<pd1->Phi()*TMath::RadToDeg()<<" eta "<<pd1->Eta()
 391 //            <<" mother label "<<pd1->GetMother(0)<<" n daug "<<pd1->GetNDaughters() <<endl;
 392 //          cout<<" 2) label "<<pmotherpi0->GetDaughter(1)<<" pdg "<<pd2->GetPdgCode()<<" E  "<<pd2->Energy()
 393 //              <<" phi "<<pd2->Phi()*TMath::RadToDeg()<<" eta "<<pd2->Eta()<<" mother label "
 394 //              <<pd2->GetMother(0)<<" n daug "<<pd2->GetNDaughters() <<endl;
 395             //}
 396           if(pd1->GetPdgCode() == 22 && pd2->GetPdgCode() == 22){
 397             TLorentzVector lv1 , lv2 ;
 398             pd1->Momentum(lv1);
 399             pd2->Momentum(lv2);
 400             Double_t angle = lv1.Angle(lv2.Vect());
 401 //          if(pmotherpi0->Energy()> 10 )
 402 //            cout<<"angle*ipdist "<<angle*ipdist<<" mindist "<< minangle <<endl;
 403             if (angle < minangle){
 404               //There is overlapping, pass the pi0 as mother
 405
 406               AliDebug(4,Form(">>>> %s cluster with E %2.2f is a overlapped pi0, E %2.2f, angle %2.4f < anglemin %2.4f\n",
 407                               calo.Data(), momentum.E(), pmotherpi0->Energy(), angle, minangle));
 408
 409               pmother = pmotherpi0 ;
 410
 411             }
 412             else mostatus = 2 ; // gamma decay not overlapped
 413           }// daughters are photons
 414           else mostatus = 2; // daughters are e-gamma or e-e, no overlapped, or charged cluster
 415         }//2 daughters
 416         else AliDebug(4,Form("pi0 has %d daughters",pmotherpi0->GetNDaughters()));
 417       }//pi0 decay photon?
 418     }//photon
 419
 420     pmother->SetStatusCode(mostatus); // if status = 2, decay photon.
 421
 422   }//label >=0
 423   else AliInfo(Form("Negative Kinematic label of PHOS cluster:  %d",label));
 424
 425   return pmother ;
 426
 427 }