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 "AliGenEventHeader.h"
  46 #include "AliGenPythiaEventHeader.h"
  47 #include "AliGenerator.h"
  48 #include "AliHeader.h"
  49
  50
  51 ClassImp(AliEMCALJetFinderInputSimPrep)
  52
  53 AliEMCALJetFinderInputSimPrep::AliEMCALJetFinderInputSimPrep()
  54 {
  55 if (fDebug > 0) Info("AliEMCALJetFinderInputSimPrep","Beginning Constructor");
  56
  57   fDebug = 0;
  58   fSmearType = kSmearEffic;
  59   fEMCALType = kHits;
  60   fTrackType = kCharged;
  61   fEfficiency = 0.90;
  62
  63 }
  64 AliEMCALJetFinderInputSimPrep::~AliEMCALJetFinderInputSimPrep()
  65 {
  66 if (fDebug > 0) Info("~AliEMCALJetFinderInputSimPrep","Beginning Destructor");
  67
  68 }
  69
  70 void AliEMCALJetFinderInputSimPrep::Reset(AliEMCALJetFinderResetType_t resettype)
  71 {
  72 if (fDebug > 1) Info("Reset","Beginning Reset");
  73         switch (resettype){
  74
  75         case kResetData:
  76                 fInputObject.Reset(resettype);
  77                 break;
  78         case kResetTracks:
  79                 fInputObject.Reset(resettype);
  80                 break;
  81         case kResetDigits:
  82                 fInputObject.Reset(resettype);
  83                 break;
  84         case kResetParameters:
  85                 fSmearType = kSmearEffic;
  86                 fEMCALType = kHits;
  87                 fTrackType = kCharged;
  88                 break;
  89         case kResetAll:
  90                 fSmearType = kSmearEffic;
  91                 fEMCALType = kHits;
  92                 fTrackType = kCharged;
  93                 fInputObject.Reset(resettype);
  94                 break;
  95         case kResetPartons:
  96           Warning("FillFromFile", "kResetPartons not implemented") ;
  97           break;
  98         case kResetParticles:
  99           Warning("FillFromFile", "kResetParticles not implemented") ;
 100           break;
 101         case kResetJets:
 102           Warning("FillFromFile", "kResetJets not implemented") ;
 103           break;
 104         }// end switch
 105
 106 }
 107
 108 Int_t AliEMCALJetFinderInputSimPrep::FillFromFile(TString *filename, AliEMCALJetFinderFileType_t filetype,Int_t EventNumber)
 109 {
 110   // gObjectTable->Print();
 111         // Test that file exists - Getter doesn't like bogus filenames
 112 if (fDebug > 1) Info("FillFromFile","Beginning FillFromFile");
 113    fFileType = filetype;
 114    TFile file(filename->Data());
 115    if (!file.IsOpen()){
 116            Error("FillFromFile","Could not open file in FillFromFile");
 117            return -1;
 118    }
 119    file.Close();
 120    /*
 121    gAlice =  static_cast<AliRun *>(file.Get("gAlice"));
 122    AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 123    if (gAlice->GetEvent(EventNumber) < 0){
 124            Error("FillFromFile","Could not open event in FillFromFile");
 125            file.Close();
 126            return -1;
 127    }*/
 128    AliEMCALGetter *gime = AliEMCALGetter::Instance(filename->Data());
 129    gime->Event(EventNumber) ;
 130
 131    if ( fEMCALType == kHits ||
 132         fEMCALType == kTimeCut )  FillHits();
 133    if ( fTrackType != kNoTracks  )  FillTracks();
 134    if ( fFileType  != kData){
 135            FillPartons();
 136 //         FillParticles();
 137    }
 138    //gime->CloseFile();
 139 //   gObjectTable->Print();
 140    delete gime;
 141    return 0;
 142 }
 143
 144 void AliEMCALJetFinderInputSimPrep::FillHits()          // Fill from the hits to input object from simulation
 145 {
 146 if (fDebug > 1) Info("FillHits","Beginning FillHits");
 147
 148 // Access hit information
 149     AliEMCALHit *mHit;
 150     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 151     Info("AliEMCALJetFinderInputSimPrep","Title of the geometry is %s",pEMCAL->GetTitle());
 152     AliEMCALGeometry* geom =  AliEMCALGetter::Instance()->EMCALGeometry();
 153 //    Float_t samplingF = geom->GetSampling();
 154     Float_t samplingF = 11.6;
 155     Info("AliEMCALJetFinderInputSimPrep","Sampling fraction is %f",samplingF);
 156     TTree *treeH = AliEMCALGetter::Instance()->TreeH();
 157     Int_t ntracks = (Int_t) treeH->GetEntries();
 158 //
 159 //   Loop over tracks
 160 //
 161     Int_t nbytes = 0;
 162     Double_t etH = 0.0;
 163
 164     for (Int_t track=0; track<ntracks;track++) {
 165         gAlice->ResetHits();
 166         nbytes += treeH->GetEvent(track);
 167 //
 168 //  Loop over hits
 169 //
 170         for(mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
 171             mHit;
 172             mHit=(AliEMCALHit*) pEMCAL->NextHit())
 173         {
 174                 if (fEMCALType == kTimeCut &&
 175                     (mHit->GetTime()>fTimeCut) ) continue;
 176             Float_t x      =    mHit->X();         // x-pos of hit
 177             Float_t y      =    mHit->Y();         // y-pos
 178             Float_t z      =    mHit->Z();         // z-pos
 179             Float_t eloss  =    mHit->GetEnergy(); // deposited energy
 180
 181             Float_t r      =    TMath::Sqrt(x*x+y*y);
 182             Float_t theta  =    TMath::ATan2(r,z);
 183             Float_t eta    =   -TMath::Log(TMath::Tan(theta/2.));
 184             Float_t phi    =    TMath::ATan2(y,x);
 185             etH = samplingF*eloss*TMath::Sin(theta);
 186             if (fDebug > 10) Info("FillHits","Values of hit energy %i",Int_t(1e7*etH));
 187             if (geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi) == -1)
 188             {
 189                     if (fDebug >5)
 190                     {
 191                             Error("FillHits","Hit not inside EMCAL!");
 192                             mHit->Dump();
 193                     }
 194                     continue;
 195             }
 196             fInputObject.AddEnergyToDigit(geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi)-1,Int_t(1e7*etH));
 197         } // Hit Loop
 198     } // Track Loop
 199
 200
 201 }
 202 void AliEMCALJetFinderInputSimPrep::FillTracks()        // Fill from particles simulating a TPC to input object from simulation
 203 {
 204
 205     if (fDebug > 1) Info("FillTracks","Beginning FillTracks");
 206
 207     TParticlePDG* pdgP = 0;
 208     TParticle *MPart;
 209     Int_t npart = (gAlice->GetHeader())->GetNprimary();
 210     Float_t bfield,rEMCAL;
 211
 212     if (fDebug > 1) Info("FillTracks","Defining the geometry");
 213
 214     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 215     AliEMCALGeometry* geom =  AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
 216     fEtaMax = geom->GetArm1EtaMax();
 217     fEtaMin = geom->GetArm1EtaMin();
 218     fPhiMax = TMath::Pi()/180.0*geom->GetArm1PhiMax();
 219     fPhiMin = TMath::Pi()/180.0*geom->GetArm1PhiMin();
 220
 221     if (fDebug > 1) Info("FillTracks","Starting particle loop");
 222
 223     for (Int_t part = 0; part < npart; part++) {
 224         MPart = gAlice->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 (MPart->Eta() > fEtaMax || MPart->Eta() < fEtaMin)    continue;
 240         if (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() == kPositron ){
 285                       if (fDebug > 5) Info("FillTracks","Storing electron or positron");
 286                       if (fSmearType == kSmear ||
 287                             fSmearType == kSmearEffic ){
 288                               Smear(MPart);/*
 289                               TParticle *tmp = Smear(MPart);
 290                               fInputObject.AddTrack(tmp);
 291                               delete tmp;*/
 292                       }else{
 293                           fInputObject.AddTrack(*MPart);
 294                       }
 295                 }
 296                 if ( MPart->GetPdgCode() == kGamma ){
 297                         fInputObject.AddTrack(*MPart);
 298                         if (fDebug > 5) Info("FillTracks","Storing gamma");
 299                 }
 300
 301            break;
 302            case kCharged: // All Particles with non-zero charge
 303                 if (pdgP->Charge() != 0) {
 304                         if (fDebug > 5) Info("FillTracks","Storing charged track");
 305                         if (fSmearType == kSmear ||
 306                             fSmearType == kSmearEffic ){
 307                                 Smear(MPart);/*
 308                                 TParticle *tmp = Smear(MPart);
 309                                 fInputObject.AddTrack(tmp);
 310                                 delete tmp;*/
 311                         }else{
 312                          fInputObject.AddTrack(*MPart);
 313                         }
 314                 }
 315            break;
 316            case kNeutral: // All particles with no charge
 317                 if (pdgP->Charge() == 0){
 318                         fInputObject.AddTrack(*MPart);
 319                         if (fDebug > 5) Info("FillTracks","Storing neutral");
 320                 }
 321            break;
 322            case kHadron: //All hadrons
 323                 if (MPart->GetPdgCode() != kElectron &&
 324                     MPart->GetPdgCode() != kPositron &&
 325                     MPart->GetPdgCode() != kGamma )
 326                 {
 327                         if (fDebug > 5) Info("FillTracks","Storing hadron");
 328                         if (pdgP->Charge() == 0){
 329                             fInputObject.AddTrack(*MPart);
 330                         }else{
 331                                 if (fSmearType == kSmear ||
 332                                     fSmearType == kSmearEffic ){
 333                                         Smear(MPart);/*
 334                                         TParticle *tmp = Smear(MPart);
 335                                         fInputObject.AddTrack(tmp);
 336                                         delete tmp;*/
 337                                 }else{
 338                                     fInputObject.AddTrack(*MPart);
 339                                 }
 340                         }
 341                 }
 342            break;
 343            case kChargedHadron:  // only charged hadrons
 344                 if (MPart->GetPdgCode() != kElectron &&
 345                     MPart->GetPdgCode() != kPositron &&
 346                     MPart->GetPdgCode() != kGamma    &&
 347                     pdgP->Charge()      != 0       ){
 348                         if (fDebug > 5) Info("FillTracks","Storing charged hadron");
 349                         if (fSmearType == kSmear ||
 350                             fSmearType == kSmearEffic ){
 351                                 Smear(MPart);/*
 352                                 TParticle *tmp = Smear(MPart);
 353                                 fInputObject.AddTrack(tmp);
 354                                 delete tmp;*/
 355                         }else{
 356                                fInputObject.AddTrack(*MPart);
 357                         }
 358                 }
 359            break;
 360            case kNoTracks:
 361            break;
 362            default:
 363            break;
 364            delete MPart;
 365         }       //end of switch
 366 //      Info("FillTracks","After Particle Storage");
 367         //if (part%10) gObjectTable->Print();
 368
 369    }    //end of particle loop
 370    //gObjectTable->Print();
 371 }
 372
 373 void AliEMCALJetFinderInputSimPrep::FillPartons()               // Fill partons to input object from simulation
 374 {
 375 if (fDebug > 1) Info("FillPartons","Beginning FillPartons");
 376
 377   AliGenEventHeader* evHeader = ((AliHeader*)(gAlice->GetHeader()))->GenEventHeader();
 378   Float_t triggerJetValues[4];
 379   AliEMCALParton tempParton;
 380
 381   if (fFileType == kPythia)
 382   {
 383         Int_t ntriggerjets = ((AliGenPythiaEventHeader*)evHeader)->NTriggerJets();
 384         if (fDebug > 1) Info("FillPartons","Number of trigger jets --> %i",ntriggerjets);
 385         for (Int_t jetcounter = 0 ; jetcounter < ntriggerjets; jetcounter++){
 386         ((AliGenPythiaEventHeader*)evHeader)->TriggerJet(jetcounter,triggerJetValues);
 387         TLorentzVector tempLParton;
 388         tempLParton.SetPxPyPzE(triggerJetValues[0],triggerJetValues[1],triggerJetValues[2],triggerJetValues[3]);
 389         tempParton.SetEnergy(tempLParton.Energy());
 390         tempParton.SetEta(tempLParton.Eta());
 391         tempParton.SetPhi(tempLParton.Phi());
 392         FillPartonTracks(&tempParton);
 393         fInputObject.AddParton(&tempParton);
 394         }
 395
 396   }
 397 }
 398
 399 void AliEMCALJetFinderInputSimPrep::FillParticles()             // Fill particles to input object from simulation
 400 {
 401 if (fDebug > 1) Info("FillParticles","Beginning FillParticles");
 402
 403     Int_t npart = (gAlice->GetHeader())->GetNprimary();
 404     TParticlePDG* pdgP = 0;
 405
 406     for (Int_t part = 0; part < npart; part++) {
 407         TParticle *MPart = gAlice->Particle(part);
 408         pdgP = MPart->GetPDG();
 409
 410         if (fDebug > 10) Info("FillParticles","Checking if particle is a primary");
 411
 412         if (fFileType == kPythia) {
 413             if (MPart->GetStatusCode() != 1) continue;
 414         } else if (fFileType == kHijing) {
 415             if (MPart->GetFirstDaughter() >= 0 && MPart->GetFirstDaughter() < npart) continue;
 416         }
 417
 418         if (fDebug > 10) Info("FillParticles","Checking if particle is in acceptance");
 419
 420         if (MPart->Eta() > fEtaMax || MPart->Eta() < fEtaMin)    continue;
 421         if (MPart->Phi() > fPhiMax || MPart->Phi() < fPhiMin )   continue;
 422
 423
 424 /*
 425         {kProton, kProtonBar, kElectron, kPositron,
 426          kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
 427          kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
 428          kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
 429          kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
 430          0,kNuMu,kNuMuBar
 431 */
 432         switch(fTrackType)
 433         {
 434
 435            case kAll:  // All Stable particles to be included
 436                 if (fDebug > 5) Info("FillParticles","Storing particle");
 437                 fInputObject.AddParticle(MPart);
 438            break;
 439            case kEM:   // All Electromagnetic particles
 440                 if (MPart->GetPdgCode() == kElectron ||
 441                     MPart->GetPdgCode() == kPositron ||
 442                     MPart->GetPdgCode() == kGamma){
 443                         if (fDebug > 5) Info("FillParticles","Storing electromagnetic particle");
 444                         fInputObject.AddParticle(MPart);
 445                 }
 446            break;
 447            case kCharged: // All Particles with non-zero charge
 448                 if (pdgP->Charge() != 0) {
 449                         if (fDebug > 5) Info("FillParticles","Storing charged particle");
 450                         fInputObject.AddParticle(MPart);
 451                 }
 452            break;
 453            case kNeutral: // All particles with no charge
 454                 if (pdgP->Charge() == 0){
 455                         if (fDebug > 5) Info("FillParticles","Storing neutral particle");
 456                         fInputObject.AddParticle(MPart);
 457                 }
 458            break;
 459            case kHadron: //All hadrons
 460                 if (MPart->GetPdgCode() != kElectron &&
 461                     MPart->GetPdgCode() != kPositron &&
 462                     MPart->GetPdgCode() != kGamma )
 463                 {
 464
 465                 if (fDebug > 5) Info("FillParticles","Storing hadron");
 466                     fInputObject.AddParticle(MPart);
 467                 }
 468            break;
 469            case kChargedHadron:  // only charged hadrons
 470                 if (MPart->GetPdgCode() != kElectron &&
 471                     MPart->GetPdgCode() != kPositron &&
 472                     MPart->GetPdgCode() != kGamma    &&
 473                     pdgP->Charge()      != 0       ){
 474                 if (fDebug > 5) Info("FillParticles","Storing charged hadron");
 475                        fInputObject.AddParticle(MPart);
 476                 }
 477            break;
 478            case kNoTracks:
 479            break;
 480            default:
 481            break;
 482         }       //end of switch
 483     }// end of loop over particles
 484
 485 }
 486 void AliEMCALJetFinderInputSimPrep::FillDigits()                // Fill digits to input object
 487 {
 488
 489 }
 490
 491 void AliEMCALJetFinderInputSimPrep::Smear(TParticle *particle)
 492 {
 493 if (fDebug > 5) Info("Smear","Beginning Smear");
 494
 495         Float_t tmp = AliEMCALFast::SmearMomentum(1,particle->P());
 496         Float_t tmpt = (tmp/particle->P()) * particle->Pt();
 497         if (fDebug > 10) Info("Smear","Smearing changed P from %f to %f",particle->Pt(),tmpt);
 498         TLorentzVector tmpvector;
 499         tmpvector.SetPtEtaPhiM(tmpt,particle->Eta(), particle->Phi(), particle->GetMass());
 500 // create a new particle with smearing momentum - and no daughter or parent information and the same
 501 // vertex
 502
 503         TParticle tmparticle(particle->GetPdgCode(), 1, 0, 0,0,0,
 504                         tmpvector.Px()  , tmpvector.Py(), tmpvector.Pz(), tmpvector.Energy(),
 505                         particle->Vx(), particle->Vy(), particle->Vz(), particle->T());
 506         fInputObject.AddTrack(tmparticle);
 507         return;
 508 }
 509
 510 void AliEMCALJetFinderInputSimPrep::FillPartonTracks(AliEMCALParton *parton)
 511 {
 512
 513         if (fDebug>1) Info("FillPartonTracks","Beginning FillPartonTracks");
 514         Int_t npart = (gAlice->GetHeader())->GetNprimary();
 515         Int_t ntracks =0;
 516         TParticlePDG *getpdg;
 517         TParticle *tempPart;
 518         for (Int_t part = 0; part < npart; part++)
 519         {
 520                 tempPart = gAlice->Particle(part);
 521                 if (tempPart->GetStatusCode() != 1) continue;
 522                 if (tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
 523                     tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin ){
 524                         if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
 525                         continue;
 526                 }
 527                 getpdg = tempPart->GetPDG();
 528                 if (getpdg->Charge() == 0.0  ) {
 529                         if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
 530                         continue;
 531                 }
 532                 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
 533                 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 ) ntracks++;
 534         }
 535         Float_t *Energy = new Float_t[ntracks];
 536         Float_t *Eta    = new Float_t[ntracks];
 537         Float_t *Phi    = new Float_t[ntracks];
 538         Int_t   *Pdg    = new Int_t[ntracks];
 539         ntracks=0;
 540         for (Int_t part = 0; part < npart; part++)
 541         {
 542                 tempPart = gAlice->Particle(part);
 543                 if (tempPart->GetStatusCode() != 1) continue;
 544                 if (tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
 545                     tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin ){
 546                         if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
 547                         continue;
 548                 }
 549                 if (tempPart->GetStatusCode() != 1) continue;
 550                 getpdg = tempPart->GetPDG();
 551                 if (getpdg->Charge() == 0.0  ) {
 552                         if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
 553                         continue;
 554                 }
 555                 if (   (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
 556                 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 )
 557                 {
 558                         Energy[ntracks] = tempPart->Pt();
 559                         Eta[ntracks] = tempPart->Eta();
 560                         Phi[ntracks] = tempPart->Phi();
 561                         Pdg[ntracks] = tempPart->GetPdgCode();
 562                         ntracks++;
 563                 }
 564         }
 565         parton->SetTrackList(ntracks,Energy,Eta,Phi,Pdg);
 566 }
 567