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,TString data="H")
 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    if (data.Contains("S"))
 120    {
 121            gime->Event(EventNumber,"XS") ;
 122    }else
 123    {
 124            gime->Event(EventNumber,"XH") ;
 125    }
 126    if (fDebug > 1) Info("FillFromFile","Got event %i with option \"XH\"",EventNumber);
 127
 128    if ( fEMCALType == kHits ||
 129         fEMCALType == kTimeCut )
 130    {
 131            if (data.Contains("S"))
 132            {
 133                    FillSDigits();
 134            }else
 135            {
 136                    FillHits();
 137            }
 138    }
 139    if ( fTrackType != kNoTracks  )  FillTracks();
 140    if ( fFileType  != kData){
 141            FillPartons();
 142    }
 143 /*   gime->EmcalLoader()->UnloadRecParticles();
 144    gime->EmcalLoader()->UnloadTracks();
 145    gime->Reset();*/
 146    return 0;
 147 }
 148
 149 void AliEMCALJetFinderInputSimPrep::FillHits()
 150 {
 151         // Fill from the hits to input object from simulation
 152 if (fDebug > 1) Info("FillHits","Beginning FillHits");
 153
 154 // Access hit information
 155     AliEMCALHit *mHit;
 156 //  AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 157 //  Info("AliEMCALJetFinderInputSimPrep","Title of the geometry is %s",pEMCAL->GetTitle());
 158     AliEMCALGeometry* geom =  AliEMCALGetter::Instance()->EMCALGeometry();
 159     AliEMCALGetter *gime = AliEMCALGetter::Instance();
 160
 161     //    TTree *treeH = AliEMCALGetter::Instance()->TreeH();
 162 //    Int_t ntracks = (Int_t) treeH->GetEntries();
 163 //
 164 //   Loop over tracks
 165 //
 166     Double_t etH = 0.0;
 167
 168 /*    for (Int_t track=0; track<ntracks;track++) {
 169         gAlice->ResetHits();
 170         nbytes += treeH->GetEvent(track);*/
 171 //
 172 //  Loop over hits
 173 //
 174         for(Int_t i =0; i < gime->Hits()->GetEntries() ;i++)
 175         {
 176             mHit = gime->Hit(i);
 177                 if (fEMCALType == kTimeCut &&
 178                     (mHit->GetTime()>fTimeCut) ) continue;
 179             Float_t x      =    mHit->X();         // x-pos of hit
 180             Float_t y      =    mHit->Y();         // y-pos
 181             Float_t z      =    mHit->Z();         // z-pos
 182             Float_t eloss  =    mHit->GetEnergy(); // deposited energy
 183
 184             Float_t r      =    TMath::Sqrt(x*x+y*y);
 185             Float_t theta  =    TMath::ATan2(r,z);
 186             Float_t eta    =   -TMath::Log(TMath::Tan(theta/2.));
 187             Float_t phi    =    TMath::ATan2(y,x);
 188             etH = eloss*TMath::Sin(theta);
 189             if (fDebug > 10) Info("FillHits","Values of hit energy %i",Int_t(1e7*etH));
 190             if (geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi) == -1)
 191             {
 192                     if (fDebug >5)
 193                     {
 194                             Error("FillHits","Hit not inside EMCAL!");
 195                             mHit->Dump();
 196                     }
 197                     continue;
 198             }
 199             fInputObject.AddEnergyToDigit(geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi)-1,Int_t(1e7*etH));
 200         } // Hit Loop
 201 //    } // Track Loop
 202
 203
 204 }
 205 void AliEMCALJetFinderInputSimPrep::FillTracks()
 206 {
 207         // Fill from particles simulating a TPC to input object from simulation
 208
 209     if (fDebug > 1) Info("FillTracks","Beginning FillTracks");
 210
 211     TParticlePDG* pdgP = 0;
 212     TParticle *mPart;
 213     Int_t npart = (gAlice->GetHeader())->GetNprimary();
 214     if (fDebug > 1) Info("FillTracks","Header says there are %i primaries",npart);
 215     //AliEMCALGetter *gime = AliEMCALGetter::Instance();
 216     Float_t bfield,rEMCAL;
 217 /*      for (Int_t i =0; i<gime->NPrimaries(); i++)
 218         {
 219                 cout<<gime->Primary(i)->GetFirstMother()<<endl;
 220                 if (gime->Primary(i)->GetFirstMother() == -1) {
 221                         counterforus++;
 222                         gime->Primary(i)->Dump();
 223                 }
 224         }*/
 225
 226     if (fDebug > 1) Info("FillTracks","Defining the geometry");
 227
 228     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 229     AliEMCALGeometry* geom =  AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
 230     fEtaMax = geom->GetArm1EtaMax();
 231     fEtaMin = geom->GetArm1EtaMin();
 232     fPhiMax = TMath::Pi()/180.0*geom->GetArm1PhiMax();
 233     fPhiMin = TMath::Pi()/180.0*geom->GetArm1PhiMin();
 234
 235     if (fDebug > 1) Info("FillTracks","Starting particle loop");
 236
 237     if (fDebug > 1) Info("FillTracks","Particle loop of %i",npart);
 238     for (Int_t part = 0; part < npart; part++) {
 239         mPart = gAlice->GetMCApp()->Particle(part);
 240         //if (part%10) gObjectTable->Print();
 241         pdgP = mPart->GetPDG();
 242
 243         if (fDebug > 5) Info("FillTracks","Checking if track is a primary");
 244
 245         if (fFileType == kPythia) {
 246             if (mPart->GetStatusCode() != 1) continue;
 247         } else if (fFileType == kHijing) {
 248             if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
 249         }
 250
 251         if (fDebug > 15) Info("FillTracks","Checking if track (eta - %f, phi - %f) is in acceptance",mPart->Eta(),mPart->Phi());
 252         if (fDebug > 10) Info("FillTracks","Checking if EMCAL acceptance  ( %f < eta < %f, %f < phi < %f) is in acceptance",fEtaMin,fEtaMax,fPhiMin,fPhiMax);
 253
 254         if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin))   continue;
 255         if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin))   continue;
 256
 257 /*
 258         {kProton, kProtonBar, kElectron, kPositron,
 259          kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
 260          kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
 261          kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
 262          kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
 263          0,kNuMu,kNuMuBar
 264 */
 265
 266         if (fDebug > 15) Info("FillTracks","Checking if track is rejected");
 267
 268         if ((fSmearType == kEfficiency  ||
 269              fSmearType == kSmearEffic) &&
 270              pdgP->Charge()!=0) {
 271             if (AliEMCALFast::RandomReject(fEfficiency)) {
 272                 continue;
 273             }
 274         }
 275
 276         bfield = gAlice->Field()->SolenoidField();
 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   AliGenEventHeader* evHeader = ((AliHeader*)(gAlice->GetHeader()))->GenEventHeader();
 440   Float_t triggerJetValues[4];
 441   AliEMCALParton tempParton;
 442
 443   if (fFileType == kPythia)
 444   {
 445         Int_t ntriggerjets = ((AliGenPythiaEventHeader*)evHeader)->NTriggerJets();
 446         if (fDebug > 1) Info("FillPartons","Number of trigger jets --> %i",ntriggerjets);
 447         for (Int_t jetcounter = 0 ; jetcounter < ntriggerjets; jetcounter++){
 448         ((AliGenPythiaEventHeader*)evHeader)->TriggerJet(jetcounter,triggerJetValues);
 449         TLorentzVector tempLParton;
 450         tempLParton.SetPxPyPzE(triggerJetValues[0],triggerJetValues[1],triggerJetValues[2],triggerJetValues[3]);
 451         tempParton.SetEnergy(tempLParton.Energy());
 452         tempParton.SetEta(tempLParton.Eta());
 453         tempParton.SetPhi(tempLParton.Phi());
 454         FillPartonTracks(&tempParton);
 455         fInputObject.AddParton(&tempParton);
 456         }
 457
 458   }
 459 }
 460
 461 void AliEMCALJetFinderInputSimPrep::FillParticles()
 462 {
 463         // Fill particles to input object from simulation
 464 if (fDebug > 1) Info("FillParticles","Beginning FillParticles");
 465
 466     Int_t npart = (gAlice->GetHeader())->GetNprimary();
 467     TParticlePDG* pdgP = 0;
 468
 469     for (Int_t part = 0; part < npart; part++) {
 470         TParticle *mPart = gAlice->GetMCApp()->Particle(part);
 471         pdgP = mPart->GetPDG();
 472
 473         if (fDebug > 10) Info("FillParticles","Checking if particle is a primary");
 474
 475         if (fFileType == kPythia) {
 476             if (mPart->GetStatusCode() != 1) continue;
 477         } else if (fFileType == kHijing) {
 478             if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
 479         }
 480
 481         if (fDebug > 10) Info("FillParticles","Checking if particle is in acceptance");
 482
 483         if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin))   continue;
 484         if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin))   continue;
 485
 486
 487 /*
 488         {kProton, kProtonBar, kElectron, kPositron,
 489          kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
 490          kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
 491          kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
 492          kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
 493          0,kNuMu,kNuMuBar
 494 */
 495         switch(fTrackType)
 496         {
 497
 498            case kAllP:  // All Stable particles to be included
 499                 if (fDebug > 5) Info("FillParticles","Storing particle");
 500                 fInputObject.AddParticle(mPart);
 501            break;
 502            case kEM:   // All Electromagnetic particles
 503                 if (mPart->GetPdgCode() == kElectron ||
 504                     mPart->GetPdgCode() == kPositron ||
 505                     mPart->GetPdgCode() == kGamma){
 506                         if (fDebug > 5) Info("FillParticles","Storing electromagnetic particle");
 507                         fInputObject.AddParticle(mPart);
 508                 }
 509            break;
 510            case kCharged: // All Particles with non-zero charge
 511                 if (pdgP->Charge() != 0) {
 512                         if (fDebug > 5) Info("FillParticles","Storing charged particle");
 513                         fInputObject.AddParticle(mPart);
 514                 }
 515            break;
 516            case kNeutral: // All particles with no charge
 517                 if (pdgP->Charge() == 0){
 518                         if (fDebug > 5) Info("FillParticles","Storing neutral particle");
 519                         fInputObject.AddParticle(mPart);
 520                 }
 521            break;
 522            case kHadron: //All hadrons
 523                 if (mPart->GetPdgCode() != kElectron &&
 524                     mPart->GetPdgCode() != kPositron &&
 525                     mPart->GetPdgCode() != kGamma )
 526                 {
 527
 528                 if (fDebug > 5) Info("FillParticles","Storing hadron");
 529                     fInputObject.AddParticle(mPart);
 530                 }
 531            break;
 532            case kChargedHadron:  // only charged hadrons
 533                 if (mPart->GetPdgCode() != kElectron &&
 534                     mPart->GetPdgCode() != kPositron &&
 535                     mPart->GetPdgCode() != kGamma    &&
 536                     pdgP->Charge()      != 0       ){
 537                 if (fDebug > 5) Info("FillParticles","Storing charged hadron");
 538                        fInputObject.AddParticle(mPart);
 539                 }
 540            break;
 541            case kNoTracks:
 542            break;
 543            case kEMChargedPi0:
 544                 if (pdgP->Charge() != 0 || mPart->GetPdgCode() == kPi0  ||
 545                         mPart->GetPdgCode() == kGamma     )
 546                 {
 547                         if (fDebug > 5) Info("FillTracks","Storing kEMChargedPi0 track");
 548                         if (fSmearType == kSmear ||
 549                             fSmearType == kSmearEffic ){
 550                                 Smear(mPart);/*
 551                                 TParticle *tmp = Smear(MPart);
 552                                 fInputObject.AddTrack(tmp);
 553                                 delete tmp;*/
 554                         }else{
 555                          fInputObject.AddTrack(*mPart);
 556                         }
 557                 }
 558            break;
 559            case kNoNeutronNeutrinoKlong:
 560                 if ( mPart->GetPdgCode() != kNeutron    &&
 561                      mPart->GetPdgCode() != kNeutronBar &&
 562                      mPart->GetPdgCode() != kK0Long     &&
 563                      mPart->GetPdgCode() != kNuE        &&
 564                      mPart->GetPdgCode() != kNuEBar     &&
 565                      mPart->GetPdgCode() != kNuMu       &&
 566                      mPart->GetPdgCode() != kNuMuBar    &&
 567                      mPart->GetPdgCode() != kNuTau      &&
 568                      mPart->GetPdgCode() != kNuTauBar   )
 569                 {
 570                         if (fDebug > 5) Info("FillTracks","Storing kNoNeutronNeutrinoKlong track");
 571                         if (fSmearType == kSmear ||
 572                                         fSmearType == kSmearEffic ){
 573                                 Smear(mPart);/*
 574                                                 TParticle *tmp = Smear(MPart);
 575                                                 fInputObject.AddTrack(tmp);
 576                                                 delete tmp;*/
 577                         }else{
 578                                 fInputObject.AddTrack(*mPart);
 579                         }
 580                 }
 581            break;
 582            default:
 583            break;
 584         }       //end of switch
 585     }// end of loop over particles
 586
 587 }
 588 void AliEMCALJetFinderInputSimPrep::FillDigits()
 589 {
 590         // Fill digits to input object
 591
 592 }
 593 void AliEMCALJetFinderInputSimPrep::FillSDigits()
 594 {
 595 Info("FillSDigits","Beginning FillSDigits");
 596         AliEMCALGetter *gime = AliEMCALGetter::Instance();
 597
 598         // Fill digits to input object
 599         for (Int_t towerid=0; towerid < gime->SDigits()->GetEntries(); towerid++)
 600         {
 601                 fInputObject.AddEnergyToDigit(gime->SDigit(towerid)->GetId(), gime->SDigit(towerid)->GetAmp());
 602         }
 603
 604 }
 605
 606 void AliEMCALJetFinderInputSimPrep::Smear(TParticle *particle)
 607 {
 608         // Smear particle momentum
 609 if (fDebug > 5) Info("Smear","Beginning Smear");
 610
 611         Float_t tmp = AliEMCALFast::SmearMomentum(1,particle->P());
 612         Float_t tmpt = (tmp/particle->P()) * particle->Pt();
 613         if (fDebug > 10) Info("Smear","Smearing changed P from %f to %f",particle->Pt(),tmpt);
 614         TLorentzVector tmpvector;
 615         tmpvector.SetPtEtaPhiM(tmpt,particle->Eta(), particle->Phi(), particle->GetMass());
 616 // create a new particle with smearing momentum - and no daughter or parent information and the same
 617 // vertex
 618
 619         TParticle tmparticle(particle->GetPdgCode(), 1, 0, 0,0,0,
 620                         tmpvector.Px()  , tmpvector.Py(), tmpvector.Pz(), tmpvector.Energy(),
 621                         particle->Vx(), particle->Vy(), particle->Vz(), particle->T());
 622         fInputObject.AddTrack(tmparticle);
 623         return;
 624 }
 625
 626 void AliEMCALJetFinderInputSimPrep::FillPartonTracks(AliEMCALParton *parton)
 627 {
 628         // Populate parton tracks for input distributions
 629         if (fDebug>1) Info("FillPartonTracks","Beginning FillPartonTracks");
 630         Int_t npart = (gAlice->GetHeader())->GetNprimary();
 631         Int_t ntracks =0;
 632         TParticlePDG *getpdg;
 633         TParticle *tempPart;
 634         for (Int_t part = 0; part < npart; part++)
 635         {
 636                 tempPart = gAlice->GetMCApp()->Particle(part);
 637                 if (tempPart->GetStatusCode() != 1) continue;
 638                 if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
 639                     tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
 640                         if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
 641                         continue;
 642                 }
 643                 getpdg = tempPart->GetPDG();
 644                 if (getpdg->Charge() == 0.0  ) {
 645                         if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
 646                         continue;
 647                 }
 648                 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
 649                 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 ) ntracks++;
 650         }
 651         Float_t *energy = new Float_t[ntracks];
 652         Float_t *eta    = new Float_t[ntracks];
 653         Float_t *phi    = new Float_t[ntracks];
 654         Int_t   *pdg    = new Int_t[ntracks];
 655         ntracks=0;
 656         for (Int_t part = 0; part < npart; part++)
 657         {
 658                 tempPart = gAlice->GetMCApp()->Particle(part);
 659                 if (tempPart->GetStatusCode() != 1) continue;
 660                 if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
 661                     tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
 662                         if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
 663                         continue;
 664                 }
 665                 if (tempPart->GetStatusCode() != 1) continue;
 666                 getpdg = tempPart->GetPDG();
 667                 if (getpdg->Charge() == 0.0  ) {
 668                         if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
 669                         continue;
 670                 }
 671                 if (   (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
 672                 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 )
 673                 {
 674                         energy[ntracks] = tempPart->Pt();
 675                         eta[ntracks] = tempPart->Eta();
 676                         phi[ntracks] = tempPart->Phi();
 677                         pdg[ntracks] = tempPart->GetPdgCode();
 678                         ntracks++;
 679                 }
 680         }
 681         parton->SetTrackList(ntracks,energy,eta,phi,pdg);
 682         delete[] energy;
 683         delete[] eta;
 684         delete[] phi;
 685         delete[] pdg;
 686 }
 687