Modified plots and made jet finder use SDigits
[u/mrichter/AliRoot.git] / 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