cfcfc52a54bb89d044f4901b1df3604a0ec80134
[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)
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 kAllP:  // 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 kAllP:  // 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