EMCAL/jetfinder/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             /*  have to be tune for TRD1; May 31,06
 192             if (geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi) == -1 )
 193             {
 194                     if (fDebug >5)
 195                     {
 196                             Error("FillHits","Hit not inside EMCAL!");
 197                             mHit->Dump();
 198                     }
 199                     continue;
 200             }
 201             fInputObject.AddEnergyToDigit(geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi)-1,Int_t(1e7*etH));
 202             */
 203         } // Hit Loop
 204 //    } // Track Loop
 205
 206
 207 }
 208 void AliEMCALJetFinderInputSimPrep::FillTracks()
 209 {
 210         // Fill from particles simulating a TPC to input object from simulation
 211     if (fDebug > 1) Info("FillTracks","Beginning FillTracks");
 212
 213     AliRunLoader *rl = AliRunLoader::GetRunLoader();
 214     TParticlePDG* pdgP = 0;
 215     TParticle *mPart;
 216     Int_t npart = rl->Stack()->GetNprimary();
 217     if (fDebug > 1) Info("FillTracks","Header says there are %i primaries",npart);
 218     Float_t bfield,rEMCAL;
 219
 220     if (fDebug > 1) Info("FillTracks","Defining the geometry");
 221
 222     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 223     AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
 224     if (geom == 0 && pEMCAL)
 225       geom = AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
 226     if (geom == 0)
 227       geom = AliEMCALGeometry::GetInstance("EMCAL_55_25","");//","SHISH_77_TRD1_2X2_FINAL_110DEG");
 228     fEtaMax = geom->GetArm1EtaMax();
 229     fEtaMin = geom->GetArm1EtaMin();
 230     fPhiMax = TMath::Pi()/180.0*geom->GetArm1PhiMax();
 231     fPhiMin = TMath::Pi()/180.0*geom->GetArm1PhiMin();
 232
 233     if (fDebug > 1) Info("FillTracks","Starting particle loop");
 234
 235     if (fDebug > 1) Info("FillTracks","Particle loop of %i",npart);
 236     for (Int_t part = 0; part < npart; part++) {
 237         mPart = rl->Stack()->Particle(part);
 238         pdgP = mPart->GetPDG();
 239
 240         if (fDebug > 5) Info("FillTracks","Checking if track is a primary");
 241
 242         if (fFileType == kPythia) {
 243             if (mPart->GetStatusCode() != 1) continue;
 244         } else if (fFileType == kHijing) {
 245             if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
 246         }
 247
 248         if (fDebug > 15) Info("FillTracks","Checking if track (eta - %f, phi - %f) is in acceptance",mPart->Eta(),mPart->Phi());
 249         if (fDebug > 10) Info("FillTracks","Checking if EMCAL acceptance  ( %f < eta < %f, %f < phi < %f) is in acceptance",fEtaMin,fEtaMax,fPhiMin,fPhiMax);
 250
 251         if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin))   continue;
 252         if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin))   continue;
 253
 254 /*
 255         {kProton, kProtonBar, kElectron, kPositron,
 256          kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
 257          kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
 258          kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
 259          kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
 260          0,kNuMu,kNuMuBar
 261 */
 262
 263         if (fDebug > 15) Info("FillTracks","Checking if track is rejected");
 264
 265         if ((fSmearType == kEfficiency  ||
 266              fSmearType == kSmearEffic) &&
 267              pdgP->Charge()!=0) {
 268             if (AliEMCALFast::RandomReject(fEfficiency)) {
 269                 continue;
 270             }
 271         }
 272
 273         if (gAlice && gAlice->Field())
 274           bfield = gAlice->Field()->SolenoidField();
 275         else
 276           bfield = 0.4;
 277         rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
 278         Float_t rB = 3335.6 * mPart->Pt() / bfield;  // [cm]  (case of |charge|=1)
 279         if (2.*rB < rEMCAL) continue;  // track curls away
 280
 281         //if (part%10) gObjectTable->Print();
 282         switch(fTrackType)
 283         {
 284
 285            case kAllP:  // All Stable particles to be included
 286                 if (fDebug > 5) Info("FillTracks","Storing track");
 287                 if (fSmearType == kSmear ||
 288                     fSmearType == kSmearEffic ){
 289                         Smear(mPart);/*
 290                         TParticle *tmp = Smear(MPart);
 291                         fInputObject.AddTrack(Smear(MPart));
 292                         delete tmp;*/
 293                 }else{
 294                         fInputObject.AddTrack(*mPart);
 295                 }
 296            break;
 297            case kEM:   // All Electromagnetic particles
 298                 if (mPart->GetPdgCode() == kElectron  ||
 299                     mPart->GetPdgCode() == kMuonPlus  ||
 300                     mPart->GetPdgCode() == kMuonMinus ||
 301                     mPart->GetPdgCode() == kPositron ){
 302                       if (fDebug > 5) Info("FillTracks","Storing electron or positron");
 303                       if (fSmearType == kSmear ||
 304                             fSmearType == kSmearEffic ){
 305                               Smear(mPart);/*
 306                               TParticle *tmp = Smear(MPart);
 307                               fInputObject.AddTrack(tmp);
 308                               delete tmp;*/
 309                       }else{
 310                           fInputObject.AddTrack(*mPart);
 311                       }
 312                 }
 313                 if ( mPart->GetPdgCode() == kGamma ){
 314                         fInputObject.AddTrack(*mPart);
 315                         if (fDebug > 5) Info("FillTracks","Storing gamma");
 316                 }
 317
 318            break;
 319            case kCharged: // All Particles with non-zero charge
 320                 if (pdgP->Charge() != 0) {
 321                         if (fDebug > 5) Info("FillTracks","Storing charged track");
 322                         if (fSmearType == kSmear ||
 323                             fSmearType == kSmearEffic ){
 324                                 Smear(mPart);/*
 325                                 TParticle *tmp = Smear(MPart);
 326                                 fInputObject.AddTrack(tmp);
 327                                 delete tmp;*/
 328                         }else{
 329                          fInputObject.AddTrack(*mPart);
 330                         }
 331                 }
 332            break;
 333            case kNeutral: // All particles with no charge
 334                 if (pdgP->Charge() == 0){
 335                         fInputObject.AddTrack(*mPart);
 336                         if (fDebug > 5) Info("FillTracks","Storing neutral");
 337                 }
 338            break;
 339            case kHadron: //All hadrons
 340                 if (mPart->GetPdgCode() != kElectron  &&
 341                     mPart->GetPdgCode() != kPositron  &&
 342                     mPart->GetPdgCode() != kMuonPlus  &&
 343                     mPart->GetPdgCode() != kMuonMinus &&
 344                     mPart->GetPdgCode() != kGamma )
 345                 {
 346                         if (fDebug > 5) Info("FillTracks","Storing hadron");
 347                         if (pdgP->Charge() == 0){
 348                             fInputObject.AddTrack(*mPart);
 349                         }else{
 350                                 if (fSmearType == kSmear ||
 351                                     fSmearType == kSmearEffic ){
 352                                         Smear(mPart);/*
 353                                         TParticle *tmp = Smear(MPart);
 354                                         fInputObject.AddTrack(tmp);
 355                                         delete tmp;*/
 356                                 }else{
 357                                     fInputObject.AddTrack(*mPart);
 358                                 }
 359                         }
 360                 }
 361            break;
 362            case kChargedHadron:  // only charged hadrons
 363                 if (mPart->GetPdgCode() != kElectron  &&
 364                     mPart->GetPdgCode() != kPositron  &&
 365                     mPart->GetPdgCode() != kGamma     &&
 366                     mPart->GetPdgCode() != kMuonPlus  &&
 367                     mPart->GetPdgCode() != kMuonMinus &&
 368                     pdgP->Charge()      != 0       ){
 369                         if (fDebug > 5) Info("FillTracks","Storing charged hadron");
 370                         if (fSmearType == kSmear ||
 371                             fSmearType == kSmearEffic ){
 372                                 Smear(mPart);/*
 373                                 TParticle *tmp = Smear(MPart);
 374                                 fInputObject.AddTrack(tmp);
 375                                 delete tmp;*/
 376                         }else{
 377                                fInputObject.AddTrack(*mPart);
 378                         }
 379                 }
 380            break;
 381            case kNoTracks:
 382            break;
 383            case kEMChargedPi0:
 384                 if (pdgP->Charge() != 0 || mPart->GetPdgCode() == kPi0  ||
 385                         mPart->GetPdgCode() == kGamma     )
 386                 {
 387                         if (fDebug > 5) Info("FillTracks","Storing charged track");
 388                         if (fSmearType == kSmear ||
 389                             fSmearType == kSmearEffic ){
 390                                 Smear(mPart);/*
 391                                 TParticle *tmp = Smear(MPart);
 392                                 fInputObject.AddTrack(tmp);
 393                                 delete tmp;*/
 394                         }else{
 395                          fInputObject.AddTrack(*mPart);
 396                         }
 397                 }
 398            break;
 399            case kNoNeutronNeutrinoKlong:
 400                 if ( mPart->GetPdgCode() != kNeutron    &&
 401                      mPart->GetPdgCode() != kNeutronBar &&
 402                      mPart->GetPdgCode() != kK0Long     &&
 403                      mPart->GetPdgCode() != kNuE        &&
 404                      mPart->GetPdgCode() != kNuEBar     &&
 405                      mPart->GetPdgCode() != kNuMu       &&
 406                      mPart->GetPdgCode() != kNuMuBar    &&
 407                      mPart->GetPdgCode() != kNuTau      &&
 408                      mPart->GetPdgCode() != kNuTauBar   )
 409                 {
 410                         if (fDebug > 5) Info("FillTracks","Storing charged track");
 411                         if (fSmearType == kSmear ||
 412                                         fSmearType == kSmearEffic ){
 413                                 Smear(mPart);/*
 414                                                 TParticle *tmp = Smear(MPart);
 415                                                 fInputObject.AddTrack(tmp);
 416                                                 delete tmp;*/
 417                         }else{
 418                                 fInputObject.AddTrack(*mPart);
 419                         }
 420                 }
 421            break;
 422            default:
 423            break;
 424            //      delete mPart;
 425         }       //end of switch
 426 //      Info("FillTracks","After Particle Storage");
 427         //if (part%10) gObjectTable->Print();
 428
 429    }    //end of particle loop
 430    //gObjectTable->Print();
 431 }
 432
 433 void AliEMCALJetFinderInputSimPrep::FillPartons()
 434 {
 435         // Fill partons to
 436         // input object from simulation
 437 if (fDebug > 1) Info("FillPartons","Beginning FillPartons");
 438
 439   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 440   AliGenEventHeader* evHeader = rl->GetHeader()->GenEventHeader();
 441   Float_t triggerJetValues[4];
 442   AliEMCALParton tempParton;
 443
 444   if (fFileType == kPythia)
 445   {
 446         Int_t ntriggerjets = ((AliGenPythiaEventHeader*)evHeader)->NTriggerJets();
 447         if (fDebug > 1) Info("FillPartons","Number of trigger jets --> %i",ntriggerjets);
 448         for (Int_t jetcounter = 0 ; jetcounter < ntriggerjets; jetcounter++){
 449         ((AliGenPythiaEventHeader*)evHeader)->TriggerJet(jetcounter,triggerJetValues);
 450         TLorentzVector tempLParton;
 451         tempLParton.SetPxPyPzE(triggerJetValues[0],triggerJetValues[1],triggerJetValues[2],triggerJetValues[3]);
 452         tempParton.SetEnergy(tempLParton.Energy());
 453         tempParton.SetEta(tempLParton.Eta());
 454         tempParton.SetPhi(tempLParton.Phi());
 455         FillPartonTracks(&tempParton);
 456         fInputObject.AddParton(&tempParton);
 457         }
 458
 459   }
 460 }
 461
 462 void AliEMCALJetFinderInputSimPrep::FillParticles()
 463 {
 464         // Fill particles to input object from simulation
 465 if (fDebug > 1) Info("FillParticles","Beginning FillParticles");
 466
 467     Int_t npart = (gAlice->GetHeader())->GetNprimary();
 468     TParticlePDG* pdgP = 0;
 469
 470     for (Int_t part = 0; part < npart; part++) {
 471         TParticle *mPart = gAlice->GetMCApp()->Particle(part);
 472         pdgP = mPart->GetPDG();
 473
 474         if (fDebug > 10) Info("FillParticles","Checking if particle is a primary");
 475
 476         if (fFileType == kPythia) {
 477             if (mPart->GetStatusCode() != 1) continue;
 478         } else if (fFileType == kHijing) {
 479             if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
 480         }
 481
 482         if (fDebug > 10) Info("FillParticles","Checking if particle is in acceptance");
 483
 484         if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin))   continue;
 485         if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin))   continue;
 486
 487
 488 /*
 489         {kProton, kProtonBar, kElectron, kPositron,
 490          kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
 491          kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
 492          kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
 493          kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
 494          0,kNuMu,kNuMuBar
 495 */
 496         switch(fTrackType)
 497         {
 498
 499            case kAllP:  // All Stable particles to be included
 500                 if (fDebug > 5) Info("FillParticles","Storing particle");
 501                 fInputObject.AddParticle(mPart);
 502            break;
 503            case kEM:   // All Electromagnetic particles
 504                 if (mPart->GetPdgCode() == kElectron ||
 505                     mPart->GetPdgCode() == kPositron ||
 506                     mPart->GetPdgCode() == kGamma){
 507                         if (fDebug > 5) Info("FillParticles","Storing electromagnetic particle");
 508                         fInputObject.AddParticle(mPart);
 509                 }
 510            break;
 511            case kCharged: // All Particles with non-zero charge
 512                 if (pdgP->Charge() != 0) {
 513                         if (fDebug > 5) Info("FillParticles","Storing charged particle");
 514                         fInputObject.AddParticle(mPart);
 515                 }
 516            break;
 517            case kNeutral: // All particles with no charge
 518                 if (pdgP->Charge() == 0){
 519                         if (fDebug > 5) Info("FillParticles","Storing neutral particle");
 520                         fInputObject.AddParticle(mPart);
 521                 }
 522            break;
 523            case kHadron: //All hadrons
 524                 if (mPart->GetPdgCode() != kElectron &&
 525                     mPart->GetPdgCode() != kPositron &&
 526                     mPart->GetPdgCode() != kGamma )
 527                 {
 528
 529                 if (fDebug > 5) Info("FillParticles","Storing hadron");
 530                     fInputObject.AddParticle(mPart);
 531                 }
 532            break;
 533            case kChargedHadron:  // only charged hadrons
 534                 if (mPart->GetPdgCode() != kElectron &&
 535                     mPart->GetPdgCode() != kPositron &&
 536                     mPart->GetPdgCode() != kGamma    &&
 537                     pdgP->Charge()      != 0       ){
 538                 if (fDebug > 5) Info("FillParticles","Storing charged hadron");
 539                        fInputObject.AddParticle(mPart);
 540                 }
 541            break;
 542            case kNoTracks:
 543            break;
 544            case kEMChargedPi0:
 545                 if (pdgP->Charge() != 0 || mPart->GetPdgCode() == kPi0  ||
 546                         mPart->GetPdgCode() == kGamma     )
 547                 {
 548                         if (fDebug > 5) Info("FillTracks","Storing kEMChargedPi0 track");
 549                         if (fSmearType == kSmear ||
 550                             fSmearType == kSmearEffic ){
 551                                 Smear(mPart);/*
 552                                 TParticle *tmp = Smear(MPart);
 553                                 fInputObject.AddTrack(tmp);
 554                                 delete tmp;*/
 555                         }else{
 556                          fInputObject.AddTrack(*mPart);
 557                         }
 558                 }
 559            break;
 560            case kNoNeutronNeutrinoKlong:
 561                 if ( mPart->GetPdgCode() != kNeutron    &&
 562                      mPart->GetPdgCode() != kNeutronBar &&
 563                      mPart->GetPdgCode() != kK0Long     &&
 564                      mPart->GetPdgCode() != kNuE        &&
 565                      mPart->GetPdgCode() != kNuEBar     &&
 566                      mPart->GetPdgCode() != kNuMu       &&
 567                      mPart->GetPdgCode() != kNuMuBar    &&
 568                      mPart->GetPdgCode() != kNuTau      &&
 569                      mPart->GetPdgCode() != kNuTauBar   )
 570                 {
 571                         if (fDebug > 5) Info("FillTracks","Storing kNoNeutronNeutrinoKlong track");
 572                         if (fSmearType == kSmear ||
 573                                         fSmearType == kSmearEffic ){
 574                                 Smear(mPart);/*
 575                                                 TParticle *tmp = Smear(MPart);
 576                                                 fInputObject.AddTrack(tmp);
 577                                                 delete tmp;*/
 578                         }else{
 579                                 fInputObject.AddTrack(*mPart);
 580                         }
 581                 }
 582            break;
 583            default:
 584            break;
 585         }       //end of switch
 586     }// end of loop over particles
 587
 588 }
 589 void AliEMCALJetFinderInputSimPrep::FillDigits()
 590 {
 591         // Fill digits to input object
 592
 593 }
 594 void AliEMCALJetFinderInputSimPrep::FillSDigits()
 595 {
 596 Info("FillSDigits","Beginning FillSDigits");
 597         AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 598
 599         // Fill digits to input object
 600         for (Int_t towerid=0; towerid < emcalLoader->SDigits()->GetEntries(); towerid++)
 601         {
 602                 fInputObject.AddEnergyToDigit(emcalLoader->SDigit(towerid)->GetId(), emcalLoader->SDigit(towerid)->GetAmp());
 603         }
 604
 605 }
 606
 607 void AliEMCALJetFinderInputSimPrep::Smear(TParticle *particle)
 608 {
 609         // Smear particle momentum
 610 if (fDebug > 5) Info("Smear","Beginning Smear");
 611
 612         Float_t tmp = AliEMCALFast::SmearMomentum(1,particle->P());
 613         Float_t tmpt = (tmp/particle->P()) * particle->Pt();
 614         if (fDebug > 10) Info("Smear","Smearing changed P from %f to %f",particle->Pt(),tmpt);
 615         TLorentzVector tmpvector;
 616         tmpvector.SetPtEtaPhiM(tmpt,particle->Eta(), particle->Phi(), particle->GetMass());
 617 // create a new particle with smearing momentum - and no daughter or parent information and the same
 618 // vertex
 619
 620         TParticle tmparticle(particle->GetPdgCode(), 1, 0, 0,0,0,
 621                         tmpvector.Px()  , tmpvector.Py(), tmpvector.Pz(), tmpvector.Energy(),
 622                         particle->Vx(), particle->Vy(), particle->Vz(), particle->T());
 623         fInputObject.AddTrack(tmparticle);
 624         return;
 625 }
 626
 627 void AliEMCALJetFinderInputSimPrep::FillPartonTracks(AliEMCALParton *parton)
 628 {
 629         // Populate parton tracks for input distributions
 630   if (fDebug>1) Info("FillPartonTracks","Beginning FillPartonTracks");
 631   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 632
 633         Int_t npart = rl->Stack()->GetNprimary();
 634         Int_t ntracks =0;
 635         TParticlePDG *getpdg;
 636         TParticle *tempPart;
 637         for (Int_t part = 0; part < npart; part++)
 638         {
 639                 tempPart = rl->Stack()->Particle(part);
 640                 if (tempPart->GetStatusCode() != 1) continue;
 641                 if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
 642                     tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
 643                         if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
 644                         continue;
 645                 }
 646                 getpdg = tempPart->GetPDG();
 647                 if (getpdg->Charge() == 0.0  ) {
 648                         if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
 649                         continue;
 650                 }
 651                 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
 652                 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 ) ntracks++;
 653         }
 654         Float_t *energy = new Float_t[ntracks];
 655         Float_t *eta    = new Float_t[ntracks];
 656         Float_t *phi    = new Float_t[ntracks];
 657         Int_t   *pdg    = new Int_t[ntracks];
 658         ntracks=0;
 659         for (Int_t part = 0; part < npart; part++)
 660         {
 661                 tempPart = rl->Stack()->Particle(part);
 662                 if (tempPart->GetStatusCode() != 1) continue;
 663                 if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
 664                     tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
 665                         if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
 666                         continue;
 667                 }
 668                 if (tempPart->GetStatusCode() != 1) continue;
 669                 getpdg = tempPart->GetPDG();
 670                 if (getpdg->Charge() == 0.0  ) {
 671                         if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
 672                         continue;
 673                 }
 674                 if (   (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
 675                 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 )
 676                 {
 677                         energy[ntracks] = tempPart->Pt();
 678                         eta[ntracks] = tempPart->Eta();
 679                         phi[ntracks] = tempPart->Phi();
 680                         pdg[ntracks] = tempPart->GetPdgCode();
 681                         ntracks++;
 682                 }
 683         }
 684         parton->SetTrackList(ntracks,energy,eta,phi,pdg);
 685         delete[] energy;
 686         delete[] eta;
 687         delete[] phi;
 688         delete[] pdg;
 689 }
 690