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