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 "AliVEvent.h"
  15 #include "AliMCEvent.h"
  16 #include "AliESDEvent.h"
  17 #include "TH2F.h"
  18 #include "TParticle.h"
  19 #include "AliGenHijingEventHeader.h"
  20 #include "AliGenPythiaEventHeader.h"
  21 #include "TList.h"
  22 #include "AliESDCaloCluster.h"
  23 #include "AliLog.h"
  24
  25   using namespace std;
  26
  27 ClassImp(AliAnalysisEtMonteCarlo);
  28
  29
  30 // ctor
  31 AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
  32                                                   ,fImpactParameter(0)
  33                                                   ,fNcoll(0)
  34                                                   ,fNpart(0)
  35                                                   ,fHistPrimElectronEtaEET(0)
  36                                                   ,fHistSecElectronEtaEET(0)
  37                                                   ,fHistConvElectronEtaEET(0)
  38                                                   ,fHistPrimGammaEtaEET(0)
  39                                                   ,fHistPion0GammaEtaEET(0)
  40                                                   ,fHistEtaGammaEtaEET(0)
  41                                                   ,fHistOmega0GammaEtaEET(0)
  42                                                   ,fHistSecGammaEtaEET(0)
  43
  44                                                   ,fHistPrimElectronEtaE(0)
  45                                                   ,fHistSecElectronEtaE(0)
  46                                                   ,fHistConvElectronEtaE(0)
  47                                                   ,fHistPrimGammaEtaE(0)
  48                                                   ,fHistPion0GammaEtaE(0)
  49                                                   ,fHistEtaGammaEtaE(0)
  50                                                   ,fHistOmega0GammaEtaE(0)
  51                                                   ,fHistSecGammaEtaE(0)
  52
  53                                                   ,fHistPrimElectronEtaERec(0)
  54                                                   ,fHistSecElectronEtaERec(0)
  55                                                   ,fHistConvElectronEtaERec(0)
  56                                                   ,fHistPrimGammaEtaERec(0)
  57                                                   ,fHistSecGammaEtaERec(0)
  58                                                   ,fHistPion0GammaEtaERec(0)
  59                                                   ,fHistEtaGammaEtaERec(0)
  60                                                   ,fHistOmega0GammaEtaERec(0)
  61
  62                                                   ,fHistAllERecEMC(0)
  63                                                   ,fHistGammaERecEMC(0)
  64                                                   ,fHistElectronERecEMC(0)
  65
  66                                                   ,fHistElectronFirstMother(0)
  67                                                   ,fHistElectronLastMother(0)
  68                                                   ,fHistElectronFirstMotherEtaAcc(0)
  69                                                   ,fHistElectronFirstMotherNPP(0)
  70                                                   ,fHistElectronFirstMotherNPPAcc(0)
  71
  72                                                   ,fHistGammaFirstMother(0)
  73                                                   ,fHistGammaLastMother(0)
  74                                                   ,fHistGammaFirstMotherEtaAcc(0)
  75                                                   ,fHistGammaFirstMotherNPP(0)
  76                                                   ,fHistGammaFirstMotherNPPAcc(0)
  77 {
  78 }
  79
  80 // dtor
  81 AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
  82 {
  83 }
  84
  85 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
  86 { // analyse MC event
  87   ResetEventValues();
  88
  89   // Get us an mc event
  90   if(!ev){
  91     AliFatal("ERROR: Event does not exist");
  92     return 0;
  93   }
  94   AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
  95   if(!event){
  96     AliFatal("ERROR: MC Event does not exist");
  97     return 0;
  98   }
  99
 100   Double_t protonMass =fgProtonMass;
 101
 102   // Hijing header
 103   AliGenEventHeader* genHeader = event->GenEventHeader();
 104   if(!genHeader){
 105     Printf("ERROR: Event generation header does not exist");
 106     return 0;
 107   }
 108   AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
 109   if (hijingGenHeader) {
 110     fImpactParameter = hijingGenHeader->ImpactParameter();
 111     fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
 112     fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
 113     /*
 114       printf("Hijing: ImpactParameter %g ReactionPlaneAngle %g \n",
 115       hijingGenHeader->ImpactParameter(), hijingGenHeader->ReactionPlaneAngle());
 116       printf("HardScatters %d ProjecileParticipants %d TargetParticipants %d\n",
 117       hijingGenHeader->HardScatters(), hijingGenHeader->ProjectileParticipants(), hijingGenHeader->TargetParticipants());
 118       printf("ProjSpectatorsn %d ProjSpectatorsp %d TargSpectatorsn %d TargSpectatorsp %d\n",
 119       hijingGenHeader->ProjSpectatorsn(), hijingGenHeader->ProjSpectatorsp(), hijingGenHeader->TargSpectatorsn(), hijingGenHeader->TargSpectatorsp());
 120       printf("NN %d NNw %d NwN %d, NwNw %d\n",
 121       hijingGenHeader->NN(), hijingGenHeader->NNw(), hijingGenHeader->NwN(), hijingGenHeader->NwNw());
 122     */
 123   }
 124
 125   /* // placeholder if we want to get some Pythia info later
 126      AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 127      if (pythiaGenHeader) { // not Hijing; try with Pythia
 128      printf("Pythia: ProcessType %d  GetPtHard %g \n",
 129      pythiaGenHeader->ProcessType(), pythiaGenHeader->GetPtHard());
 130      }
 131   */
 132
 133   // Let's play with the stack!
 134   AliStack *stack = event->Stack();
 135
 136   Int_t nPrim = stack->GetNtrack();
 137
 138   for (Int_t iPart = 0; iPart < nPrim; iPart++)
 139     {
 140
 141       TParticle *part = stack->Particle(iPart);
 142       TParticle *partMom = 0;
 143       TParticle *partMomLast = 0;
 144
 145       if (!part)
 146         {
 147           Printf("ERROR: Could not get particle %d", iPart);
 148           continue;
 149         }
 150
 151       Int_t iPartMom = part->GetMother(0);
 152       Int_t iPartLastMom = part->GetMother(1);
 153
 154       TParticlePDG *pdg = part->GetPDG(0);
 155       TParticlePDG *pdgMom = 0;
 156       TParticlePDG *pdgMomLast = 0;
 157
 158       if (!pdg)
 159         {
 160           //Printf("ERROR: Could not get particle PDG %d", iPart);
 161           continue;
 162         }
 163
 164       if (iPartMom>0)
 165         {
 166           partMom = stack->Particle(iPartMom);
 167           pdgMom = partMom->GetPDG(0);
 168         }
 169
 170       if (iPartLastMom>0)
 171         {
 172           partMomLast = stack->Particle(iPartLastMom);
 173           pdgMomLast = partMomLast->GetPDG(0);
 174         }
 175
 176
 177       Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
 178
 179       // Check if it is a primary particle
 180       //if (!stack->IsPhysicalPrimary(iPart)) continue;
 181
 182       // it goes to the next particle in case it is not a physical primary or not a electron or gamma (want to check the contribution from secondary)
 183       if (!stack->IsPhysicalPrimary(iPart))
 184         {
 185           // this part is used only to check the origin of secondary electrons and gammas
 186           if (iPartMom>0)
 187             {
 188               //if (!stack->IsPhysicalPrimary(iPartMom)) continue;
 189
 190               if (pdgMom)
 191                 {
 192                   if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
 193                     {
 194                       fHistElectronFirstMotherNPP->Fill(pdgMom->PdgCode());
 195                       // inside EMCal acceptance
 196                       if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
 197                         fHistElectronFirstMotherNPPAcc->Fill(pdgMom->PdgCode());
 198                     }
 199                   else if (pdg->PdgCode() == fgGammaCode)
 200                     {
 201                       fHistGammaFirstMotherNPP->Fill(pdgMom->PdgCode());
 202                       if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
 203                         fHistGammaFirstMotherNPPAcc->Fill(pdgMom->PdgCode());
 204                     }
 205                 }
 206               else
 207                 {
 208                   if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
 209                     fHistElectronFirstMotherNPP->Fill(-598);
 210                   else if (pdg->PdgCode() == fgGammaCode)
 211                     fHistGammaFirstMotherNPP->Fill(-598);
 212
 213                   continue;
 214                 }
 215             }
 216           else
 217             {
 218               if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
 219                 fHistElectronFirstMotherNPP->Fill(-599);
 220               else if (pdg->PdgCode() == fgGammaCode)
 221                 fHistGammaFirstMotherNPP->Fill(-599);
 222
 223               continue;
 224             }
 225
 226           if (!((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode) || (pdg->PdgCode() == fgGammaCode)))
 227             continue;
 228         }
 229       // this part is used only to check the origin of physical primary electrons and gammas
 230       else
 231         {
 232           if (iPartMom>0)
 233             {
 234               if (pdgMom)
 235                 {
 236                   if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
 237                     fHistElectronFirstMother->Fill(pdgMom->PdgCode());
 238                   else if (pdg->PdgCode() == fgGammaCode)
 239                     fHistGammaFirstMother->Fill(pdgMom->PdgCode());
 240
 241                   if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
 242                     {
 243                       if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
 244                         fHistElectronFirstMotherEtaAcc->Fill(pdgMom->PdgCode());
 245                       else if (pdg->PdgCode() == fgGammaCode)
 246                         fHistGammaFirstMotherEtaAcc->Fill(pdgMom->PdgCode());
 247                     }
 248                 }
 249               else
 250                 {
 251                   if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
 252                     fHistElectronFirstMother->Fill(-598);
 253                   else if (pdg->PdgCode() == fgGammaCode)
 254                     fHistGammaFirstMother->Fill(-598);
 255                 }
 256             }
 257           else
 258             {
 259               if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
 260                 fHistElectronFirstMother->Fill(-599);
 261               else if (pdg->PdgCode() == fgGammaCode)
 262                 fHistGammaFirstMother->Fill(-599);
 263             }
 264
 265           if (iPartLastMom>0)
 266             {
 267               if (pdgMomLast)
 268                 {
 269                   if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
 270                     fHistElectronLastMother->Fill(pdgMomLast->PdgCode());
 271                   else if (pdg->PdgCode() == fgGammaCode)
 272                     fHistGammaLastMother->Fill(pdgMomLast->PdgCode());
 273                 }
 274               else
 275                 {
 276                   if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
 277                     fHistElectronLastMother->Fill(-598);
 278                   else if (pdg->PdgCode() == fgGammaCode)
 279                     fHistGammaLastMother->Fill(-598);
 280                 }
 281             }
 282           else
 283             {
 284               if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
 285                 fHistElectronLastMother->Fill(-599);
 286               else if (pdg->PdgCode() == fgGammaCode)
 287                 fHistGammaLastMother->Fill(-599);
 288             }
 289         }
 290
 291       //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
 292
 293       // Check for reasonable (for now neutral and singly charged) charge on the particle
 294       //TODO:Maybe not only singly charged?
 295       if (TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;
 296
 297       fMultiplicity++;
 298
 299       if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
 300         {
 301           // calculate E_T
 302           if (
 303               TMath::Abs(pdg->PdgCode()) == fgProtonCode ||
 304               TMath::Abs(pdg->PdgCode()) == fgNeutronCode ||
 305               TMath::Abs(pdg->PdgCode()) == fgLambdaCode ||
 306               TMath::Abs(pdg->PdgCode()) == fgXiCode ||
 307               TMath::Abs(pdg->PdgCode()) == fgXi0Code ||
 308               TMath::Abs(pdg->PdgCode()) == fgOmegaCode
 309               )
 310             {
 311               if (pdg->PdgCode() > 0) { particleMassPart = - protonMass;}
 312               if (pdg->PdgCode() < 0) { particleMassPart = protonMass;}
 313             }
 314           Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
 315
 316           // Fill up total E_T counters for each particle species
 317           if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
 318             {
 319               fProtonEt += et;
 320             }
 321           if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
 322             {
 323               fPionEt += et;
 324             }
 325           if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
 326             {
 327               fChargedKaonEt += et;
 328             }
 329           if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
 330             {
 331               fMuonEt += et;
 332             }
 333           if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
 334             {
 335               fElectronEt += et;
 336             }
 337
 338           // some neutrals also
 339           if(pdg->PdgCode() == fgNeutronCode)
 340             {
 341               fNeutronEt += et;
 342             }
 343           if(pdg->PdgCode() == fgAntiNeutronCode)
 344             {
 345               fAntiNeutronEt += et;
 346             }
 347           if(pdg->PdgCode() == fgGammaCode)
 348             {
 349               fGammaEt += et;
 350             }
 351
 352           // Neutral particles
 353           if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
 354             {
 355               fNeutralMultiplicity++;
 356               fTotNeutralEt += et;
 357
 358               // inside EMCal acceptance
 359               if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
 360                 {
 361                   fTotNeutralEtAcc += et;
 362                   fTotEtAcc += et;
 363
 364                   if(pdg->PdgCode() == fgGammaCode)
 365                     {
 366                       if (!stack->IsPhysicalPrimary(iPart))
 367                         {
 368                           fHistSecGammaEtaE->Fill(part->Energy(),part->Eta());
 369                           fHistSecGammaEtaEET->Fill(part->Energy(),part->Eta(),et);
 370                         }
 371                       else
 372                         {
 373                           if (pdgMom)
 374                             {
 375                               if (pdgMom->PdgCode() == fgPi0Code)
 376                                 {
 377                                   fHistPion0GammaEtaE->Fill(part->Energy(),part->Eta());
 378                                   fHistPion0GammaEtaEET->Fill(part->Energy(),part->Eta(),et);
 379                                 }
 380                               else if (pdgMom->PdgCode() == fgOmega0Code)
 381                                 {
 382                                   fHistOmega0GammaEtaE->Fill(part->Energy(),part->Eta());
 383                                   fHistOmega0GammaEtaEET->Fill(part->Energy(),part->Eta(),et);
 384                                 }
 385                               else if (pdgMom->PdgCode() == fgEtaCode)
 386                                 {
 387                                   fHistEtaGammaEtaE->Fill(part->Energy(),part->Eta());
 388                                   fHistEtaGammaEtaEET->Fill(part->Energy(),part->Eta(),et);
 389                                 }
 390                               else
 391                                 {
 392                                   fHistPrimGammaEtaE->Fill(part->Energy(),part->Eta());
 393                                   fHistPrimGammaEtaEET->Fill(part->Energy(),part->Eta(),et);
 394                                 }
 395                             }
 396                           else
 397                             {
 398                               fHistPrimGammaEtaE->Fill(part->Energy(),part->Eta());
 399                               fHistPrimGammaEtaEET->Fill(part->Energy(),part->Eta(),et);
 400                             }
 401                         }
 402                     }
 403                 }
 404             }
 405           //Charged particles
 406           else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
 407             {
 408               fChargedMultiplicity++;
 409               fTotChargedEt += et;
 410
 411               // inside EMCal acceptance
 412               if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
 413                 {
 414                   fTotChargedEtAcc += et;
 415                   fTotEtAcc += et;
 416
 417                   if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
 418                     {
 419                       fProtonEtAcc += et;
 420                     }
 421                   if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
 422                     {
 423                       fPionEtAcc += et;
 424                     }
 425                   if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
 426                     {
 427                       fChargedKaonEtAcc += et;
 428                     }
 429                   if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
 430                     {
 431                       fMuonEtAcc += et;
 432                     }
 433                   /*
 434                     if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
 435                     {
 436                     fElectronEtAcc += et;
 437                     }
 438                   */
 439                   if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
 440                     {
 441                       if (!stack->IsPhysicalPrimary(iPart))
 442                         {
 443                           if (pdgMom)
 444                             {
 445                               if ((pdgMom->PdgCode() == fgGammaCode) && (stack->IsPhysicalPrimary(iPartMom)))
 446                                 {
 447                                   fHistConvElectronEtaE->Fill(part->Energy(),part->Eta());
 448                                   fHistConvElectronEtaEET->Fill(part->Energy(),part->Eta(),et);
 449                                 }
 450                               else
 451                                 {
 452                                   fHistSecElectronEtaE->Fill(part->Energy(),part->Eta());
 453                                   fHistSecElectronEtaEET->Fill(part->Energy(),part->Eta(),et);
 454                                 }
 455                             }
 456                           else
 457                             {
 458                               fHistSecElectronEtaE->Fill(part->Energy(),part->Eta());
 459                               fHistSecElectronEtaEET->Fill(part->Energy(),part->Eta(),et);
 460                             }
 461                         }
 462                       else
 463                         {
 464                           fElectronEtAcc += et;
 465                           fHistPrimElectronEtaE->Fill(part->Energy(),part->Eta());
 466                           fHistPrimElectronEtaEET->Fill(part->Energy(),part->Eta(),et);
 467                         }
 468                     } // electron
 469                 } // inside EMCal acceptance
 470
 471               //          if (TrackHitsCalorimeter(part, event->GetMagneticField()))
 472               if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
 473                 {
 474                   if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
 475                   else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
 476                 }
 477             }
 478         }
 479     }
 480
 481   fTotEt = fTotChargedEt + fTotNeutralEt;
 482   fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
 483
 484   FillHistograms();
 485
 486   return 0;
 487 }
 488 //Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliMCEvent* mcEvent,AliESDEvent* realEvent)
 489 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
 490 { // analyse MC and real event info
 491   //if(!mcEvent || !realEvent){
 492   if(!ev || !ev2){
 493     AliFatal("ERROR: Event does not exist");
 494     return 0;
 495   }
 496   //AnalyseEvent(mcEvent);
 497   AnalyseEvent(ev);
 498   AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
 499   AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
 500   if(!mcEvent || !realEvent){
 501     AliFatal("ERROR: mcEvent or realEvent does not exist");
 502     return 0;
 503   }
 504
 505   AliStack *stack = mcEvent->Stack();
 506
 507   // get all emcal clusters
 508   TRefArray* caloClusters = new TRefArray();
 509   realEvent->GetEMCALClusters( caloClusters );
 510
 511   Int_t nCluster = caloClusters->GetEntries();
 512
 513   // loop the clusters
 514   for (int iCluster = 0; iCluster < nCluster; iCluster++ )
 515     {
 516       AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
 517       Float_t caloE = caloCluster->E();
 518
 519       UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
 520       TParticle *part  = stack->Particle(iPart);
 521       TParticle *partMom = 0;
 522
 523       if (!part)
 524         {
 525           Printf("No MC particle %d", iCluster);
 526           continue;
 527         }
 528
 529       // compare MC and Rec energies for all particles
 530       fHistAllERecEMC->Fill(part->Energy(),caloE);
 531
 532       TParticlePDG *pdg = part->GetPDG(0);
 533       TParticlePDG *pdgMom = 0;
 534
 535       if (!pdg)
 536         {
 537           Printf("ERROR: Could not get particle PDG %d", iPart);
 538           continue;
 539         }
 540
 541       Int_t iPartMom = part->GetMother(0);
 542
 543       if (iPartMom>0)
 544         {
 545           partMom = stack->Particle(iPartMom);
 546           pdgMom = partMom->GetPDG(0);
 547         }
 548
 549       // Check if it is a primary particle
 550       //if (!stack->IsPhysicalPrimary(iPart)) continue;
 551       if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
 552         {
 553           if (!((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode) || (pdg->PdgCode() == fgGammaCode)))
 554             continue;
 555         } // end of primary particle check
 556
 557       if (TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;
 558
 559       if(pdg->PdgCode() == fgGammaCode)
 560         {
 561           // compare MC and Rec energies for gammas
 562           fHistGammaERecEMC->Fill(part->Energy(),caloE);
 563
 564           if (!stack->IsPhysicalPrimary(iPart))
 565             {
 566               fHistSecGammaEtaERec->Fill(part->Energy(),part->Eta());
 567             }
 568           else
 569             {
 570               if (pdgMom)
 571                 {
 572                   if (pdgMom->PdgCode() == fgPi0Code)
 573                     {
 574                       fHistPion0GammaEtaERec->Fill(part->Energy(),part->Eta());
 575                     }
 576                   else if (partMom->GetPDG(0)->PdgCode() == fgOmega0Code)
 577                     {
 578                       fHistOmega0GammaEtaERec->Fill(part->Energy(),part->Eta());
 579                     }
 580                   else if (partMom->GetPDG(0)->PdgCode() == fgEtaCode)
 581                     {
 582                       fHistEtaGammaEtaERec->Fill(part->Energy(),part->Eta());
 583                     }
 584                   else
 585                     {
 586                       fHistPrimGammaEtaERec->Fill(part->Energy(),part->Eta());
 587                     }
 588                 }
 589               else
 590                 {
 591                   fHistPrimGammaEtaERec->Fill(part->Energy(),part->Eta());
 592                 }
 593             }
 594         } // gamma
 595
 596       if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
 597         {
 598           // compare MC and Rec energies for electrons
 599           fHistElectronERecEMC->Fill(part->Energy(),caloE);
 600
 601           if (!stack->IsPhysicalPrimary(iPart))
 602             {
 603               if (pdgMom)
 604                 {
 605                   if ((pdgMom->PdgCode() == fgGammaCode) && (stack->IsPhysicalPrimary(iPartMom)))
 606                     {
 607                       fHistConvElectronEtaERec->Fill(part->Energy(),part->Eta());
 608                     }
 609                   else
 610                     {
 611                       fHistSecElectronEtaERec->Fill(part->Energy(),part->Eta());
 612                     }
 613                 }
 614               else
 615                 {
 616                   fHistSecElectronEtaERec->Fill(part->Energy(),part->Eta());
 617                 }
 618             }
 619           else
 620             {
 621               fHistPrimElectronEtaERec->Fill(part->Energy(),part->Eta());
 622             }
 623         }
 624     } // end of loop over clusters
 625   return 0;
 626 }
 627
 628 void AliAnalysisEtMonteCarlo::Init()
 629 { // init
 630   AliAnalysisEt::Init();
 631 }
 632
 633 void AliAnalysisEtMonteCarlo::ResetEventValues()
 634 { // reset event values
 635   AliAnalysisEt::ResetEventValues();
 636
 637   // collision geometry defaults for p+p:
 638   fImpactParameter = 0;
 639   fNcoll = 1;
 640   fNpart = 2;
 641 }
 642
 643 void AliAnalysisEtMonteCarlo::CreateHistograms()
 644 { // histogram related additions
 645   AliAnalysisEt::CreateHistograms();
 646   if (fTree) {
 647     fTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
 648     fTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
 649     fTree->Branch("fNpart",&fNpart,"fNpart/I");
 650   }
 651   fHistPrimElectronEtaEET = CreateEtaEHisto2D("fHistPrimElectronEtaEET","MC E_{T}, primary electrons","E_{T}(GeV)");
 652   fHistSecElectronEtaEET = CreateEtaEHisto2D("fHistSecElectronEtaEET","MC E_{T}, secondary (no conversion) electrons","E_{T}(GeV)");
 653   fHistConvElectronEtaEET = CreateEtaEHisto2D("fHistConvElectronEtaEET","MC E_{T}, electrons from conversion","E_{T}(GeV)");
 654   fHistPrimGammaEtaEET = CreateEtaEHisto2D("fHistPrimGammaEtaEET","MC E_{T}, primary gammas","E_{T}(GeV)");
 655   fHistPion0GammaEtaEET = CreateEtaEHisto2D("fHistPion0GammaEtaEET","MC E_{T}, #pi^{0}","E_{T}(GeV)");
 656   fHistOmega0GammaEtaEET = CreateEtaEHisto2D("fHistOmega0GammaEtaEET","MC E_{T}, #omega^{0}","E_{T}(GeV)");
 657   fHistEtaGammaEtaEET = CreateEtaEHisto2D("fHistEtaGammaEtaEET","MC E_{T}, #eta","E_{T}(GeV)");
 658   fHistSecGammaEtaEET = CreateEtaEHisto2D("fHistSecGammaEtaEET","MC E_{T}, secondary (no #pi^{0}, #eta or #omega) gammas","E_{T}(GeV)");
 659
 660   fHistPrimElectronEtaE = CreateEtaEHisto2D("fHistPrimElectronEtaE","MC E_{T}, primary electrons","#");
 661   fHistSecElectronEtaE = CreateEtaEHisto2D("fHistSecElectronEtaE","MC E_{T}, secondary (no conversion) electrons","#");
 662   fHistConvElectronEtaE = CreateEtaEHisto2D("fHistConvElectronEtaE","MC E_{T}, electrons from conversion","#");
 663   fHistPrimGammaEtaE = CreateEtaEHisto2D("fHistPrimGammaEtaE","MC E_{T}, primary gammas","#");
 664   fHistPion0GammaEtaE = CreateEtaEHisto2D("fHistPion0GammaEtaE","MC E_{T}, #pi^{0}","#");
 665   fHistOmega0GammaEtaE = CreateEtaEHisto2D("fHistOmega0GammaEtaE","MC E_{T}, #omega^{0}","#");
 666   fHistEtaGammaEtaE = CreateEtaEHisto2D("fHistEtaGammaEtaE","MC E_{T}, #eta","#");
 667   fHistSecGammaEtaE = CreateEtaEHisto2D("fHistSecGammaEtaE","MC E_{T}, secondary (no #pi^{0}, #eta or #omega) gammas","#");
 668
 669   fHistPrimElectronEtaERec = CreateEtaEHisto2D("fHistPrimElectronEtaERec","MC E_{T}, primary electrons","#");
 670   fHistSecElectronEtaERec = CreateEtaEHisto2D("fHistSecElectronEtaERec","MC E_{T}, secondary (no conversion) electrons","#");
 671   fHistConvElectronEtaERec = CreateEtaEHisto2D("fHistConvElectronEtaERec","MC E_{T}, electrons from conversion","#");
 672   fHistPrimGammaEtaERec = CreateEtaEHisto2D("fHistPrimGammaEtaERec","MC E_{T}, primary gammas","#");
 673   fHistPion0GammaEtaERec = CreateEtaEHisto2D("fHistPion0GammaEtaERec","MC E_{T}, #pi^{0}","#");
 674   fHistOmega0GammaEtaERec = CreateEtaEHisto2D("fHistOmega0GammaEtaERec","MC E_{T}, #omega","#");
 675   fHistEtaGammaEtaERec = CreateEtaEHisto2D("fHistEtaGammaEtaERec","MC E_{T}, #eta","#");
 676   fHistSecGammaEtaERec = CreateEtaEHisto2D("fHistSecGammaEtaERec","MC E_{T}, secondary (no #pi^{0}, #eta or #omega) gammas","#");
 677
 678   fHistAllERecEMC = new TH2F("fHistAllERecEMC","E cluster Rec vs MC, all particles",fgNumOfEBins, fgEAxis,fgNumOfEBins, fgEAxis);
 679   fHistAllERecEMC->SetXTitle("E_{MC}(GeV)");
 680   fHistAllERecEMC->SetYTitle("E_{Rec}(GeV)");
 681
 682   fHistElectronERecEMC = new TH2F("fHistElectronERecEMC","E cluster Rec vs MC, Electrons",fgNumOfEBins, fgEAxis,fgNumOfEBins, fgEAxis);
 683   fHistElectronERecEMC->SetXTitle("E_{MC}(GeV)");
 684   fHistElectronERecEMC->SetYTitle("E_{Rec}(GeV)");
 685
 686   fHistGammaERecEMC = new TH2F("fHistGammaERecEMC","E cluster Rec vs MC, Gammas",fgNumOfEBins, fgEAxis,fgNumOfEBins, fgEAxis);
 687   fHistGammaERecEMC->SetXTitle("E_{MC}(GeV)");
 688   fHistGammaERecEMC->SetYTitle("E_{Rec}(GeV)");
 689
 690   fHistElectronFirstMother = new TH1F("fHistElectronFirstMother","Electron First Mother PDG Code Distribution",1201,-600.5,600.5);
 691   fHistElectronLastMother = new TH1F("fHistElectronLastMother","Electron Last Mother PDG Code Distribution",1201,-600.5,600.5);
 692   fHistElectronFirstMotherEtaAcc = new TH1F("fHistElectronFirstMotherEtaAcc","Electron First Mother PDG Code Distribution",1201,-600.5,600.5);
 693   fHistElectronFirstMotherNPP = new TH1F("fHistElectronFirstMotherNPP","Electron First Mother PDG Code Distribution",1201,-600.5,600.5);
 694   fHistElectronFirstMotherNPPAcc = new TH1F("fHistElectronFirstMotherNPPAcc","Electron First Mother PDG Code Distribution",1201,-600.5,600.5);
 695
 696   fHistGammaFirstMother = new TH1F("fHistGammaFirstMother","Gamma First Mother PDG Code Distribution",1201,-600.5,600.5);
 697   fHistGammaLastMother = new TH1F("fHistGammaLastMother","Gamma Last Mother PDG Code Distribution",1201,-600.5,600.5);
 698   fHistGammaFirstMotherEtaAcc = new TH1F("fHistGammaFirstMotherEtaAcc","Gamma First Mother PDG Code Distribution",1201,-600.5,600.5);
 699   fHistGammaFirstMotherNPP = new TH1F("fHistGammaFirstMotherNPP","Gamma First Mother PDG Code Distribution",1201,-600.5,600.5);
 700   fHistGammaFirstMotherNPPAcc = new TH1F("fHistGammaFirstMotherNPP","Gamma First Mother PDG Code Distribution",1201,-600.5,600.5);
 701 }
 702
 703 void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
 704 {//fill the output list
 705   AliAnalysisEt::FillOutputList(list);
 706
 707   list->Add(fHistPrimElectronEtaEET);
 708   list->Add(fHistSecElectronEtaEET);
 709   list->Add(fHistConvElectronEtaEET);
 710   list->Add(fHistPrimGammaEtaEET);
 711   list->Add(fHistPion0GammaEtaEET);
 712   list->Add(fHistOmega0GammaEtaEET);
 713   list->Add(fHistEtaGammaEtaEET);
 714   list->Add(fHistSecGammaEtaEET);
 715
 716   list->Add(fHistPrimElectronEtaE);
 717   list->Add(fHistSecElectronEtaE);
 718   list->Add(fHistConvElectronEtaE);
 719   list->Add(fHistPrimGammaEtaE);
 720   list->Add(fHistPion0GammaEtaE);
 721   list->Add(fHistOmega0GammaEtaE);
 722   list->Add(fHistEtaGammaEtaE);
 723   list->Add(fHistSecGammaEtaE);
 724
 725   list->Add(fHistPrimElectronEtaERec);
 726   list->Add(fHistSecElectronEtaERec);
 727   list->Add(fHistConvElectronEtaERec);
 728   list->Add(fHistPrimGammaEtaERec);
 729   list->Add(fHistPion0GammaEtaERec);
 730   list->Add(fHistOmega0GammaEtaERec);
 731   list->Add(fHistEtaGammaEtaERec);
 732   list->Add(fHistSecGammaEtaERec);
 733
 734   list->Add(fHistAllERecEMC);
 735   list->Add(fHistElectronERecEMC);
 736   list->Add(fHistGammaERecEMC);
 737
 738   list->Add(fHistElectronFirstMother);
 739   list->Add(fHistElectronLastMother);
 740   list->Add(fHistElectronFirstMotherEtaAcc);
 741   list->Add(fHistElectronFirstMotherNPP);
 742   list->Add(fHistElectronFirstMotherNPPAcc);
 743
 744   list->Add(fHistGammaFirstMother);
 745   list->Add(fHistGammaLastMother);
 746   list->Add(fHistGammaFirstMotherEtaAcc);
 747   list->Add(fHistGammaFirstMotherNPP);
 748   list->Add(fHistGammaFirstMotherNPPAcc);
 749 }
 750
 751
 752 bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
 753 {
 754   //  printf(" TrackHitsCalorimeter - magField %f\n", magField);
 755   AliESDtrack *esdTrack = new AliESDtrack(part);
 756   // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 757
 758   Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
 759
 760   // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 761
 762   bool status = prop &&
 763     TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
 764     esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
 765     esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
 766   delete esdTrack;
 767
 768   return status;
 769 }
 770