PWG4/totEt/AliAnalysisEtMonteCarlo.cxx

   1 //_________________________________________________________________________
   2 //  Utility Class for transverse energy studies
   3 //  Base class for MC analysis
   4 //  - MC output
   5 //  implementation file
   6 //
   7 //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
   8 //_________________________________________________________________________
   9
  10 #include "AliAnalysisEtMonteCarlo.h"
  11 #include "AliAnalysisEtCuts.h"
  12 #include "AliESDtrack.h"
  13 #include "AliStack.h"
  14 #include "AliMCEvent.h"
  15 #include "TH2F.h"
  16 #include "TParticle.h"
  17 #include "AliGenHijingEventHeader.h"
  18 #include "AliGenPythiaEventHeader.h"
  19
  20 using namespace std;
  21
  22 ClassImp(AliAnalysisEtMonteCarlo);
  23
  24
  25 // ctor
  26 AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo() :
  27   AliAnalysisEt()
  28   ,fImpactParameter(0)
  29   ,fNcoll(0)
  30   ,fNpart(0)
  31 {
  32 }
  33
  34 // dtor
  35 AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
  36 {
  37 }
  38
  39 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
  40 { // analyse MC event
  41      ResetEventValues();
  42
  43     // Get us an mc event
  44     AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
  45
  46     Double_t protonMass =fgProtonMass;
  47
  48     // Hijing header
  49     AliGenEventHeader* genHeader = event->GenEventHeader();
  50     AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
  51     if (hijingGenHeader) {
  52       fImpactParameter = hijingGenHeader->ImpactParameter();
  53       fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
  54       fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
  55       /*
  56       printf("Hijing: ImpactParameter %g ReactionPlaneAngle %g \n",
  57              hijingGenHeader->ImpactParameter(), hijingGenHeader->ReactionPlaneAngle());
  58       printf("HardScatters %d ProjecileParticipants %d TargetParticipants %d\n",
  59              hijingGenHeader->HardScatters(), hijingGenHeader->ProjectileParticipants(), hijingGenHeader->TargetParticipants());
  60       printf("ProjSpectatorsn %d ProjSpectatorsp %d TargSpectatorsn %d TargSpectatorsp %d\n",
  61              hijingGenHeader->ProjSpectatorsn(), hijingGenHeader->ProjSpectatorsp(), hijingGenHeader->TargSpectatorsn(), hijingGenHeader->TargSpectatorsp());
  62       printf("NN %d NNw %d NwN %d, NwNw %d\n",
  63              hijingGenHeader->NN(), hijingGenHeader->NNw(), hijingGenHeader->NwN(), hijingGenHeader->NwNw());
  64       */
  65     }
  66
  67     /* // placeholder if we want to get some Pythia info later
  68     AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
  69     if (pythiaGenHeader) { // not Hijing; try with Pythia
  70       printf("Pythia: ProcessType %d  GetPtHard %g \n",
  71              pythiaGenHeader->ProcessType(), pythiaGenHeader->GetPtHard());
  72     }
  73     */
  74
  75     // Let's play with the stack!
  76     AliStack *stack = event->Stack();
  77
  78     Int_t nPrim = stack->GetNtrack();
  79
  80     for (Int_t iPart = 0; iPart < nPrim; iPart++)
  81     {
  82
  83         TParticle *part = stack->Particle(iPart);
  84
  85         if (!part)
  86         {
  87             Printf("ERROR: Could not get particle %d", iPart);
  88             continue;
  89         }
  90
  91         TParticlePDG *pdg = part->GetPDG(0);
  92
  93         Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
  94
  95         // Check if it is a primary particle
  96         if (!stack->IsPhysicalPrimary(iPart)) continue;
  97
  98         // printf("MC: iPart %03d eta %4.3f phi %4.3f code %d charge %g \n", iPart, part->Eta(), part->Phi(), pdg->PdgCode(), pdg->Charge()); // tmp/debug printout
  99
 100         // Check for reasonable (for now neutral and singly charged) charge on the particle
 101         //TODO:Maybe not only singly charged?
 102         if (TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;
 103
 104         fMultiplicity++;
 105
 106         if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
 107         {
 108
 109             if (
 110                 TMath::Abs(pdg->PdgCode()) == fgProtonCode ||
 111                 TMath::Abs(pdg->PdgCode()) == fgNeutronCode ||
 112                 TMath::Abs(pdg->PdgCode()) == fgLambdaCode ||
 113                 TMath::Abs(pdg->PdgCode()) == fgXiCode ||
 114                 TMath::Abs(pdg->PdgCode()) == fgXi0Code ||
 115                 TMath::Abs(pdg->PdgCode()) == fgOmegaCode
 116                )
 117             {
 118               if (pdg->PdgCode() > 0) { particleMassPart = - protonMass;}
 119               if (pdg->PdgCode() < 0) { particleMassPart = protonMass;}
 120             }
 121             Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
 122
 123             if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
 124               {
 125                 fProtonEt += et;
 126               }
 127             if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
 128               {
 129                 fPionEt += et;
 130               }
 131             if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
 132               {
 133                 fChargedKaonEt += et;
 134               }
 135             if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
 136               {
 137                 fMuonEt += et;
 138               }
 139             if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
 140               {
 141                 fElectronEt += et;
 142               }
 143
 144             // some neutrals also
 145             if(pdg->PdgCode() == fgNeutronCode)
 146             {
 147               fNeutronEt += et;
 148             }
 149             if(pdg->PdgCode() == fgAntiNeutronCode)
 150             {
 151               fAntiNeutronEt += et;
 152             }
 153             if(pdg->PdgCode() == fgGammaCode)
 154             {
 155               fGammaEt += et;
 156             }
 157
 158             if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
 159             {
 160                fNeutralMultiplicity++;
 161                fTotNeutralEt += et;
 162
 163                 if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
 164                 {
 165                   fTotNeutralEtAcc += et;
 166                   fTotEtAcc += et;
 167                 }
 168             }
 169             else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())<1e-3 )
 170             {
 171                fChargedMultiplicity++;
 172                fTotChargedEt += et;
 173                 if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
 174                 {
 175                   fTotChargedEtAcc += et;
 176                   fTotEtAcc += et;
 177
 178
 179                     if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
 180                       {
 181                         fProtonEtAcc += et;
 182                       }
 183                     if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
 184                       {
 185                         fPionEtAcc += et;
 186                       }
 187                     if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
 188                       {
 189                         fChargedKaonEtAcc += et;
 190                       }
 191                     if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
 192                       {
 193                         fMuonEtAcc += et;
 194                       }
 195                     if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
 196                       {
 197                         fElectronEtAcc += et;
 198                       }
 199
 200                 }
 201
 202                 //        if (TrackHitsCalorimeter(part, event->GetMagneticField()))
 203                 if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
 204                   {
 205                     if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
 206                     else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
 207                   }
 208             }
 209         }
 210     }
 211
 212     fTotEt = fTotChargedEt + fTotNeutralEt;
 213     fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
 214
 215     FillHistograms();
 216
 217     return 0;
 218 }
 219
 220 void AliAnalysisEtMonteCarlo::Init()
 221 { // init
 222     AliAnalysisEt::Init();
 223 }
 224
 225 void AliAnalysisEtMonteCarlo::ResetEventValues()
 226 { // reset event values
 227   AliAnalysisEt::ResetEventValues();
 228
 229   // collision geometry defaults for p+p:
 230   fImpactParameter = 0;
 231   fNcoll = 1;
 232   fNpart = 2;
 233 }
 234
 235 void AliAnalysisEtMonteCarlo::CreateHistograms()
 236 { // histogram related additions
 237   AliAnalysisEt::CreateHistograms();
 238   if (fTree) {
 239     fTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
 240     fTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
 241     fTree->Branch("fNpart",&fNpart,"fNpart/I");
 242   }
 243 }
 244
 245 bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
 246 {
 247   //  printf(" TrackHitsCalorimeter - magField %f\n", magField);
 248    AliESDtrack *esdTrack = new AliESDtrack(part);
 249    // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 250
 251     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
 252
 253     // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 254
 255     bool status = prop &&
 256       TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
 257       esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
 258       esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
 259     delete esdTrack;
 260
 261     return status;
 262 }
 263