EMCAL/AliEMCALJetFinderInputSimPrep.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 /* $Id$ */
  18
  19 //_________________________________________________________________________
  20 //  Class for JetFinder Input preparation from simulated data
  21 //
  22 //*-- Author: Mark Horner (LBL/UCT)
  23 //
  24
  25
  26
  27 #include "AliEMCALJetFinderInputSimPrep.h"
  28
  29 #include <TParticle.h>
  30 #include <TParticlePDG.h>
  31 #include <TPDGCode.h>
  32 #include <TFile.h>
  33 #include <TTree.h>
  34 #include <TObjectTable.h>
  35 #include <TMath.h>
  36
  37
  38 #include "AliRun.h"
  39 #include "AliRunLoader.h"
  40 #include "AliStack.h"
  41 #include "AliMagF.h"
  42 #include "AliEMCALFast.h"
  43 #include "AliEMCAL.h"
  44 #include "AliEMCALHit.h"
  45 #include "AliEMCALGeometry.h"
  46 #include "AliEMCALLoader.h"
  47 #include "AliGenEventHeader.h"
  48 #include "AliGenPythiaEventHeader.h"
  49 #include "AliGenerator.h"
  50 #include "AliHeader.h"
  51 #include "AliMC.h"
  52
  53
  54
  55 ClassImp(AliEMCALJetFinderInputSimPrep)
  56
  57 AliEMCALJetFinderInputSimPrep::AliEMCALJetFinderInputSimPrep()
  58 {
  59         // Default constructor
  60 if (fDebug > 0) Info("AliEMCALJetFinderInputSimPrep","Beginning Constructor");
  61
  62   fDebug = 0;
  63   fSmearType = kSmearEffic;
  64   fEMCALType = kHits;
  65   fTrackType = kCharged;
  66   fEfficiency = 0.90;
  67
  68 }
  69 AliEMCALJetFinderInputSimPrep::~AliEMCALJetFinderInputSimPrep()
  70 {
  71 if (fDebug > 0) Info("~AliEMCALJetFinderInputSimPrep","Beginning Destructor");
  72
  73 }
  74
  75 void AliEMCALJetFinderInputSimPrep::Reset(AliEMCALJetFinderResetType_t resettype)
  76 {
  77         // Method to reset data
  78   if (fDebug > 1) Info("Reset","Beginning Reset");
  79         switch (resettype){
  80
  81         case kResetData:
  82                 fInputObject.Reset(resettype);
  83                 break;
  84         case kResetTracks:
  85                 fInputObject.Reset(resettype);
  86                 break;
  87         case kResetDigits:
  88                 fInputObject.Reset(resettype);
  89                 break;
  90         case kResetParameters:
  91                 fSmearType = kSmearEffic;
  92                 fEMCALType = kHits;
  93                 fTrackType = kCharged;
  94                 break;
  95         case kResetAll:
  96                 fSmearType = kSmearEffic;
  97                 fEMCALType = kHits;
  98                 fTrackType = kCharged;
  99                 fInputObject.Reset(resettype);
 100                 break;
 101         case kResetPartons:
 102           Warning("FillFromFile", "kResetPartons not implemented") ;
 103           break;
 104         case kResetParticles:
 105           Warning("FillFromFile", "kResetParticles not implemented") ;
 106           break;
 107         case kResetJets:
 108           Warning("FillFromFile", "kResetJets not implemented") ;
 109           break;
 110         }// end switch
 111
 112 }
 113
 114 Int_t AliEMCALJetFinderInputSimPrep::FillFromFile(TString *filename, AliEMCALJetFinderFileType_t filetype,Int_t EventNumber,TString data)
 115 {
 116   if (fDebug > 1) Info("FillFromFile","Beginning FillFromFile");
 117   fFileType = filetype;
 118   AliRunLoader *rl=AliRunLoader::Open(filename->Data());
 119   if (fDebug > 1) Info("FillFromFile","Opened file %s",filename->Data());
 120    if (data.Contains("S"))
 121    {
 122            rl->LoadKinematics();
 123            rl->LoadSDigits();
 124            rl->GetEvent(EventNumber);
 125    }else
 126    {
 127            rl->LoadKinematics();
 128            rl->LoadHits();
 129            rl->GetEvent(EventNumber);
 130    }
 131    if (fDebug > 1) Info("FillFromFile","Got event %i with option \"XH\"",EventNumber);
 132
 133    if ( fEMCALType == kHits ||
 134         fEMCALType == kTimeCut )
 135    {
 136            if (data.Contains("S"))
 137            {
 138                    FillSDigits();
 139            }else
 140            {
 141                    FillHits();
 142            }
 143    }
 144    if ( fTrackType != kNoTracks  )  FillTracks();
 145    if ( fFileType  != kData){
 146            FillPartons();
 147    }
 148    return 0;
 149 }
 150
 151 void AliEMCALJetFinderInputSimPrep::FillHits()
 152 {
 153         // Fill from the hits to input object from simulation
 154 if (fDebug > 1) Info("FillHits","Beginning FillHits");
 155
 156 // Access hit information
 157 //  AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 158 //  Info("AliEMCALJetFinderInputSimPrep","Title of the geometry is %s",pEMCAL->GetTitle());
 159     AliEMCALGeometry* geom =  AliEMCALGeometry::GetInstance();
 160     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 161
 162     //    TTree *treeH = AliEMCALGetter::Instance()->TreeH();
 163 //    Int_t ntracks = (Int_t) treeH->GetEntries();
 164 //
 165 //   Loop over tracks
 166 //
 167     Double_t etH = 0.0;
 168
 169 /*    for (Int_t track=0; track<ntracks;track++) {
 170         gAlice->ResetHits();
 171         nbytes += treeH->GetEvent(track);*/
 172 //
 173 //  Loop over hits
 174 //
 175         for(Int_t i =0; i < emcalLoader->Hits()->GetEntries() ;i++)
 176         {
 177             const AliEMCALHit *mHit = emcalLoader->Hit(i);
 178             if (fEMCALType == kTimeCut &&
 179                     (mHit->GetTime()>fTimeCut) ) continue;
 180             Float_t x      =    mHit->X();         // x-pos of hit
 181             Float_t y      =    mHit->Y();         // y-pos
 182             Float_t z      =    mHit->Z();         // z-pos
 183             Float_t eloss  =    mHit->GetEnergy(); // deposited energy
 184
 185             Float_t r      =    TMath::Sqrt(x*x+y*y);
 186             Float_t theta  =    TMath::ATan2(r,z);
 187             Float_t eta    =   -TMath::Log(TMath::Tan(theta/2.));
 188             Float_t phi    =    TMath::ATan2(y,x);
 189             etH = eloss*TMath::Sin(theta);
 190             if (fDebug > 10) Info("FillHits","Values of hit energy %i",Int_t(1e7*etH));
 191             if (geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi) == -1)
 192             {
 193                     if (fDebug >5)
 194                     {
 195                             Error("FillHits","Hit not inside EMCAL!");
 196                             mHit->Dump();
 197                     }
 198                     continue;
 199             }
 200             fInputObject.AddEnergyToDigit(geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi)-1,Int_t(1e7*etH));
 201         } // Hit Loop
 202 //    } // Track Loop
 203
 204
 205 }
 206 void AliEMCALJetFinderInputSimPrep::FillTracks()
 207 {
 208         // Fill from particles simulating a TPC to input object from simulation
 209     if (fDebug > 1) Info("FillTracks","Beginning FillTracks");
 210
 211     AliRunLoader *rl = AliRunLoader::GetRunLoader();
 212     TParticlePDG* pdgP = 0;
 213     TParticle *mPart;
 214     Int_t npart = rl->Stack()->GetNprimary();
 215     if (fDebug > 1) Info("FillTracks","Header says there are %i primaries",npart);
 216     Float_t bfield,rEMCAL;
 217
 218     if (fDebug > 1) Info("FillTracks","Defining the geometry");
 219
 220     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 221     AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
 222     if (geom == 0 && pEMCAL)
 223       geom = AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
 224     if (geom == 0)
 225       geom = AliEMCALGeometry::GetInstance("EMCAL_55_25","");//","SHISH_77_TRD1_2X2_FINAL_110DEG");
 226     fEtaMax = geom->GetArm1EtaMax();
 227     fEtaMin = geom->GetArm1EtaMin();
 228     fPhiMax = TMath::Pi()/180.0*geom->GetArm1PhiMax();
 229     fPhiMin = TMath::Pi()/180.0*geom->GetArm1PhiMin();
 230
 231     if (fDebug > 1) Info("FillTracks","Starting particle loop");
 232
 233     if (fDebug > 1) Info("FillTracks","Particle loop of %i",npart);
 234     for (Int_t part = 0; part < npart; part++) {
 235         mPart = rl->Stack()->Particle(part);
 236         pdgP = mPart->GetPDG();
 237
 238         if (fDebug > 5) Info("FillTracks","Checking if track is a primary");
 239
 240         if (fFileType == kPythia) {
 241             if (mPart->GetStatusCode() != 1) continue;
 242         } else if (fFileType == kHijing) {
 243             if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
 244         }
 245
 246         if (fDebug > 15) Info("FillTracks","Checking if track (eta - %f, phi - %f) is in acceptance",mPart->Eta(),mPart->Phi());
 247         if (fDebug > 10) Info("FillTracks","Checking if EMCAL acceptance  ( %f < eta < %f, %f < phi < %f) is in acceptance",fEtaMin,fEtaMax,fPhiMin,fPhiMax);
 248
 249         if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin))   continue;
 250         if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin))   continue;
 251
 252 /*
 253         {kProton, kProtonBar, kElectron, kPositron,
 254          kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
 255          kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
 256          kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
 257          kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
 258          0,kNuMu,kNuMuBar
 259 */
 260
 261         if (fDebug > 15) Info("FillTracks","Checking if track is rejected");
 262
 263         if ((fSmearType == kEfficiency  ||
 264              fSmearType == kSmearEffic) &&
 265              pdgP->Charge()!=0) {
 266             if (AliEMCALFast::RandomReject(fEfficiency)) {
 267                 continue;
 268             }
 269         }
 270
 271         if (gAlice && gAlice->Field())
 272           bfield = gAlice->Field()->SolenoidField();
 273         else
 274           bfield = 0.4;
 275         rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
 276         Float_t rB = 3335.6 * mPart->Pt() / bfield;  // [cm]  (case of |charge|=1)
 277         if (2.*rB < rEMCAL) continue;  // track curls away
 278
 279         //if (part%10) gObjectTable->Print();
 280         switch(fTrackType)
 281         {
 282
 283            case kAllP:  // All Stable particles to be included
 284                 if (fDebug > 5) Info("FillTracks","Storing track");
 285                 if (fSmearType == kSmear ||
 286                     fSmearType == kSmearEffic ){
 287                         Smear(mPart);/*
 288                         TParticle *tmp = Smear(MPart);
 289                         fInputObject.AddTrack(Smear(MPart));
 290                         delete tmp;*/
 291                 }else{
 292                         fInputObject.AddTrack(*mPart);
 293                 }
 294            break;
 295            case kEM:   // All Electromagnetic particles
 296                 if (mPart->GetPdgCode() == kElectron  ||
 297                     mPart->GetPdgCode() == kMuonPlus  ||
 298                     mPart->GetPdgCode() == kMuonMinus ||
 299                     mPart->GetPdgCode() == kPositron ){
 300                       if (fDebug > 5) Info("FillTracks","Storing electron or positron");
 301                       if (fSmearType == kSmear ||
 302                             fSmearType == kSmearEffic ){
 303                               Smear(mPart);/*
 304                               TParticle *tmp = Smear(MPart);
 305                               fInputObject.AddTrack(tmp);
 306                               delete tmp;*/
 307                       }else{
 308                           fInputObject.AddTrack(*mPart);
 309                       }
 310                 }
 311                 if ( mPart->GetPdgCode() == kGamma ){
 312                         fInputObject.AddTrack(*mPart);
 313                         if (fDebug > 5) Info("FillTracks","Storing gamma");
 314                 }
 315
 316            break;
 317            case kCharged: // All Particles with non-zero charge
 318                 if (pdgP->Charge() != 0) {
 319                         if (fDebug > 5) Info("FillTracks","Storing charged track");
 320                         if (fSmearType == kSmear ||
 321                             fSmearType == kSmearEffic ){
 322                                 Smear(mPart);/*
 323                                 TParticle *tmp = Smear(MPart);
 324                                 fInputObject.AddTrack(tmp);
 325                                 delete tmp;*/
 326                         }else{
 327                          fInputObject.AddTrack(*mPart);
 328                         }
 329                 }
 330            break;
 331            case kNeutral: // All particles with no charge
 332                 if (pdgP->Charge() == 0){
 333                         fInputObject.AddTrack(*mPart);
 334                         if (fDebug > 5) Info("FillTracks","Storing neutral");
 335                 }
 336            break;
 337            case kHadron: //All hadrons
 338                 if (mPart->GetPdgCode() != kElectron  &&
 339                     mPart->GetPdgCode() != kPositron  &&
 340                     mPart->GetPdgCode() != kMuonPlus  &&
 341                     mPart->GetPdgCode() != kMuonMinus &&
 342                     mPart->GetPdgCode() != kGamma )
 343                 {
 344                         if (fDebug > 5) Info("FillTracks","Storing hadron");
 345                         if (pdgP->Charge() == 0){
 346                             fInputObject.AddTrack(*mPart);
 347                         }else{
 348                                 if (fSmearType == kSmear ||
 349                                     fSmearType == kSmearEffic ){
 350                                         Smear(mPart);/*
 351                                         TParticle *tmp = Smear(MPart);
 352                                         fInputObject.AddTrack(tmp);
 353                                         delete tmp;*/
 354                                 }else{
 355                                     fInputObject.AddTrack(*mPart);
 356                                 }
 357                         }
 358                 }
 359            break;
 360            case kChargedHadron:  // only charged hadrons
 361                 if (mPart->GetPdgCode() != kElectron  &&
 362                     mPart->GetPdgCode() != kPositron  &&
 363                     mPart->GetPdgCode() != kGamma     &&
 364                     mPart->GetPdgCode() != kMuonPlus  &&
 365                     mPart->GetPdgCode() != kMuonMinus &&
 366                     pdgP->Charge()      != 0       ){
 367                         if (fDebug > 5) Info("FillTracks","Storing charged hadron");
 368                         if (fSmearType == kSmear ||
 369                             fSmearType == kSmearEffic ){
 370                                 Smear(mPart);/*
 371                                 TParticle *tmp = Smear(MPart);
 372                                 fInputObject.AddTrack(tmp);
 373                                 delete tmp;*/
 374                         }else{
 375                                fInputObject.AddTrack(*mPart);
 376                         }
 377                 }
 378            break;
 379            case kNoTracks:
 380            break;
 381            case kEMChargedPi0:
 382                 if (pdgP->Charge() != 0 || mPart->GetPdgCode() == kPi0  ||
 383                         mPart->GetPdgCode() == kGamma     )
 384                 {
 385                         if (fDebug > 5) Info("FillTracks","Storing charged track");
 386                         if (fSmearType == kSmear ||
 387                             fSmearType == kSmearEffic ){
 388                                 Smear(mPart);/*
 389                                 TParticle *tmp = Smear(MPart);
 390                                 fInputObject.AddTrack(tmp);
 391                                 delete tmp;*/
 392                         }else{
 393                          fInputObject.AddTrack(*mPart);
 394                         }
 395                 }
 396            break;
 397            case kNoNeutronNeutrinoKlong:
 398                 if ( mPart->GetPdgCode() != kNeutron    &&
 399                      mPart->GetPdgCode() != kNeutronBar &&
 400                      mPart->GetPdgCode() != kK0Long     &&
 401                      mPart->GetPdgCode() != kNuE        &&
 402                      mPart->GetPdgCode() != kNuEBar     &&
 403                      mPart->GetPdgCode() != kNuMu       &&
 404                      mPart->GetPdgCode() != kNuMuBar    &&
 405                      mPart->GetPdgCode() != kNuTau      &&
 406                      mPart->GetPdgCode() != kNuTauBar   )
 407                 {
 408                         if (fDebug > 5) Info("FillTracks","Storing charged track");
 409                         if (fSmearType == kSmear ||
 410                                         fSmearType == kSmearEffic ){
 411                                 Smear(mPart);/*
 412                                                 TParticle *tmp = Smear(MPart);
 413                                                 fInputObject.AddTrack(tmp);
 414                                                 delete tmp;*/
 415                         }else{
 416                                 fInputObject.AddTrack(*mPart);
 417                         }
 418                 }
 419            break;
 420            default:
 421            break;
 422            //      delete mPart;
 423         }       //end of switch
 424 //      Info("FillTracks","After Particle Storage");
 425         //if (part%10) gObjectTable->Print();
 426
 427    }    //end of particle loop
 428    //gObjectTable->Print();
 429 }
 430
 431 void AliEMCALJetFinderInputSimPrep::FillPartons()
 432 {
 433         // Fill partons to
 434         // input object from simulation
 435 if (fDebug > 1) Info("FillPartons","Beginning FillPartons");
 436
 437   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 438   AliGenEventHeader* evHeader = rl->GetHeader()->GenEventHeader();
 439   Float_t triggerJetValues[4];
 440   AliEMCALParton tempParton;
 441
 442   if (fFileType == kPythia)
 443   {
 444         Int_t ntriggerjets = ((AliGenPythiaEventHeader*)evHeader)->NTriggerJets();
 445         if (fDebug > 1) Info("FillPartons","Number of trigger jets --> %i",ntriggerjets);
 446         for (Int_t jetcounter = 0 ; jetcounter < ntriggerjets; jetcounter++){
 447         ((AliGenPythiaEventHeader*)evHeader)->TriggerJet(jetcounter,triggerJetValues);
 448         TLorentzVector tempLParton;
 449         tempLParton.SetPxPyPzE(triggerJetValues[0],triggerJetValues[1],triggerJetValues[2],triggerJetValues[3]);
 450         tempParton.SetEnergy(tempLParton.Energy());
 451         tempParton.SetEta(tempLParton.Eta());
 452         tempParton.SetPhi(tempLParton.Phi());
 453         FillPartonTracks(&tempParton);
 454         fInputObject.AddParton(&tempParton);
 455         }
 456
 457   }
 458 }
 459
 460 void AliEMCALJetFinderInputSimPrep::FillParticles()
 461 {
 462         // Fill particles to input object from simulation
 463 if (fDebug > 1) Info("FillParticles","Beginning FillParticles");
 464
 465     Int_t npart = (gAlice->GetHeader())->GetNprimary();
 466     TParticlePDG* pdgP = 0;
 467
 468     for (Int_t part = 0; part < npart; part++) {
 469         TParticle *mPart = gAlice->GetMCApp()->Particle(part);
 470         pdgP = mPart->GetPDG();
 471
 472         if (fDebug > 10) Info("FillParticles","Checking if particle is a primary");
 473
 474         if (fFileType == kPythia) {
 475             if (mPart->GetStatusCode() != 1) continue;
 476         } else if (fFileType == kHijing) {
 477             if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
 478         }
 479
 480         if (fDebug > 10) Info("FillParticles","Checking if particle is in acceptance");
 481
 482         if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin))   continue;
 483         if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin))   continue;
 484
 485
 486 /*
 487         {kProton, kProtonBar, kElectron, kPositron,
 488          kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
 489          kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
 490          kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
 491          kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
 492          0,kNuMu,kNuMuBar
 493 */
 494         switch(fTrackType)
 495         {
 496
 497            case kAllP:  // All Stable particles to be included
 498                 if (fDebug > 5) Info("FillParticles","Storing particle");
 499                 fInputObject.AddParticle(mPart);
 500            break;
 501            case kEM:   // All Electromagnetic particles
 502                 if (mPart->GetPdgCode() == kElectron ||
 503                     mPart->GetPdgCode() == kPositron ||
 504                     mPart->GetPdgCode() == kGamma){
 505                         if (fDebug > 5) Info("FillParticles","Storing electromagnetic particle");
 506                         fInputObject.AddParticle(mPart);
 507                 }
 508            break;
 509            case kCharged: // All Particles with non-zero charge
 510                 if (pdgP->Charge() != 0) {
 511                         if (fDebug > 5) Info("FillParticles","Storing charged particle");
 512                         fInputObject.AddParticle(mPart);
 513                 }
 514            break;
 515            case kNeutral: // All particles with no charge
 516                 if (pdgP->Charge() == 0){
 517                         if (fDebug > 5) Info("FillParticles","Storing neutral particle");
 518                         fInputObject.AddParticle(mPart);
 519                 }
 520            break;
 521            case kHadron: //All hadrons
 522                 if (mPart->GetPdgCode() != kElectron &&
 523                     mPart->GetPdgCode() != kPositron &&
 524                     mPart->GetPdgCode() != kGamma )
 525                 {
 526
 527                 if (fDebug > 5) Info("FillParticles","Storing hadron");
 528                     fInputObject.AddParticle(mPart);
 529                 }
 530            break;
 531            case kChargedHadron:  // only charged hadrons
 532                 if (mPart->GetPdgCode() != kElectron &&
 533                     mPart->GetPdgCode() != kPositron &&
 534                     mPart->GetPdgCode() != kGamma    &&
 535                     pdgP->Charge()      != 0       ){
 536                 if (fDebug > 5) Info("FillParticles","Storing charged hadron");
 537                        fInputObject.AddParticle(mPart);
 538                 }
 539            break;
 540            case kNoTracks:
 541            break;
 542            case kEMChargedPi0:
 543                 if (pdgP->Charge() != 0 || mPart->GetPdgCode() == kPi0  ||
 544                         mPart->GetPdgCode() == kGamma     )
 545                 {
 546                         if (fDebug > 5) Info("FillTracks","Storing kEMChargedPi0 track");
 547                         if (fSmearType == kSmear ||
 548                             fSmearType == kSmearEffic ){
 549                                 Smear(mPart);/*
 550                                 TParticle *tmp = Smear(MPart);
 551                                 fInputObject.AddTrack(tmp);
 552                                 delete tmp;*/
 553                         }else{
 554                          fInputObject.AddTrack(*mPart);
 555                         }
 556                 }
 557            break;
 558            case kNoNeutronNeutrinoKlong:
 559                 if ( mPart->GetPdgCode() != kNeutron    &&
 560                      mPart->GetPdgCode() != kNeutronBar &&
 561                      mPart->GetPdgCode() != kK0Long     &&
 562                      mPart->GetPdgCode() != kNuE        &&
 563                      mPart->GetPdgCode() != kNuEBar     &&
 564                      mPart->GetPdgCode() != kNuMu       &&
 565                      mPart->GetPdgCode() != kNuMuBar    &&
 566                      mPart->GetPdgCode() != kNuTau      &&
 567                      mPart->GetPdgCode() != kNuTauBar   )
 568                 {
 569                         if (fDebug > 5) Info("FillTracks","Storing kNoNeutronNeutrinoKlong track");
 570                         if (fSmearType == kSmear ||
 571                                         fSmearType == kSmearEffic ){
 572                                 Smear(mPart);/*
 573                                                 TParticle *tmp = Smear(MPart);
 574                                                 fInputObject.AddTrack(tmp);
 575                                                 delete tmp;*/
 576                         }else{
 577                                 fInputObject.AddTrack(*mPart);
 578                         }
 579                 }
 580            break;
 581            default:
 582            break;
 583         }       //end of switch
 584     }// end of loop over particles
 585
 586 }
 587 void AliEMCALJetFinderInputSimPrep::FillDigits()
 588 {
 589         // Fill digits to input object
 590
 591 }
 592 void AliEMCALJetFinderInputSimPrep::FillSDigits()
 593 {
 594 Info("FillSDigits","Beginning FillSDigits");
 595         AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 596
 597         // Fill digits to input object
 598         for (Int_t towerid=0; towerid < emcalLoader->SDigits()->GetEntries(); towerid++)
 599         {
 600                 fInputObject.AddEnergyToDigit(emcalLoader->SDigit(towerid)->GetId(), emcalLoader->SDigit(towerid)->GetAmp());
 601         }
 602
 603 }
 604
 605 void AliEMCALJetFinderInputSimPrep::Smear(TParticle *particle)
 606 {
 607         // Smear particle momentum
 608 if (fDebug > 5) Info("Smear","Beginning Smear");
 609
 610         Float_t tmp = AliEMCALFast::SmearMomentum(1,particle->P());
 611         Float_t tmpt = (tmp/particle->P()) * particle->Pt();
 612         if (fDebug > 10) Info("Smear","Smearing changed P from %f to %f",particle->Pt(),tmpt);
 613         TLorentzVector tmpvector;
 614         tmpvector.SetPtEtaPhiM(tmpt,particle->Eta(), particle->Phi(), particle->GetMass());
 615 // create a new particle with smearing momentum - and no daughter or parent information and the same
 616 // vertex
 617
 618         TParticle tmparticle(particle->GetPdgCode(), 1, 0, 0,0,0,
 619                         tmpvector.Px()  , tmpvector.Py(), tmpvector.Pz(), tmpvector.Energy(),
 620                         particle->Vx(), particle->Vy(), particle->Vz(), particle->T());
 621         fInputObject.AddTrack(tmparticle);
 622         return;
 623 }
 624
 625 void AliEMCALJetFinderInputSimPrep::FillPartonTracks(AliEMCALParton *parton)
 626 {
 627         // Populate parton tracks for input distributions
 628   if (fDebug>1) Info("FillPartonTracks","Beginning FillPartonTracks");
 629   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 630
 631         Int_t npart = rl->Stack()->GetNprimary();
 632         Int_t ntracks =0;
 633         TParticlePDG *getpdg;
 634         TParticle *tempPart;
 635         for (Int_t part = 0; part < npart; part++)
 636         {
 637                 tempPart = rl->Stack()->Particle(part);
 638                 if (tempPart->GetStatusCode() != 1) continue;
 639                 if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
 640                     tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
 641                         if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
 642                         continue;
 643                 }
 644                 getpdg = tempPart->GetPDG();
 645                 if (getpdg->Charge() == 0.0  ) {
 646                         if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
 647                         continue;
 648                 }
 649                 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
 650                 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 ) ntracks++;
 651         }
 652         Float_t *energy = new Float_t[ntracks];
 653         Float_t *eta    = new Float_t[ntracks];
 654         Float_t *phi    = new Float_t[ntracks];
 655         Int_t   *pdg    = new Int_t[ntracks];
 656         ntracks=0;
 657         for (Int_t part = 0; part < npart; part++)
 658         {
 659                 tempPart = rl->Stack()->Particle(part);
 660                 if (tempPart->GetStatusCode() != 1) continue;
 661                 if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
 662                     tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
 663                         if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
 664                         continue;
 665                 }
 666                 if (tempPart->GetStatusCode() != 1) continue;
 667                 getpdg = tempPart->GetPDG();
 668                 if (getpdg->Charge() == 0.0  ) {
 669                         if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
 670                         continue;
 671                 }
 672                 if (   (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
 673                 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 )
 674                 {
 675                         energy[ntracks] = tempPart->Pt();
 676                         eta[ntracks] = tempPart->Eta();
 677                         phi[ntracks] = tempPart->Phi();
 678                         pdg[ntracks] = tempPart->GetPdgCode();
 679                         ntracks++;
 680                 }
 681         }
 682         parton->SetTrackList(ntracks,energy,eta,phi,pdg);
 683         delete[] energy;
 684         delete[] eta;
 685         delete[] phi;
 686         delete[] pdg;
 687 }
 688