PWG/muondep/AccEffTemplates/GenPythia6Diff.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2
   3 /// !!!!!!!!!!!!!!!!!!
   4 ///
   5 ///      WARNING
   6 ///
   7 /// !!!!!!!!!!!!!!!!!!
   8 ///
   9 /// THIS IS A QUICK HACK IN ORDER TO TEST THE PYTHIA6 TUNED FOR DIFFRACTION
  10 ///
  11 /// IT'S A MERE CUT & PASTE OF THE ORIGIGNAL ALIGENPYTHIA.H/.CXX (VERSION AROUND ALIROOT V5-04-REV-18
  12 /// WITH TRIVIAL ADDITIONS TO TREAT 8 TEV PP AS 7 TEV PP ...
  13 ///
  14 /// IT'S NOT THE WAY TO DO IT ON A PERMANENT BASIS...
  15 ///
  16 ///
  17
  18 #include "AliGenMC.h"
  19 #include "AliPythia.h"
  20
  21 class AliPythia;
  22 class TParticle;
  23 class AliGenPythiaEventHeader;
  24 class AliGenEventHeader;
  25 class AliStack;
  26 class AliRunLoader;
  27 class TObjArray;
  28
  29 class AliGenPythia6Diff : public AliGenMC
  30 {
  31 public:
  32
  33   typedef enum {kFlavorSelection, kParentSelection, kHeavyFlavor} StackFillOpt_t;
  34   typedef enum {kCountAll, kCountParents, kCountTrackables} CountMode_t;
  35   typedef enum {kCluster, kCell} JetRecMode_t;
  36
  37   AliGenPythia6Diff();
  38   AliGenPythia6Diff(Int_t npart);
  39   virtual ~AliGenPythia6Diff();
  40   virtual void    Generate();
  41   virtual void    Init();
  42   // Range of events to be printed
  43   virtual void    SetEventListRange(Int_t eventFirst=-1, Int_t eventLast=-1);
  44   // Select process type
  45   virtual void    SetProcess(Process_t proc = kPyCharm) {fProcess = proc;}
  46   virtual void    SetTune(Int_t itune) {fItune = itune;}
  47
  48   // Select structure function
  49   virtual void    SetStrucFunc(StrucFunc_t func =  kCTEQ5L) {fStrucFunc = func;}
  50   // Select pt of hard scattering
  51   virtual void    SetPtHard(Float_t ptmin = 0, Float_t ptmax = 1.e10)
  52         {fPtHardMin = ptmin; fPtHardMax = ptmax; }
  53   // y of hard scattering
  54   virtual void    SetYHard(Float_t ymin = -1.e10, Float_t ymax = 1.e10)
  55         {fYHardMin = ymin; fYHardMax = ymax; }
  56   // Set initial and final state gluon radiation
  57   virtual void    SetGluonRadiation(Int_t iIn, Int_t iFin)
  58         {fGinit = iIn; fGfinal = iFin;}
  59   // Intrinsic kT
  60   virtual void    SetPtKick(Float_t kt = 1.)
  61         {fPtKick = kt;}
  62   // Use the Pythia 6.3 new multiple interations scenario
  63   virtual void    UseNewMultipleInteractionsScenario() {fNewMIS = kTRUE;}
  64   // Switch off heavy flavors
  65   virtual void    SwitchHFOff() {fHFoff = kTRUE;}
  66   // Set centre of mass energy
  67   virtual void    SetEnergyCMS(Float_t energy = 5500) {fEnergyCMS = energy;}
  68   // Treat protons as inside nuclei with mass numbers a1 and a2
  69   virtual void    SetNuclei(Int_t a1, Int_t a2, Int_t pdfset = 0);
  70   // Set colliding nuclei ("p","n",...)
  71   virtual void    SetCollisionSystem(TString projectile, TString target) { fProjectile = projectile; fTarget = target; }
  72   virtual void    SetNuclearPDF(Int_t pdf) {fNucPdf = pdf;}
  73   virtual void    SetUseNuclearPDF(Bool_t val) {fUseNuclearPDF = val;}
  74   virtual void    SetUseLorentzBoost(Bool_t val) {fUseLorentzBoost = val;}
  75   //
  76   // Trigger options
  77   //
  78   // Energy range for jet trigger
  79   virtual void    SetJetEtRange(Float_t etmin = 0., Float_t etmax = 1.e4)
  80         {fEtMinJet = etmin; fEtMaxJet = etmax;}
  81   // Eta range for jet trigger
  82   virtual void    SetJetEtaRange(Float_t etamin = -20., Float_t etamax = 20.)
  83         {fEtaMinJet = etamin; fEtaMaxJet = etamax;}
  84   // Phi range for jet trigger
  85   virtual void    SetJetPhiRange(Float_t phimin = 0., Float_t phimax = 360.)
  86         {fPhiMinJet = TMath::Pi()*phimin/180.; fPhiMaxJet = TMath::Pi()*phimax/180.;}
  87   // Jet reconstruction mode; default is cone algorithm
  88   virtual void    SetJetReconstructionMode(Int_t mode = kCell) {fJetReconstruction = mode;}
  89   // Eta range for gamma trigger
  90   virtual void    SetGammaEtaRange(Float_t etamin = -20., Float_t etamax = 20.)
  91         {fEtaMinGamma = etamin; fEtaMaxGamma = etamax;}
  92   // Phi range for gamma trigger
  93   virtual void    SetGammaPhiRange(Float_t phimin = 0., Float_t phimax = 360.)
  94         {fPhiMinGamma = TMath::Pi()*phimin/180.; fPhiMaxGamma = TMath::Pi()*phimax/180.;}
  95
  96   // Select events with fragmentation photon, decay photon, pi0 or eta going to PHOS or EMCAL and central barrel
  97   virtual Bool_t TriggerOnSelectedParticles(Int_t np);
  98
  99   virtual void  SetCheckPHOS         (Bool_t b) {fCheckPHOS    = b;}
 100   virtual void  SetCheckEMCAL        (Bool_t b) {fCheckEMCAL   = b;}
 101   virtual void  SetCheckBarrel       (Bool_t b) {fCheckBarrel  = b;}
 102
 103   //virtual void  SetElectronInEMCAL   (Bool_t b) {fEleInEMCAL   = b;}
 104   //virtual void  SetPhotonInPHOS      (Bool_t b) {fCheckPHOS    = b; fPhotonInCalo     = b;} // Not in use
 105
 106   virtual void  SetFragPhotonInCalo  (Bool_t b) {                   fFragPhotonInCalo = b;}
 107   virtual void  SetFragPhotonInBarrel(Bool_t b) {fCheckBarrel  = b; fFragPhotonInCalo = b;}
 108   virtual void  SetFragPhotonInEMCAL (Bool_t b) {fCheckEMCAL   = b; fFragPhotonInCalo = b;}
 109   virtual void  SetFragPhotonInPHOS  (Bool_t b) {fCheckPHOS    = b; fFragPhotonInCalo = b;}
 110
 111   virtual void  SetHadronInCalo      (Bool_t b) {                   fHadronInCalo     = b;}
 112   virtual void  SetHadronInBarrel    (Bool_t b) {fCheckBarrel  = b; fHadronInCalo     = b;}
 113   virtual void  SetHadronInEMCAL     (Bool_t b) {fCheckEMCAL   = b; fHadronInCalo     = b;}
 114   virtual void  SetHadronInPHOS      (Bool_t b) {fCheckPHOS    = b; fHadronInCalo     = b;}
 115
 116   virtual void  SetElectronInCalo    (Bool_t b) {                   fEleInCalo        = b;}
 117   virtual void  SetElectronInBarrel  (Bool_t b) {fCheckBarrel  = b; fEleInCalo        = b;}
 118   virtual void  SetElectronInEMCAL   (Bool_t b) {fCheckEMCAL   = b; fEleInCalo        = b;}
 119   virtual void  SetElectronInPHOS    (Bool_t b) {fCheckPHOS    = b; fEleInCalo        = b;}
 120
 121   virtual void  SetDecayPhotonInCalo (Bool_t d)  {fDecayPhotonInCalo = d;}
 122   virtual void  SetDecayPhotonInBarrel(Bool_t d) {fDecayPhotonInCalo = d; fCheckBarrel  = d;}
 123   virtual void  SetDecayPhotonInEMCAL(Bool_t d)  {fDecayPhotonInCalo = d; fCheckEMCAL   = d;}
 124   virtual void  SetDecayPhotonInPHOS (Bool_t d)  {fDecayPhotonInCalo = d; fCheckPHOS    = d;}
 125
 126   virtual void  SetPi0InCalo         (Bool_t b, Bool_t f = kFALSE) {fPi0InCalo = b; fForceNeutralMeson2PhotonDecay = f;}
 127   virtual void  SetPi0InBarrel       (Bool_t b, Bool_t f = kFALSE) {fPi0InCalo = b; fForceNeutralMeson2PhotonDecay = f; fCheckBarrel= b; }
 128   virtual void  SetPi0InEMCAL        (Bool_t b, Bool_t f = kFALSE) {fPi0InCalo = b; fForceNeutralMeson2PhotonDecay = f; fCheckEMCAL = b; }
 129   virtual void  SetPi0InPHOS         (Bool_t b, Bool_t f = kFALSE) {fPi0InCalo = b; fForceNeutralMeson2PhotonDecay = f; fCheckPHOS  = b; }
 130
 131   virtual void  SetEtaInCalo         (Bool_t b, Bool_t f = kFALSE) {fEtaInCalo = b; fForceNeutralMeson2PhotonDecay = f;}
 132   virtual void  SetEtaInBarrel       (Bool_t b, Bool_t f = kFALSE) {fEtaInCalo = b; fForceNeutralMeson2PhotonDecay = f; fCheckBarrel= b; }
 133   virtual void  SetEtaInEMCAL        (Bool_t b, Bool_t f = kFALSE) {fEtaInCalo = b; fForceNeutralMeson2PhotonDecay = f; fCheckEMCAL = b; }
 134   virtual void  SetEtaInPHOS         (Bool_t b, Bool_t f = kFALSE) {fEtaInCalo = b; fForceNeutralMeson2PhotonDecay = f; fCheckPHOS  = b; }
 135
 136   virtual void  SetPi0PhotonDecayInBarrel(Bool_t b, Bool_t f = kFALSE) {fPi0InCalo = b; fDecayPhotonInCalo = b; fForceNeutralMeson2PhotonDecay = f; fCheckBarrel  = b; }
 137   virtual void  SetPi0PhotonDecayInEMCAL (Bool_t b, Bool_t f = kFALSE) {fPi0InCalo = b; fDecayPhotonInCalo = b; fForceNeutralMeson2PhotonDecay = f; fCheckEMCAL   = b; }
 138   virtual void  SetPi0PhotonDecayInPHOS  (Bool_t b, Bool_t f = kFALSE) {fPi0InCalo = b; fDecayPhotonInCalo = b; fForceNeutralMeson2PhotonDecay = f; fCheckPHOS    = b; }
 139
 140   virtual void  SetEtaPhotonDecayInBarrel(Bool_t b, Bool_t f = kFALSE) {fEtaInCalo = b; fDecayPhotonInCalo = b; fForceNeutralMeson2PhotonDecay = f; fCheckBarrel  = b; }
 141   virtual void  SetEtaPhotonDecayInEMCAL (Bool_t b, Bool_t f = kFALSE) {fEtaInCalo = b; fDecayPhotonInCalo = b; fForceNeutralMeson2PhotonDecay = f; fCheckEMCAL   = b; }
 142   virtual void  SetEtaPhotonDecayInPHOS  (Bool_t b, Bool_t f = kFALSE) {fEtaInCalo = b; fDecayPhotonInCalo = b; fForceNeutralMeson2PhotonDecay = f; fCheckPHOS    = b; }
 143
 144
 145   // Trigger on a minimum multiplicity
 146   virtual void  SetTriggerChargedMultiplicity(Int_t multiplicity, Float_t etamax = 0, Float_t ptmin = -1.)
 147   {fTriggerMultiplicity = multiplicity; fTriggerMultiplicityEta = etamax;
 148     fTriggerMultiplicityPtMin = ptmin;}
 149
 150   // Calorimeters acceptance
 151   // Set Phi in degrees, and Eta coverage, should not be negative
 152   virtual void  SetBarrelAcceptance(Float_t deta) {fTriggerEta = deta ;}
 153   virtual void  SetEMCALAcceptance (Float_t phimin, Float_t phimax, Float_t deta) {fEMCALMinPhi = phimin ; fEMCALMaxPhi = phimax ; fEMCALEta = deta ; }
 154   virtual void  SetPHOSAcceptance  (Float_t phimin, Float_t phimax, Float_t deta) {fPHOSMinPhi  = phimin ; fPHOSMaxPhi  = phimax ; fPHOSEta  = deta ; }
 155   virtual void  SetRotateParticleInPHOSeta(Bool_t b) {fCheckPHOSeta = b;}
 156
 157   virtual void  SetTriggerParticleMinPt(Float_t pt) {fTriggerParticleMinPt = pt;}
 158   //    virtual void  SetPhotonMinPt(Float_t pt)          {fPhotonMinPt = pt;}
 159   //    virtual void  SetElectronMinPt(Float_t pt)        {fElectronMinPt = pt;}
 160   // Trigger and rotate event
 161   void RotatePhi(Bool_t& okdd);
 162
 163   // Trigger on a single particle (not related to calorimeter trigger above)
 164   virtual void    SetTriggerParticle(Int_t particle = 0, Float_t etamax = 0.9, Float_t ptmin = -1, Float_t ptmax = 1000)
 165         {fTriggerParticle = particle; fTriggerEta = etamax; fTriggerMinPt = ptmin; fTriggerMaxPt = ptmax;}
 166
 167   //
 168   // Heavy flavor options
 169   //
 170   // Set option for feed down from higher family
 171   virtual void SetFeedDownHigherFamily(Bool_t opt) {
 172     fFeedDownOpt = opt;
 173   }
 174   // Set option for selecting particles kept in stack according to flavor
 175   // or to parent selection
 176   virtual void SetStackFillOpt(StackFillOpt_t opt) {
 177     fStackFillOpt = opt;
 178   }
 179   // Set fragmentation option
 180   virtual void SetFragmentation(Bool_t opt) {
 181     fFragmentation = opt;
 182   }
 183   // Set counting mode
 184   virtual void SetCountMode(CountMode_t mode) {
 185     fCountMode = mode;
 186   }
 187   //
 188   // Quenching
 189   //
 190   // Set quenching mode 0 = no, 1 = AM, 2 = IL,  3 = NA, 4 = ACS
 191   virtual void SetQuench(Int_t flag = 0) {fQuench = flag;}
 192   // Set transport coefficient.
 193   void SetQhat(Float_t qhat) {fQhat = qhat;}
 194   //Set initial medium length.
 195   void SetLength(Float_t length) {fLength = length;}
 196   //set parameters for pyquen afterburner
 197   virtual void SetPyquenPar(Float_t t0=1., Float_t tau0=0.1, Int_t nf=0,Int_t iengl=0, Int_t iangl=3)
 198   {fpyquenT = t0; fpyquenTau = tau0; fpyquenNf=nf;fpyquenEloss=iengl;fpyquenAngle=iangl;}
 199   virtual void SetHadronisation(Int_t flag = 1) {fHadronisation = flag;}
 200   virtual void SetPatchOmegaDalitz(Int_t flag = 1) {fPatchOmegaDalitz = flag;}
 201   virtual void SetReadFromFile(const Text_t *filname) {fkFileName = filname;  fReadFromFile = 1;}
 202   virtual void SetReadLHEF(const Text_t *filename) {fkNameLHEF = filename; fReadLHEF = 1;}
 203
 204   //
 205   // Pile-up
 206   //
 207   // Get interaction rate for pileup studies
 208   virtual void    SetInteractionRate(Float_t rate,Float_t timewindow = 90.e-6);
 209   virtual Float_t GetInteractionRate() const {return fInteractionRate;}
 210   // get cross section of process
 211   virtual Float_t GetXsection() const {return fXsection;}
 212   // get triggered jets
 213   void GetJets(Int_t& njets, Int_t& ntrig, Float_t jets[4][10]);
 214 //  void RecJetsUA1(Int_t& njets, Float_t jets[4][50]);
 215   void SetPycellParameters(Float_t etamax = 2., Int_t neta = 274, Int_t nphi = 432,
 216                            Float_t thresh = 0., Float_t etseed = 4.,
 217                            Float_t minet = 10., Float_t r = 1.);
 218
 219   void LoadEvent(AliStack* stack, Int_t flag = 0, Int_t reHadr = 0);
 220   void LoadEvent(const TObjArray* stack, Int_t flag = 0, Int_t reHadr = 0);
 221   // Getters
 222   virtual Process_t    GetProcess() const {return fProcess;}
 223   virtual StrucFunc_t  GetStrucFunc() const {return fStrucFunc;}
 224   virtual void         GetPtHard(Float_t& ptmin, Float_t& ptmax) const
 225         {ptmin = fPtHardMin; ptmax = fPtHardMax;}
 226   virtual void         GetNuclei(Int_t&  a1, Int_t& a2) const
 227         {a1 = fAProjectile; a2 = fATarget;}
 228   virtual void         GetJetEtRange(Float_t& etamin, Float_t& etamax) const
 229         {etamin = fEtaMinJet; etamax = fEtaMaxJet;}
 230   virtual void         GetJetPhiRange(Float_t& phimin, Float_t& phimax) const
 231         {phimin = fPhiMinJet*180./TMath::Pi(); phimax = fPhiMaxJet*180/TMath::Pi();}
 232   virtual void         GetGammaEtaRange(Float_t& etamin, Float_t& etamax) const
 233         {etamin = fEtaMinGamma; etamax = fEtaMaxGamma;}
 234   virtual void         GetGammaPhiRange(Float_t& phimin, Float_t& phimax) const
 235         {phimin = fPhiMinGamma*180./TMath::Pi(); phimax = fPhiMaxGamma*180./TMath::Pi();}
 236   //
 237   Bool_t CheckDetectorAcceptance(Float_t phi, Float_t eta, Int_t iparticle);
 238   Bool_t IsInEMCAL (Float_t phi, Float_t eta) const;
 239   Bool_t IsInPHOS  (Float_t phi, Float_t eta, Int_t iparticle) ;
 240   Bool_t IsInBarrel(Float_t eta) const;
 241   Bool_t IsFromHeavyFlavor(Int_t ipart);
 242   //
 243   virtual void FinishRun();
 244   Bool_t CheckTrigger(const TParticle* jet1, const TParticle* jet2);
 245   //Used in some processes to selected child properties
 246   Bool_t CheckKinematicsOnChild();
 247   void     GetSubEventTime();
 248
 249   void SetTuneForDiff(Bool_t a=kTRUE) {fkTuneForDiff=a;}
 250   AliDecayer * GetDecayer(){return fDecayer;}
 251
 252 protected:
 253   // adjust the weight from kinematic cuts
 254   void     AdjustWeights() const;
 255   Int_t    GenerateMB();
 256   void     MakeHeader();
 257   void     GeneratePileup();
 258   Process_t   fProcess;           //Process type
 259   Int_t       fItune;             // Pythia tune > 6.4
 260   StrucFunc_t fStrucFunc;         //Structure Function
 261   Float_t     fKineBias;          //!Bias from kinematic selection
 262   Int_t       fTrials;            //!Number of trials for current event
 263   Int_t       fTrialsRun;         //!Number of trials for run
 264   Float_t     fQ;                 //Mean Q
 265   Float_t     fX1;                //Mean x1
 266   Float_t     fX2;                //Mean x2
 267   Float_t     fEventTime;         //Time of the subevent
 268   Float_t     fInteractionRate;   //Interaction rate (set by user)
 269   Float_t     fTimeWindow;        //Time window for pileup events (set by user)
 270   Int_t       fCurSubEvent;       //Index of the current sub-event
 271   TArrayF     *fEventsTime;       //Subevents time for pileup
 272   Int_t       fNev;               //Number of events
 273   Int_t       fFlavorSelect;      //Heavy Flavor Selection
 274   Float_t     fXsection;          //Cross-section
 275   AliPythia   *fPythia;           //!Pythia
 276   Float_t     fPtHardMin;         //lower pT-hard cut
 277   Float_t     fPtHardMax;         //higher pT-hard cut
 278   Float_t     fYHardMin;          //lower  y-hard cut
 279   Float_t     fYHardMax;          //higher y-hard cut
 280   Int_t       fGinit;             //initial state gluon radiation
 281   Int_t       fGfinal;            //final state gluon radiation
 282   Int_t       fHadronisation;     //hadronisation
 283   Bool_t      fPatchOmegaDalitz;  //flag for omega dalitz decay patch
 284   Int_t       fNpartons;          //Number of partons before hadronisation
 285   Int_t       fReadFromFile;      //read partons from file
 286   Int_t       fReadLHEF;          //read lhef file
 287   Int_t       fQuench;            //Flag for quenching
 288   Float_t     fQhat;              //Transport coefficient (GeV^2/fm)
 289   Float_t     fLength;            //Medium length (fm)
 290   Float_t     fpyquenT;           //Pyquen initial temperature
 291   Float_t     fpyquenTau;         //Pyquen initial proper time
 292   Int_t       fpyquenNf;          //Pyquen number of flavours into the game
 293   Int_t       fpyquenEloss;       //Pyquen type of energy loss
 294   Int_t       fpyquenAngle;       //Pyquen radiation angle for gluons
 295   Float_t     fImpact;            //Impact parameter for quenching simulation (q-pythia)
 296   Float_t     fPtKick;            //Transverse momentum kick
 297   Bool_t      fFullEvent;         //!Write Full event if true
 298   AliDecayer  *fDecayer;          //!Pointer to the decayer instance
 299   Int_t       fDebugEventFirst;   //!First event to debug
 300   Int_t       fDebugEventLast;    //!Last  event to debug
 301   Float_t     fEtMinJet;          //Minimum et of triggered Jet
 302   Float_t     fEtMaxJet;          //Maximum et of triggered Jet
 303   Float_t     fEtaMinJet;         //Minimum eta of triggered Jet
 304   Float_t     fEtaMaxJet;         //Maximum eta of triggered Jet
 305   Float_t     fPhiMinJet;         //Minimum phi of triggered Jet
 306   Float_t     fPhiMaxJet;         //Maximum phi of triggered Jet
 307   Int_t       fJetReconstruction; //Jet Reconstruction mode
 308   Float_t     fEtaMinGamma;       // Minimum eta of triggered gamma
 309   Float_t     fEtaMaxGamma;       // Maximum eta of triggered gamma
 310   Float_t     fPhiMinGamma;       // Minimum phi of triggered gamma
 311   Float_t     fPhiMaxGamma;       // Maximum phi of triggered gamma
 312   Float_t     fPycellEtaMax;      // Max. eta for Pycell
 313   Int_t       fPycellNEta;        // Number of eta bins for Pycell
 314   Int_t       fPycellNPhi;        // Number of phi bins for Pycell
 315   Float_t       fPycellThreshold;   // Pycell threshold
 316   Float_t       fPycellEtSeed;      // Pycell seed
 317   Float_t       fPycellMinEtJet;    // Pycell min. jet et
 318   Float_t       fPycellMaxRadius;   // Pycell cone radius
 319   StackFillOpt_t fStackFillOpt;   // Stack filling with all particles with
 320   // that flavour or only with selected
 321   // parents and their decays
 322   Bool_t fFeedDownOpt;            // Option to set feed down from higher
 323   // quark families (e.g. b->c)
 324   Bool_t  fFragmentation;         // Option to activate fragmentation by Pythia
 325   Bool_t  fSetNuclei;             // Flag indicating that SetNuclei has been called
 326   Bool_t  fUseNuclearPDF;         // flag if nuclear pdf should be applied
 327   Bool_t  fUseLorentzBoost;       // flag if lorentz boost should be applied
 328   Bool_t  fNewMIS;                // Flag for the new multipple interactions scenario
 329   Bool_t  fHFoff;                 // Flag for switching heafy flavor production off
 330   Int_t   fNucPdf;                // Nuclear pdf 0: EKS98 1: EPS08
 331   Int_t   fTriggerParticle;       // Trigger on this particle ...
 332   Float_t fTriggerEta;            // .. within |eta| < fTriggerEta
 333   Float_t fTriggerMinPt;          // .. within pt > fTriggerMinPt
 334   Float_t fTriggerMaxPt;          // .. within pt < fTriggerMaxPt
 335   Int_t       fTriggerMultiplicity;       // Trigger on events with a minimum charged multiplicity
 336   Float_t     fTriggerMultiplicityEta;    // in a given eta range
 337   Float_t     fTriggerMultiplicityPtMin;  // above this pT
 338   CountMode_t fCountMode;         // Options for counting when the event will be finished.
 339   // fCountMode = kCountAll         --> All particles that end up in the
 340   //                                    stack are counted
 341   // fCountMode = kCountParents     --> Only selected parents are counted
 342   // fCountMode = kCountTrackabless --> Only particles flagged for tracking
 343   //                                     are counted
 344   //
 345   //
 346
 347   AliGenPythiaEventHeader* fHeader;  //! Event header
 348   AliRunLoader*            fRL;      //! Run Loader
 349   const Text_t* fkFileName;          //! Name of file to read from
 350   const Text_t* fkNameLHEF;          //! Name of lhef file to read from
 351   Bool_t fFragPhotonInCalo; // Option to ask for Fragmentation Photon in calorimeters acceptance
 352   Bool_t fHadronInCalo;     // Option to ask for hadron (not pi0) in calorimeters acceptance
 353   Bool_t fPi0InCalo;        // Option to ask for Pi0 in calorimeters acceptance
 354   Bool_t fEtaInCalo;        // Option to ask for Eta in calorimeters acceptance
 355   Bool_t fPhotonInCalo;     // Option to ask for Photon in calorimeter acceptance (not in use)
 356   Bool_t fDecayPhotonInCalo;// Option to ask for Decay Photon in calorimeter acceptance
 357   Bool_t fForceNeutralMeson2PhotonDecay; // Option to ask for Pi0/Eta in calorimeters acceptance when decay into 2 photons
 358   Bool_t fEleInCalo;        // Option to ask for Electron in EMCAL acceptance
 359   Bool_t fEleInEMCAL;       // Option to ask for Electron in EMCAL acceptance (not in use)
 360   Bool_t fCheckBarrel;      // Option to ask for FragPhoton or Pi0 or Eta or gamma decays in central barrel acceptance
 361   Bool_t fCheckEMCAL;       // Option to ask for FragPhoton or Pi0 or Eta or gamma decays in calorimeters EMCAL acceptance
 362   Bool_t fCheckPHOS;        // Option to ask for FragPhoton or Pi0 or Eta or gamma decays in calorimeters PHOS acceptance
 363   Bool_t fCheckPHOSeta;     // Option to ask for rotate event particles in phi to have in PHOS acceptance a requested particle that previously had the good eta
 364   Int_t  fPHOSRotateCandidate; // Internal member to select the particle candidate to trigger the event phi rotation, to put it in PHOS phi acceptance
 365   Float_t fTriggerParticleMinPt; // Minimum momentum of Fragmentation Photon or Pi0 or other hadron
 366   Float_t fPhotonMinPt;          // Minimum momentum of Photon  (not in use)
 367   Float_t fElectronMinPt;        // Minimum momentum of Electron (not in use)
 368   //Calorimeters eta-phi acceptance
 369   Float_t fPHOSMinPhi;           // Minimum phi PHOS, degrees
 370   Float_t fPHOSMaxPhi;           // Maximum phi PHOS, degrees
 371   Float_t fPHOSEta;              // Minimum eta PHOS, coverage delta eta
 372   Float_t fEMCALMinPhi;          // Minimum phi EMCAL, degrees
 373   Float_t fEMCALMaxPhi;          // Maximum phi EMCAL, degrees
 374   Float_t fEMCALEta;             // Maximum eta EMCAL, coverage delta eta
 375
 376   Bool_t fkTuneForDiff;    // Pythia tune
 377   Int_t  fProcDiff;
 378 private:
 379   AliGenPythia6Diff(const AliGenPythia6Diff &Pythia);
 380   AliGenPythia6Diff & operator=(const AliGenPythia6Diff & rhs);
 381
 382
 383   Bool_t CheckDiffraction();
 384   Bool_t GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax,
 385                                Double_t &wSD, Double_t &wDD, Double_t &wND);
 386
 387   ClassDef(AliGenPythia6Diff, 14) // AliGenerator interface to Pythia
 388 };
 389
 390 #endif
 391
 392 #include <TClonesArray.h>
 393 #include <TDatabasePDG.h>
 394 #include <TParticle.h>
 395 #include <TPDGCode.h>
 396 #include <TObjArray.h>
 397 #include <TSystem.h>
 398 #include <TTree.h>
 399 #include "AliConst.h"
 400 #include "AliDecayerPythia.h"
 401 #include "AliFastGlauber.h"
 402 #include "AliHeader.h"
 403 #include "AliGenPythiaEventHeader.h"
 404 #include "AliPythia.h"
 405 #include "AliPythiaRndm.h"
 406 #include "AliRun.h"
 407 #include "AliStack.h"
 408 #include "AliRunLoader.h"
 409 #include "AliMC.h"
 410 #include "AliLog.h"
 411 #include "PyquenCommon.h"
 412
 413 ClassImp(AliGenPythia6Diff)
 414
 415
 416 AliGenPythia6Diff::AliGenPythia6Diff():
 417 AliGenMC(),
 418 fProcess(kPyCharm),
 419 fItune(-1),
 420 fStrucFunc(kCTEQ5L),
 421 fKineBias(0.),
 422 fTrials(0),
 423 fTrialsRun(0),
 424 fQ(0.),
 425 fX1(0.),
 426 fX2(0.),
 427 fEventTime(0.),
 428 fInteractionRate(0.),
 429 fTimeWindow(0.),
 430 fCurSubEvent(0),
 431 fEventsTime(0),
 432 fNev(0),
 433 fFlavorSelect(0),
 434 fXsection(0.),
 435 fPythia(0),
 436 fPtHardMin(0.),
 437 fPtHardMax(1.e4),
 438 fYHardMin(-1.e10),
 439 fYHardMax(1.e10),
 440 fGinit(1),
 441 fGfinal(1),
 442 fHadronisation(1),
 443 fPatchOmegaDalitz(0),
 444 fNpartons(0),
 445 fReadFromFile(0),
 446 fReadLHEF(0),
 447 fQuench(0),
 448 fQhat(0.),
 449 fLength(0.),
 450 fpyquenT(1.),
 451 fpyquenTau(0.1),
 452 fpyquenNf(0),
 453 fpyquenEloss(0),
 454 fpyquenAngle(0),
 455 fImpact(0.),
 456 fPtKick(1.),
 457 fFullEvent(kTRUE),
 458 fDecayer(new AliDecayerPythia()),
 459 fDebugEventFirst(-1),
 460 fDebugEventLast(-1),
 461 fEtMinJet(0.),
 462 fEtMaxJet(1.e4),
 463 fEtaMinJet(-20.),
 464 fEtaMaxJet(20.),
 465 fPhiMinJet(0.),
 466 fPhiMaxJet(2.* TMath::Pi()),
 467 fJetReconstruction(kCell),
 468 fEtaMinGamma(-20.),
 469 fEtaMaxGamma(20.),
 470 fPhiMinGamma(0.),
 471 fPhiMaxGamma(2. * TMath::Pi()),
 472 fPycellEtaMax(2.),
 473 fPycellNEta(274),
 474 fPycellNPhi(432),
 475 fPycellThreshold(0.),
 476 fPycellEtSeed(4.),
 477 fPycellMinEtJet(10.),
 478 fPycellMaxRadius(1.),
 479 fStackFillOpt(kFlavorSelection),
 480 fFeedDownOpt(kTRUE),
 481 fFragmentation(kTRUE),
 482 fSetNuclei(kFALSE),
 483 fUseNuclearPDF(kFALSE),
 484 fUseLorentzBoost(kTRUE),
 485 fNewMIS(kFALSE),
 486 fHFoff(kFALSE),
 487 fNucPdf(0),
 488 fTriggerParticle(0),
 489 fTriggerEta(0.9),
 490 fTriggerMinPt(-1),
 491 fTriggerMaxPt(1000),
 492 fTriggerMultiplicity(0),
 493 fTriggerMultiplicityEta(0),
 494 fTriggerMultiplicityPtMin(0),
 495 fCountMode(kCountAll),
 496 fHeader(0),
 497 fRL(0),
 498 fkFileName(0),
 499 fkNameLHEF(0),
 500 fFragPhotonInCalo(kFALSE),
 501 fHadronInCalo(kFALSE) ,
 502 fPi0InCalo(kFALSE) ,
 503 fEtaInCalo(kFALSE) ,
 504 fPhotonInCalo(kFALSE), // not in use
 505 fDecayPhotonInCalo(kFALSE),
 506 fForceNeutralMeson2PhotonDecay(kFALSE),
 507 fEleInCalo(kFALSE),
 508 fEleInEMCAL(kFALSE), // not in use
 509 fCheckBarrel(kFALSE),
 510 fCheckEMCAL(kFALSE),
 511 fCheckPHOS(kFALSE),
 512 fCheckPHOSeta(kFALSE),
 513 fPHOSRotateCandidate(-1),
 514 fTriggerParticleMinPt(0),
 515 fPhotonMinPt(0), // not in use
 516 fElectronMinPt(0), // not in use
 517 fPHOSMinPhi(219.),
 518 fPHOSMaxPhi(321.),
 519 fPHOSEta(0.13),
 520 fEMCALMinPhi(79.),
 521 fEMCALMaxPhi(191.),
 522 fEMCALEta(0.71),
 523 fkTuneForDiff(0),
 524 fProcDiff(0)
 525 {
 526   // Default Constructor
 527   fEnergyCMS = 5500.;
 528   if (!AliPythiaRndm::GetPythiaRandom())
 529     AliPythiaRndm::SetPythiaRandom(GetRandom());
 530 }
 531
 532 AliGenPythia6Diff::AliGenPythia6Diff(Int_t npart)
 533 :AliGenMC(npart),
 534 fProcess(kPyCharm),
 535 fItune(-1),
 536 fStrucFunc(kCTEQ5L),
 537 fKineBias(0.),
 538 fTrials(0),
 539 fTrialsRun(0),
 540 fQ(0.),
 541 fX1(0.),
 542 fX2(0.),
 543 fEventTime(0.),
 544 fInteractionRate(0.),
 545 fTimeWindow(0.),
 546 fCurSubEvent(0),
 547 fEventsTime(0),
 548 fNev(0),
 549 fFlavorSelect(0),
 550 fXsection(0.),
 551 fPythia(0),
 552 fPtHardMin(0.),
 553 fPtHardMax(1.e4),
 554 fYHardMin(-1.e10),
 555 fYHardMax(1.e10),
 556 fGinit(kTRUE),
 557 fGfinal(kTRUE),
 558 fHadronisation(kTRUE),
 559 fPatchOmegaDalitz(0),
 560 fNpartons(0),
 561 fReadFromFile(kFALSE),
 562 fReadLHEF(0),
 563 fQuench(kFALSE),
 564 fQhat(0.),
 565 fLength(0.),
 566 fpyquenT(1.),
 567 fpyquenTau(0.1),
 568 fpyquenNf(0),
 569 fpyquenEloss(0),
 570 fpyquenAngle(0),
 571 fImpact(0.),
 572 fPtKick(1.),
 573 fFullEvent(kTRUE),
 574 fDecayer(new AliDecayerPythia()),
 575 fDebugEventFirst(-1),
 576 fDebugEventLast(-1),
 577 fEtMinJet(0.),
 578 fEtMaxJet(1.e4),
 579 fEtaMinJet(-20.),
 580 fEtaMaxJet(20.),
 581 fPhiMinJet(0.),
 582 fPhiMaxJet(2.* TMath::Pi()),
 583 fJetReconstruction(kCell),
 584 fEtaMinGamma(-20.),
 585 fEtaMaxGamma(20.),
 586 fPhiMinGamma(0.),
 587 fPhiMaxGamma(2. * TMath::Pi()),
 588 fPycellEtaMax(2.),
 589 fPycellNEta(274),
 590 fPycellNPhi(432),
 591 fPycellThreshold(0.),
 592 fPycellEtSeed(4.),
 593 fPycellMinEtJet(10.),
 594 fPycellMaxRadius(1.),
 595 fStackFillOpt(kFlavorSelection),
 596 fFeedDownOpt(kTRUE),
 597 fFragmentation(kTRUE),
 598 fSetNuclei(kFALSE),
 599 fUseNuclearPDF(kFALSE),
 600 fUseLorentzBoost(kTRUE),
 601 fNewMIS(kFALSE),
 602 fHFoff(kFALSE),
 603 fNucPdf(0),
 604 fTriggerParticle(0),
 605 fTriggerEta(0.9),
 606 fTriggerMinPt(-1),
 607 fTriggerMaxPt(1000),
 608 fTriggerMultiplicity(0),
 609 fTriggerMultiplicityEta(0),
 610 fTriggerMultiplicityPtMin(0),
 611 fCountMode(kCountAll),
 612 fHeader(0),
 613 fRL(0),
 614 fkFileName(0),
 615 fkNameLHEF(0),
 616 fFragPhotonInCalo(kFALSE),
 617 fHadronInCalo(kFALSE) ,
 618 fPi0InCalo(kFALSE) ,
 619 fEtaInCalo(kFALSE) ,
 620 fPhotonInCalo(kFALSE), // not in use
 621 fDecayPhotonInCalo(kFALSE),
 622 fForceNeutralMeson2PhotonDecay(kFALSE),
 623 fEleInCalo(kFALSE),
 624 fEleInEMCAL(kFALSE), // not in use
 625 fCheckBarrel(kFALSE),
 626 fCheckEMCAL(kFALSE),
 627 fCheckPHOS(kFALSE),
 628 fCheckPHOSeta(kFALSE),
 629 fPHOSRotateCandidate(-1),
 630 fTriggerParticleMinPt(0),
 631 fPhotonMinPt(0), // not in use
 632 fElectronMinPt(0), // not in use
 633 fPHOSMinPhi(219.),
 634 fPHOSMaxPhi(321.),
 635 fPHOSEta(0.13),
 636 fEMCALMinPhi(79.),
 637 fEMCALMaxPhi(191.),
 638 fEMCALEta(0.71),
 639 fkTuneForDiff(0),
 640 fProcDiff(0)
 641 {
 642   // default charm production at 5. 5 TeV
 643   // semimuonic decay
 644   // structure function GRVHO
 645   //
 646   fEnergyCMS = 5500.;
 647   fName = "Pythia";
 648   fTitle= "Particle Generator using PYTHIA";
 649   SetForceDecay();
 650   // Set random number generator
 651   if (!AliPythiaRndm::GetPythiaRandom())
 652     AliPythiaRndm::SetPythiaRandom(GetRandom());
 653 }
 654
 655 AliGenPythia6Diff::~AliGenPythia6Diff()
 656 {
 657   // Destructor
 658   if(fEventsTime) delete fEventsTime;
 659 }
 660
 661 void AliGenPythia6Diff::SetInteractionRate(Float_t rate,Float_t timewindow)
 662 {
 663   // Generate pileup using user specified rate
 664   fInteractionRate = rate;
 665   fTimeWindow = timewindow;
 666   GeneratePileup();
 667 }
 668
 669 void AliGenPythia6Diff::GeneratePileup()
 670 {
 671   // Generate sub events time for pileup
 672   fEventsTime = 0;
 673   if(fInteractionRate == 0.) {
 674     Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
 675     return;
 676   }
 677
 678   Int_t npart = NumberParticles();
 679   if(npart < 0) {
 680     Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
 681     return;
 682   }
 683
 684   if(fEventsTime) delete fEventsTime;
 685   fEventsTime = new TArrayF(npart);
 686   TArrayF &array = *fEventsTime;
 687   for(Int_t ipart = 0; ipart < npart; ipart++)
 688     array[ipart] = 0.;
 689
 690   Float_t eventtime = 0.;
 691   while(1)
 692   {
 693     eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
 694     if(eventtime > fTimeWindow) break;
 695     array.Set(array.GetSize()+1);
 696     array[array.GetSize()-1] = eventtime;
 697   }
 698
 699   eventtime = 0.;
 700   while(1)
 701   {
 702     eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
 703     if(TMath::Abs(eventtime) > fTimeWindow) break;
 704     array.Set(array.GetSize()+1);
 705     array[array.GetSize()-1] = eventtime;
 706   }
 707
 708   SetNumberParticles(fEventsTime->GetSize());
 709 }
 710
 711 void AliGenPythia6Diff::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
 712                                             Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
 713 {
 714   // Set pycell parameters
 715   fPycellEtaMax    =  etamax;
 716   fPycellNEta      =  neta;
 717   fPycellNPhi      =  nphi;
 718   fPycellThreshold =  thresh;
 719   fPycellEtSeed    =  etseed;
 720   fPycellMinEtJet  =  minet;
 721   fPycellMaxRadius =  r;
 722 }
 723
 724
 725
 726 void AliGenPythia6Diff::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 727 {
 728   // Set a range of event numbers, for which a table
 729   // of generated particle will be printed
 730   fDebugEventFirst = eventFirst;
 731   fDebugEventLast  = eventLast;
 732   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
 733 }
 734
 735 void AliGenPythia6Diff::Init()
 736 {
 737   // Initialisation
 738
 739   SetMC(AliPythia::Instance());
 740   fPythia=(AliPythia*) fMCEvGen;
 741
 742   //
 743   fParentWeight=1./Float_t(fNpart);
 744   //
 745
 746
 747   fPythia->SetCKIN(3,fPtHardMin);
 748   fPythia->SetCKIN(4,fPtHardMax);
 749   fPythia->SetCKIN(7,fYHardMin);
 750   fPythia->SetCKIN(8,fYHardMax);
 751
 752   if (fProjectile != "p" || fTarget != "p") fPythia->SetCollisionSystem(fProjectile,fTarget);
 753
 754   if(fUseNuclearPDF)
 755     fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
 756   // Fragmentation?
 757   if (fFragmentation) {
 758     fPythia->SetMSTP(111,1);
 759   } else {
 760     fPythia->SetMSTP(111,0);
 761   }
 762
 763
 764   //  initial state radiation
 765   fPythia->SetMSTP(61,fGinit);
 766   //  final state radiation
 767   fPythia->SetMSTP(71,fGfinal);
 768   //  pt - kick
 769   if (fPtKick > 0.) {
 770     fPythia->SetMSTP(91,1);
 771     fPythia->SetPARP(91,fPtKick);
 772     fPythia->SetPARP(93, 4. * fPtKick);
 773   } else {
 774     fPythia->SetMSTP(91,0);
 775   }
 776
 777   if (fReadLHEF) fPythia->OpenFortranFile(97, const_cast<char*>(fkNameLHEF));
 778
 779   if (fReadFromFile) {
 780     fRL  =  AliRunLoader::Open(fkFileName, "Partons");
 781     fRL->LoadKinematics();
 782     fRL->LoadHeader();
 783   } else {
 784     fRL = 0x0;
 785   }
 786   //
 787   fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
 788   //  Forward Paramters to the AliPythia object
 789   fDecayer->SetForceDecay(fForceDecay);
 790   // Switch off Heavy Flavors on request
 791   if (fHFoff) {
 792     // Maximum number of quark flavours used in pdf
 793     fPythia->SetMSTP(58, 3);
 794     // Maximum number of flavors that can be used in showers
 795     fPythia->SetMSTJ(45, 3);
 796     // Switch off g->QQbar splitting in decay table
 797     ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
 798   }
 799
 800   fDecayer->Init();
 801
 802
 803   //  Parent and Children Selection
 804   switch (fProcess)
 805   {
 806     case kPyOldUEQ2ordered:
 807     case kPyOldUEQ2ordered2:
 808     case kPyOldPopcorn:
 809       break;
 810     case kPyCharm:
 811     case kPyCharmUnforced:
 812     case kPyCharmPbPbMNR:
 813     case kPyCharmpPbMNR:
 814     case kPyCharmppMNR:
 815     case kPyCharmppMNRwmi:
 816       break;
 817     case kPyD0PbPbMNR:
 818     case kPyD0pPbMNR:
 819     case kPyD0ppMNR:
 820       fParentSelect[0] =   421;
 821       fFlavorSelect    =   4;
 822       break;
 823     case kPyDPlusPbPbMNR:
 824     case kPyDPluspPbMNR:
 825     case kPyDPlusppMNR:
 826       fParentSelect[0] =   411;
 827       fFlavorSelect    =   4;
 828       break;
 829     case kPyDPlusStrangePbPbMNR:
 830     case kPyDPlusStrangepPbMNR:
 831     case kPyDPlusStrangeppMNR:
 832       fParentSelect[0] =   431;
 833       fFlavorSelect    =   4;
 834       break;
 835     case kPyLambdacppMNR:
 836       fParentSelect[0] =  4122;
 837       fFlavorSelect    =   4;
 838       break;
 839     case kPyBeauty:
 840     case kPyBeautyJets:
 841     case kPyBeautyPbPbMNR:
 842     case kPyBeautypPbMNR:
 843     case kPyBeautyppMNR:
 844     case kPyBeautyppMNRwmi:
 845        break;
 846     case kPyBeautyUnforced:
 847       fParentSelect[0] =  511;
 848       fParentSelect[1] =  521;
 849       fParentSelect[2] =  531;
 850       fParentSelect[3] = 5122;
 851       fParentSelect[4] = 5132;
 852       fParentSelect[5] = 5232;
 853       fParentSelect[6] = 5332;
 854       fFlavorSelect    = 5;
 855       break;
 856     case kPyJpsiChi:
 857     case kPyJpsi:
 858       fParentSelect[0] = 443;
 859       break;
 860     case kPyMbDefault:
 861     case kPyMbAtlasTuneMC09:
 862     case kPyMb:
 863     case kPyMbWithDirectPhoton:
 864     case kPyMbNonDiffr:
 865     case kPyMbMSEL1:
 866     case kPyJets:
 867     case kPyDirectGamma:
 868     case kPyLhwgMb:
 869       break;
 870     case kPyW:
 871     case kPyZ:
 872     case kPyMBRSingleDiffraction:
 873     case kPyMBRDoubleDiffraction:
 874     case kPyMBRCentralDiffraction:
 875       break;
 876   }
 877   //
 878   //
 879   //  JetFinder for Trigger
 880   //
 881   //  Configure detector (EMCAL like)
 882   //
 883   fPythia->SetPARU(51, fPycellEtaMax);
 884   fPythia->SetMSTU(51, fPycellNEta);
 885   fPythia->SetMSTU(52, fPycellNPhi);
 886   //
 887   //  Configure Jet Finder
 888   //
 889   fPythia->SetPARU(58,  fPycellThreshold);
 890   fPythia->SetPARU(52,  fPycellEtSeed);
 891   fPythia->SetPARU(53,  fPycellMinEtJet);
 892   fPythia->SetPARU(54,  fPycellMaxRadius);
 893   fPythia->SetMSTU(54,  2);
 894   //
 895   //  This counts the total number of calls to Pyevnt() per run.
 896   fTrialsRun = 0;
 897   fQ         = 0.;
 898   fX1        = 0.;
 899   fX2        = 0.;
 900   fNev       = 0 ;
 901   //
 902   //
 903   //
 904   AliGenMC::Init();
 905
 906   // Reset Lorentz boost if demanded
 907   if(!fUseLorentzBoost) {
 908     fDyBoost = 0;
 909     Warning("Init","Demand to discard Lorentz boost.\n");
 910   }
 911   //
 912   //
 913   //
 914   if (fSetNuclei) {
 915     fDyBoost = 0;
 916     Warning("Init","Deprecated function SetNuclei() used (nPDFs + no boost). Use SetProjectile + SetTarget + SetUseNuclearPDF + SetUseLorentzBoost instead.\n");
 917   }
 918   fPythia->SetPARJ(200, 0.0);
 919   fPythia->SetPARJ(199, 0.0);
 920   fPythia->SetPARJ(198, 0.0);
 921   fPythia->SetPARJ(197, 0.0);
 922
 923   if (fQuench == 1) {
 924     fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
 925   }
 926
 927   if(fQuench ==2){fPythia->SetPyquenParameters(fpyquenT,fpyquenTau,fpyquenNf,fpyquenEloss,fpyquenAngle);}
 928
 929   if (fQuench == 3) {
 930     // Nestor's change of the splittings
 931     fPythia->SetPARJ(200, 0.8);
 932     fPythia->SetMSTJ(41, 1);  // QCD radiation only
 933     fPythia->SetMSTJ(42, 2);  // angular ordering
 934     fPythia->SetMSTJ(44, 2);  // option to run alpha_s
 935     fPythia->SetMSTJ(47, 0);  // No correction back to hard scattering element
 936     fPythia->SetMSTJ(50, 0);  // No coherence in first branching
 937     fPythia->SetPARJ(82, 1.); // Cut off for parton showers
 938   } else if (fQuench == 4) {
 939     // Armesto-Cunqueiro-Salgado change of the splittings.
 940     AliFastGlauber* glauber = AliFastGlauber::Instance();
 941     glauber->Init(2);
 942     //read and store transverse almonds corresponding to differnt
 943     //impact parameters.
 944     glauber->SetCentralityClass(0.,0.1);
 945     fPythia->SetPARJ(200, 1.);
 946     fPythia->SetPARJ(198, fQhat);
 947     fPythia->SetPARJ(199, fLength);
 948     fPythia->SetMSTJ(42, 2);  // angular ordering
 949     fPythia->SetMSTJ(44, 2);  // option to run alpha_s
 950     fPythia->SetPARJ(82, 1.); // Cut off for parton showers
 951   }
 952
 953   if ( AliLog::GetDebugLevel("","AliGenPythia6Diff") >= 1 ) {
 954     fPythia->Pystat(4);
 955     fPythia->Pystat(5);
 956   }
 957 }
 958
 959 void AliGenPythia6Diff::Generate()
 960 {
 961   // Generate one event
 962   if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
 963   fDecayer->ForceDecay();
 964
 965   Double_t polar[3]   =   {0,0,0};
 966   Double_t origin[3]  =   {0,0,0};
 967   Double_t p[4];
 968   //  converts from mm/c to s
 969   const Float_t kconv=0.001/2.999792458e8;
 970   //
 971   Int_t nt=0;
 972   Int_t jev=0;
 973   Int_t j, kf;
 974   fTrials=0;
 975   fEventTime = 0.;
 976
 977
 978
 979   //  Set collision vertex position
 980   if (fVertexSmear == kPerEvent) Vertex();
 981
 982   //  event loop
 983   while(1)
 984   {
 985     //
 986     // Produce event
 987     //
 988     //
 989     // Switch hadronisation off
 990     //
 991     fPythia->SetMSTJ(1, 0);
 992
 993     if (fQuench ==4){
 994             Double_t bimp;
 995             // Quenching comes through medium-modified splitting functions.
 996             AliFastGlauber::Instance()->GetRandomBHard(bimp);
 997             fPythia->SetPARJ(197, bimp);
 998             fImpact = bimp;
 999             fPythia->Qpygin0();
1000     }
1001     //
1002     // Either produce new event or read partons from file
1003     //
1004     if (!fReadFromFile) {
1005             if (!fNewMIS) {
1006         fPythia->Pyevnt();
1007             } else {
1008         fPythia->Pyevnw();
1009             }
1010             fNpartons = fPythia->GetN();
1011     } else {
1012             printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
1013             fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
1014             fPythia->SetN(0);
1015             LoadEvent(fRL->Stack(), 0 , 1);
1016             fPythia->Pyedit(21);
1017     }
1018
1019     //
1020     //  Run quenching routine
1021     //
1022     if (fQuench == 1) {
1023             fPythia->Quench();
1024     } else if (fQuench == 2){
1025             fPythia->Pyquen(208., 0, 0.);
1026     } else if (fQuench == 3) {
1027             // Quenching is via multiplicative correction of the splittings
1028     }
1029
1030     //
1031     // Switch hadronisation on
1032     //
1033     if (fHadronisation) {
1034             fPythia->SetMSTJ(1, 1);
1035       //
1036       // .. and perform hadronisation
1037       //        printf("Calling hadronisation %d\n", fPythia->GetN());
1038
1039             if (fPatchOmegaDalitz) {
1040               fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
1041               fPythia->Pyexec();
1042               fPythia->DalitzDecays();
1043               fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
1044             }
1045             fPythia->Pyexec();
1046     }
1047     fTrials++;
1048     fPythia->ImportParticles(&fParticles,"All");
1049
1050     if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
1051     if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
1052
1053     //
1054     //
1055     //
1056     Int_t i;
1057
1058     fNprimaries = 0;
1059     Int_t np = fParticles.GetEntriesFast();
1060
1061     if (np == 0) continue;
1062     //
1063
1064     //
1065     Int_t* pParent   = new Int_t[np];
1066     Int_t* pSelected = new Int_t[np];
1067     Int_t* trackIt   = new Int_t[np];
1068     for (i = 0; i < np; i++) {
1069             pParent[i]   = -1;
1070             pSelected[i] =  0;
1071             trackIt[i]   =  0;
1072     }
1073
1074     Int_t nc = 0;        // Total n. of selected particles
1075     Int_t nParents = 0;  // Selected parents
1076     Int_t nTkbles = 0;   // Trackable particles
1077     if (fProcess != kPyMbDefault &&
1078         fProcess != kPyMb &&
1079         fProcess != kPyMbAtlasTuneMC09 &&
1080         fProcess != kPyMbWithDirectPhoton &&
1081         fProcess != kPyJets &&
1082         fProcess != kPyDirectGamma &&
1083         fProcess != kPyMbNonDiffr  &&
1084         fProcess != kPyMbMSEL1     &&
1085         fProcess != kPyW &&
1086         fProcess != kPyZ &&
1087         fProcess != kPyCharmppMNRwmi &&
1088         fProcess != kPyBeautyppMNRwmi &&
1089         fProcess != kPyBeautyJets) {
1090
1091             for (i = 0; i < np; i++) {
1092         TParticle* iparticle = (TParticle *) fParticles.At(i);
1093         Int_t ks = iparticle->GetStatusCode();
1094         kf = CheckPDGCode(iparticle->GetPdgCode());
1095         // No initial state partons
1096         if (ks==21) continue;
1097         //
1098         // Heavy Flavor Selection
1099         //
1100         // quark ?
1101         kf = TMath::Abs(kf);
1102         Int_t kfl = kf;
1103         // Resonance
1104
1105         if (kfl > 100000) kfl %= 100000;
1106         if (kfl > 10000)  kfl %= 10000;
1107         // meson ?
1108         if  (kfl > 10) kfl/=100;
1109         // baryon
1110         if (kfl > 10) kfl/=10;
1111         Int_t ipa = iparticle->GetFirstMother()-1;
1112         Int_t kfMo = 0;
1113         //
1114         // Establish mother daughter relation between heavy quarks and mesons
1115         //
1116         if (kf >= fFlavorSelect && kf <= 6) {
1117           Int_t idau = iparticle->GetFirstDaughter() - 1;
1118           if (idau > -1) {
1119             TParticle* daughter = (TParticle *) fParticles.At(idau);
1120             Int_t pdgD = daughter->GetPdgCode();
1121             if (pdgD == 91 || pdgD == 92) {
1122               Int_t jmin = daughter->GetFirstDaughter() - 1;
1123               Int_t jmax = daughter->GetLastDaughter()  - 1;
1124               for (Int_t jp = jmin; jp <= jmax; jp++)
1125                 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
1126             } // is string or cluster
1127           } // has daughter
1128         } // heavy quark
1129
1130
1131         if (ipa > -1) {
1132           TParticle *  mother = (TParticle *) fParticles.At(ipa);
1133           kfMo = TMath::Abs(mother->GetPdgCode());
1134         }
1135
1136         // What to keep in Stack?
1137         Bool_t flavorOK = kFALSE;
1138         Bool_t selectOK = kFALSE;
1139         if (fFeedDownOpt) {
1140           if (kfl >= fFlavorSelect) flavorOK = kTRUE;
1141         } else {
1142           if (kfl > fFlavorSelect) {
1143             nc = -1;
1144             break;
1145           }
1146           if (kfl == fFlavorSelect) flavorOK = kTRUE;
1147         }
1148         switch (fStackFillOpt) {
1149           case kHeavyFlavor:
1150           case kFlavorSelection:
1151             selectOK = kTRUE;
1152             break;
1153           case kParentSelection:
1154             if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
1155             break;
1156         }
1157         if (flavorOK && selectOK) {
1158           //
1159           // Heavy flavor hadron or quark
1160           //
1161           // Kinematic seletion on final state heavy flavor mesons
1162           if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
1163           {
1164             continue;
1165           }
1166           pSelected[i] = 1;
1167           if (ParentSelected(kf)) ++nParents; // Update parent count
1168           //                printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
1169         } else {
1170           // Kinematic seletion on decay products
1171           if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
1172               && !KinematicSelection(iparticle, 1))
1173           {
1174             continue;
1175           }
1176           //
1177           // Decay products
1178           // Select if mother was selected and is not tracked
1179
1180           if (pSelected[ipa] &&
1181               !trackIt[ipa]  &&     // mother will be  tracked ?
1182               kfMo !=  5 &&         // mother is b-quark, don't store fragments
1183               kfMo !=  4 &&         // mother is c-quark, don't store fragments
1184               kf   != 92)           // don't store string
1185           {
1186             //
1187             // Semi-stable or de-selected: diselect decay products:
1188             //
1189             //
1190             if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
1191             {
1192               Int_t ipF = iparticle->GetFirstDaughter();
1193               Int_t ipL = iparticle->GetLastDaughter();
1194               if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
1195             }
1196             //                  printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
1197             pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
1198           }
1199         }
1200         if (pSelected[i] == -1) pSelected[i] = 0;
1201         if (!pSelected[i]) continue;
1202         // Count quarks only if you did not include fragmentation
1203         if (fFragmentation && kf <= 10) continue;
1204
1205         nc++;
1206         // Decision on tracking
1207         trackIt[i] = 0;
1208         //
1209         // Track final state particle
1210         if (ks == 1) trackIt[i] = 1;
1211         // Track semi-stable particles
1212         if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
1213         // Track particles selected by process if undecayed.
1214         if (fForceDecay == kNoDecay) {
1215           if (ParentSelected(kf)) trackIt[i] = 1;
1216         } else {
1217           if (ParentSelected(kf)) trackIt[i] = 0;
1218         }
1219         if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
1220         //
1221         //
1222
1223       } // particle selection loop
1224             if (nc > 0) {
1225         for (i = 0; i<np; i++) {
1226           if (!pSelected[i]) continue;
1227           TParticle *  iparticle = (TParticle *) fParticles.At(i);
1228           kf = CheckPDGCode(iparticle->GetPdgCode());
1229           Int_t ks = iparticle->GetStatusCode();
1230           p[0] = iparticle->Px();
1231           p[1] = iparticle->Py();
1232           p[2] = iparticle->Pz();
1233           p[3] = iparticle->Energy();
1234
1235           origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1236           origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1237           origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1238
1239           Float_t tof   = fTime + kconv*iparticle->T();
1240           Int_t ipa     = iparticle->GetFirstMother()-1;
1241           Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
1242
1243           PushTrack(fTrackIt*trackIt[i], iparent, kf,
1244                     p[0], p[1], p[2], p[3],
1245                     origin[0], origin[1], origin[2], tof,
1246                     polar[0], polar[1], polar[2],
1247                     kPPrimary, nt, 1., ks);
1248           pParent[i] = nt;
1249           KeepTrack(nt);
1250           fNprimaries++;
1251         } //  PushTrack loop
1252             }
1253         } else {
1254             nc = GenerateMB();
1255     } // mb ?
1256
1257     GetSubEventTime();
1258
1259     delete[] pParent;
1260     delete[] pSelected;
1261     delete[] trackIt;
1262
1263     if (nc > 0) {
1264       switch (fCountMode) {
1265         case kCountAll:
1266           // printf(" Count all \n");
1267           jev += nc;
1268           break;
1269         case kCountParents:
1270           // printf(" Count parents \n");
1271           jev += nParents;
1272           break;
1273         case kCountTrackables:
1274           // printf(" Count trackable \n");
1275           jev += nTkbles;
1276           break;
1277       }
1278             if (jev >= fNpart || fNpart == -1) {
1279         fKineBias=Float_t(fNpart)/Float_t(fTrials);
1280
1281         fQ  += fPythia->GetVINT(51);
1282         fX1 += fPythia->GetVINT(41);
1283         fX2 += fPythia->GetVINT(42);
1284         fTrialsRun += fTrials;
1285         fNev++;
1286         MakeHeader();
1287         break;
1288             }
1289     }
1290   } // event loop
1291   SetHighWaterMark(nt);
1292   //  adjust weight due to kinematic selection
1293   //    AdjustWeights();
1294   //  get cross-section
1295   fXsection=fPythia->GetPARI(1);
1296 }
1297
1298 Int_t  AliGenPythia6Diff::GenerateMB()
1299 {
1300   //
1301   // Min Bias selection and other global selections
1302   //
1303   Int_t i, kf, nt, iparent;
1304   Int_t nc = 0;
1305   Double_t p[4];
1306   Double_t polar[3]   =   {0,0,0};
1307   Double_t origin[3]  =   {0,0,0};
1308   //  converts from mm/c to s
1309   const Float_t kconv=0.001/2.999792458e8;
1310
1311
1312   Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1313
1314
1315
1316   Int_t* pParent = new Int_t[np];
1317   for (i=0; i< np; i++) pParent[i] = -1;
1318   if ((fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi )
1319       && fEtMinJet > 0.) {
1320     TParticle* jet1 = (TParticle *) fParticles.At(6);
1321     TParticle* jet2 = (TParticle *) fParticles.At(7);
1322
1323     if (!jet1 || ! jet2 || !CheckTrigger(jet1, jet2)) {
1324       delete [] pParent;
1325       return 0;
1326     }
1327   }
1328
1329
1330   // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
1331   // implemented primaryly for kPyJets, but extended to any kind of process.
1332   if ((fFragPhotonInCalo || fPi0InCalo || fEtaInCalo || fEleInCalo || fHadronInCalo || fDecayPhotonInCalo) &&
1333       (fCheckPHOS || fCheckEMCAL || fCheckBarrel) ) {
1334     Bool_t ok = TriggerOnSelectedParticles(np);
1335
1336     if(!ok) {
1337       delete[] pParent;
1338       return 0;
1339     }
1340   }
1341
1342   // Check for diffraction
1343   if(fkTuneForDiff) {
1344     if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1) || ( TMath::Abs(fEnergyCMS - 8000) < 1 ) ) ) {
1345       if(!CheckDiffraction()) {
1346         delete [] pParent;
1347         return 0;
1348       }
1349     }
1350   }
1351
1352   // Check for minimum multiplicity
1353   if (fTriggerMultiplicity > 0) {
1354     Int_t multiplicity = 0;
1355     for (i = 0; i < np; i++) {
1356       TParticle *  iparticle = (TParticle *) fParticles.At(i);
1357
1358       Int_t statusCode = iparticle->GetStatusCode();
1359
1360       // Initial state particle
1361       if (statusCode != 1)
1362         continue;
1363       // eta cut
1364       if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
1365         continue;
1366       // pt cut
1367       if (iparticle->Pt() < fTriggerMultiplicityPtMin)
1368         continue;
1369
1370       TParticlePDG* pdgPart = iparticle->GetPDG();
1371       if (pdgPart && pdgPart->Charge() == 0)
1372         continue;
1373
1374       ++multiplicity;
1375     }
1376
1377     if (multiplicity < fTriggerMultiplicity) {
1378       delete [] pParent;
1379       return 0;
1380     }
1381     Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
1382   }
1383
1384
1385   if (fTriggerParticle) {
1386     Bool_t triggered = kFALSE;
1387     for (i = 0; i < np; i++) {
1388             TParticle *  iparticle = (TParticle *) fParticles.At(i);
1389             kf = CheckPDGCode(iparticle->GetPdgCode());
1390             if (kf != fTriggerParticle) continue;
1391             if (iparticle->Pt() == 0.) continue;
1392             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1393             if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
1394             triggered = kTRUE;
1395             break;
1396     }
1397     if (!triggered) {
1398       delete [] pParent;
1399       return 0;
1400     }
1401   }
1402
1403
1404   // Check if there is a ccbar or bbbar pair with at least one of the two
1405   // in fYMin < y < fYMax
1406
1407   if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
1408     TParticle *partCheck;
1409     TParticle *mother;
1410     Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
1411     Bool_t  theChild=kFALSE;
1412     Bool_t  triggered=kFALSE;
1413     Float_t y;
1414     Int_t   pdg,mpdg,mpdgUpperFamily;
1415     for(i=0; i<np; i++) {
1416       partCheck = (TParticle*)fParticles.At(i);
1417       pdg = partCheck->GetPdgCode();
1418       if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1419         if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1420         y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1421                            (partCheck->Energy()-partCheck->Pz()+1.e-13));
1422         if(y>fYMin && y<fYMax) inYcut=kTRUE;
1423       }
1424       if(fTriggerParticle) {
1425         if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1426       }
1427       if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1428         Int_t mi = partCheck->GetFirstMother() - 1;
1429         if(mi<0) continue;
1430         mother = (TParticle*)fParticles.At(mi);
1431         mpdg=TMath::Abs(mother->GetPdgCode());
1432         mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
1433         if ( ParentSelected(mpdg) ||
1434             (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1435           if (KinematicSelection(partCheck,1)) {
1436             theChild=kTRUE;
1437           }
1438         }
1439       }
1440     }
1441     if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
1442       delete[] pParent;
1443       return 0;
1444     }
1445     if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1446       delete[] pParent;
1447       return 0;
1448     }
1449     if(fTriggerParticle && !triggered) { // particle requested is not produced
1450       delete[] pParent;
1451       return 0;
1452     }
1453
1454   }
1455
1456   //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
1457   if ( (
1458
1459         fProcess == kPyW ||
1460         fProcess == kPyZ ||
1461         fProcess == kPyMbDefault ||
1462         fProcess == kPyMb ||
1463         fProcess == kPyMbAtlasTuneMC09 ||
1464         fProcess == kPyMbWithDirectPhoton ||
1465         fProcess == kPyMbNonDiffr)
1466       && (fCutOnChild == 1) ) {
1467     if ( !CheckKinematicsOnChild() ) {
1468       delete[] pParent;
1469       return 0;
1470     }
1471   }
1472
1473
1474   for (i = 0; i < np; i++) {
1475     Int_t trackIt = 0;
1476     TParticle *  iparticle = (TParticle *) fParticles.At(i);
1477     kf = CheckPDGCode(iparticle->GetPdgCode());
1478     Int_t ks = iparticle->GetStatusCode();
1479     Int_t km = iparticle->GetFirstMother();
1480     if (
1481         (((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
1482         ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)
1483         )
1484           {
1485       nc++;
1486             if (ks == 1) trackIt = 1;
1487             Int_t ipa = iparticle->GetFirstMother()-1;
1488
1489             iparent = (ipa > -1) ? pParent[ipa] : -1;
1490
1491       //
1492       // store track information
1493             p[0] = iparticle->Px();
1494             p[1] = iparticle->Py();
1495             p[2] = iparticle->Pz();
1496             p[3] = iparticle->Energy();
1497
1498
1499             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1500             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1501             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1502
1503             Float_t tof = fTime + fEventTime + kconv * iparticle->T();
1504
1505             PushTrack(fTrackIt*trackIt, iparent, kf,
1506                 p[0], p[1], p[2], p[3],
1507                 origin[0], origin[1], origin[2], tof,
1508                 polar[0], polar[1], polar[2],
1509                 kPPrimary, nt, 1., ks);
1510             fNprimaries++;
1511             KeepTrack(nt);
1512             pParent[i] = nt;
1513             SetHighWaterMark(nt);
1514
1515     } // select particle
1516   } // particle loop
1517
1518   delete[] pParent;
1519
1520   return 1;
1521 }
1522
1523
1524 void AliGenPythia6Diff::FinishRun()
1525 {
1526   // Print x-section summary
1527   fPythia->Pystat(1);
1528
1529   if (fNev > 0.) {
1530     fQ  /= fNev;
1531     fX1 /= fNev;
1532     fX2 /= fNev;
1533   }
1534
1535   printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1536   printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1537 }
1538
1539 void AliGenPythia6Diff::AdjustWeights() const
1540 {
1541   // Adjust the weights after generation of all events
1542   //
1543   if (gAlice) {
1544     TParticle *part;
1545     Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1546     for (Int_t i=0; i<ntrack; i++) {
1547             part= gAlice->GetMCApp()->Particle(i);
1548             part->SetWeight(part->GetWeight()*fKineBias);
1549     }
1550   }
1551 }
1552
1553 void AliGenPythia6Diff::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1554 {
1555   // Treat protons as inside nuclei with mass numbers a1 and a2
1556
1557   fAProjectile   = a1;
1558   fATarget       = a2;
1559   fNucPdf        = pdfset;  // 0 EKS98 9 EPS09LO 19 EPS09NLO
1560   fUseNuclearPDF = kTRUE;
1561   fSetNuclei     = kTRUE;
1562 }
1563
1564
1565 void AliGenPythia6Diff::MakeHeader()
1566 {
1567   //
1568   // Make header for the simulated event
1569   //
1570   if (gAlice) {
1571     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1572         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1573   }
1574
1575   // Builds the event header, to be called after each event
1576   if (fHeader) delete fHeader;
1577   fHeader = new AliGenPythiaEventHeader("Pythia");
1578   //
1579   // Event type
1580   if(fProcDiff>0){
1581     //      if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1582     //      printf("fPythia->GetMSTI(1) = %d   fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1583     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1584   }
1585   else
1586     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1587   //
1588   // Number of trials
1589   ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1590   //
1591   // Event Vertex
1592   fHeader->SetPrimaryVertex(fVertex);
1593   fHeader->SetInteractionTime(fTime+fEventTime);
1594   //
1595   // Number of primaries
1596   fHeader->SetNProduced(fNprimaries);
1597   //
1598   // Event weight
1599   fHeader->SetEventWeight(fPythia->GetVINT(97));
1600   //
1601   // Jets that have triggered
1602
1603   //Need to store jets for b-jet studies too!
1604   if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
1605   {
1606     Int_t ntrig, njet;
1607     Float_t jets[4][10];
1608     GetJets(njet, ntrig, jets);
1609
1610
1611     for (Int_t i = 0; i < ntrig; i++) {
1612             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1613                                                    jets[3][i]);
1614     }
1615   }
1616   //
1617   // Copy relevant information from external header, if present.
1618   //
1619   Float_t uqJet[4];
1620
1621   if (fRL) {
1622     AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1623     for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1624     {
1625             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1626
1627
1628             exHeader->TriggerJet(i, uqJet);
1629             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1630     }
1631   }
1632   //
1633   // Store quenching parameters
1634   //
1635   if (fQuench){
1636     Double_t z[4] = {0.};
1637     Double_t xp = 0.;
1638     Double_t yp = 0.;
1639
1640     if (fQuench == 1) {
1641             // Pythia::Quench()
1642             fPythia->GetQuenchingParameters(xp, yp, z);
1643     } else if (fQuench == 2){
1644             // Pyquen
1645             Double_t r1 = PARIMP.rb1;
1646             Double_t r2 = PARIMP.rb2;
1647             Double_t b  = PARIMP.b1;
1648             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1649             Double_t phi = PARIMP.psib1;
1650             xp = r * TMath::Cos(phi);
1651             yp = r * TMath::Sin(phi);
1652
1653     } else if (fQuench == 4) {
1654             // QPythia
1655             Double_t xy[2];
1656             Double_t i0i1[2];
1657             AliFastGlauber::Instance()->GetSavedXY(xy);
1658             AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1659             xp = xy[0];
1660             yp = xy[1];
1661             ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
1662     }
1663
1664     ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1665     ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1666   }
1667   //
1668   // Store pt^hard and cross-section
1669   ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1670 //  ((AliGenPythiaEventHeader*) fHeader)->SetXsection(fPythia->GetPARI(1));
1671
1672   //
1673   // Store Event Weight
1674   ((AliGenPythiaEventHeader*) fHeader)->SetEventWeight(fPythia->GetPARI(7));
1675
1676   //
1677   //  Pass header
1678   //
1679   AddHeader(fHeader);
1680   fHeader = 0x0;
1681 }
1682
1683 Bool_t AliGenPythia6Diff::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
1684 {
1685   // Check the kinematic trigger condition
1686   //
1687   Double_t eta[2];
1688   eta[0] = jet1->Eta();
1689   eta[1] = jet2->Eta();
1690   Double_t phi[2];
1691   phi[0] = jet1->Phi();
1692   phi[1] = jet2->Phi();
1693   Int_t    pdg[2];
1694   pdg[0] = jet1->GetPdgCode();
1695   pdg[1] = jet2->GetPdgCode();
1696   Bool_t   triggered = kFALSE;
1697
1698   if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess ==  kPyBeautyppMNRwmi) {
1699     Int_t njets = 0;
1700     Int_t ntrig = 0;
1701     Float_t jets[4][10];
1702     //
1703     // Use Pythia clustering on parton level to determine jet axis
1704     //
1705     GetJets(njets, ntrig, jets);
1706
1707     if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1708     //
1709   } else {
1710     Int_t ij = 0;
1711     Int_t ig = 1;
1712     if (pdg[0] == kGamma) {
1713             ij = 1;
1714             ig = 0;
1715     }
1716     //Check eta range first...
1717     if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1718         (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1719     {
1720             //Eta is okay, now check phi range
1721             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1722           (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1723             {
1724         triggered = kTRUE;
1725             }
1726     }
1727   }
1728   return triggered;
1729 }
1730
1731
1732
1733 Bool_t AliGenPythia6Diff::CheckKinematicsOnChild(){
1734   //
1735   //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1736   //
1737   Bool_t checking = kFALSE;
1738   Int_t j, kcode, ks, km;
1739   Int_t nPartAcc = 0; //number of particles in the acceptance range
1740   Int_t numberOfAcceptedParticles = 1;
1741   if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1742   Int_t npart = fParticles.GetEntriesFast();
1743
1744   for (j = 0; j<npart; j++) {
1745     TParticle *  jparticle = (TParticle *) fParticles.At(j);
1746     kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1747     ks = jparticle->GetStatusCode();
1748     km = jparticle->GetFirstMother();
1749
1750     if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1751             nPartAcc++;
1752     }
1753     if( numberOfAcceptedParticles <= nPartAcc){
1754       checking = kTRUE;
1755       break;
1756     }
1757   }
1758
1759   return checking;
1760 }
1761
1762 void  AliGenPythia6Diff::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1763 {
1764   //
1765   // Load event into Pythia Common Block
1766   //
1767
1768   Int_t npart = stack -> GetNprimary();
1769   Int_t n0 = 0;
1770
1771   if (!flag) {
1772     (fPythia->GetPyjets())->N = npart;
1773   } else {
1774     n0 = (fPythia->GetPyjets())->N;
1775     (fPythia->GetPyjets())->N = n0 + npart;
1776   }
1777
1778
1779   for (Int_t part = 0; part < npart; part++) {
1780     TParticle *mPart = stack->Particle(part);
1781
1782     Int_t kf     =  mPart->GetPdgCode();
1783     Int_t ks     =  mPart->GetStatusCode();
1784     Int_t idf    =  mPart->GetFirstDaughter();
1785     Int_t idl    =  mPart->GetLastDaughter();
1786
1787     if (reHadr) {
1788             if (ks == 11 || ks == 12) {
1789         ks  -= 10;
1790         idf  = -1;
1791         idl  = -1;
1792             }
1793     }
1794
1795     Float_t px = mPart->Px();
1796     Float_t py = mPart->Py();
1797     Float_t pz = mPart->Pz();
1798     Float_t e  = mPart->Energy();
1799     Float_t m  = mPart->GetCalcMass();
1800
1801
1802     (fPythia->GetPyjets())->P[0][part+n0] = px;
1803     (fPythia->GetPyjets())->P[1][part+n0] = py;
1804     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1805     (fPythia->GetPyjets())->P[3][part+n0] = e;
1806     (fPythia->GetPyjets())->P[4][part+n0] = m;
1807
1808     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1809     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1810     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1811     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1812     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1813   }
1814 }
1815
1816 void  AliGenPythia6Diff::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1817 {
1818   //
1819   // Load event into Pythia Common Block
1820   //
1821
1822   Int_t npart = stack -> GetEntries();
1823   Int_t n0 = 0;
1824
1825   if (!flag) {
1826     (fPythia->GetPyjets())->N = npart;
1827   } else {
1828     n0 = (fPythia->GetPyjets())->N;
1829     (fPythia->GetPyjets())->N = n0 + npart;
1830   }
1831
1832
1833   for (Int_t part = 0; part < npart; part++) {
1834     TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1835     if (!mPart) continue;
1836
1837     Int_t kf     =  mPart->GetPdgCode();
1838     Int_t ks     =  mPart->GetStatusCode();
1839     Int_t idf    =  mPart->GetFirstDaughter();
1840     Int_t idl    =  mPart->GetLastDaughter();
1841
1842     if (reHadr) {
1843       if (ks == 11 || ks == 12) {
1844         ks  -= 10;
1845         idf  = -1;
1846         idl  = -1;
1847       }
1848     }
1849
1850     Float_t px = mPart->Px();
1851     Float_t py = mPart->Py();
1852     Float_t pz = mPart->Pz();
1853     Float_t e  = mPart->Energy();
1854     Float_t m  = mPart->GetCalcMass();
1855
1856
1857     (fPythia->GetPyjets())->P[0][part+n0] = px;
1858     (fPythia->GetPyjets())->P[1][part+n0] = py;
1859     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1860     (fPythia->GetPyjets())->P[3][part+n0] = e;
1861     (fPythia->GetPyjets())->P[4][part+n0] = m;
1862
1863     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1864     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1865     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1866     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1867     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1868   }
1869 }
1870
1871
1872 //void AliGenPythia6Diff::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1873 //{
1874 //  //
1875 //  //  Calls the Pythia jet finding algorithm to find jets in the current event
1876 //  //
1877 //  //
1878 //  //
1879 //  //  Save jets
1880 //  Int_t n     = fPythia->GetN();
1881 //
1882 //  //
1883 //  //  Run Jet Finder
1884 //  fPythia->Pycell(njets);
1885 //  Int_t i;
1886 //  for (i = 0; i < njets; i++) {
1887 //    Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1888 //    Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1889 //    Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1890 //    Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1891 //
1892 //    jets[0][i] = px;
1893 //    jets[1][i] = py;
1894 //    jets[2][i] = pz;
1895 //    jets[3][i] = e;
1896 //  }
1897 //}
1898 //
1899 //
1900
1901 void  AliGenPythia6Diff::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1902 {
1903   //
1904   //  Calls the Pythia clustering algorithm to find jets in the current event
1905   //
1906   Int_t n     = fPythia->GetN();
1907   nJets       = 0;
1908   nJetsTrig   = 0;
1909   if (fJetReconstruction == kCluster) {
1910     //
1911     //  Configure cluster algorithm
1912     //
1913     fPythia->SetPARU(43, 2.);
1914     fPythia->SetMSTU(41, 1);
1915     //
1916     //  Call cluster algorithm
1917     //
1918     fPythia->Pyclus(nJets);
1919     //
1920     //  Loading jets from common block
1921     //
1922   } else {
1923
1924     //
1925     //  Run Jet Finder
1926     fPythia->Pycell(nJets);
1927   }
1928
1929   Int_t i;
1930   for (i = 0; i < nJets; i++) {
1931     Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1932     Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1933     Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1934     Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1935     Float_t pt    = TMath::Sqrt(px * px + py * py);
1936     Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1937     Float_t theta = TMath::ATan2(pt,pz);
1938     Float_t et    = e * TMath::Sin(theta);
1939     Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1940     if (
1941         eta > fEtaMinJet && eta < fEtaMaxJet &&
1942         phi > fPhiMinJet && phi < fPhiMaxJet &&
1943         et  > fEtMinJet  && et  < fEtMaxJet
1944         )
1945     {
1946             jets[0][nJetsTrig] = px;
1947             jets[1][nJetsTrig] = py;
1948             jets[2][nJetsTrig] = pz;
1949             jets[3][nJetsTrig] = e;
1950             nJetsTrig++;
1951       //            printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1952     } else {
1953       //            printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1954     }
1955   }
1956 }
1957
1958 void AliGenPythia6Diff::GetSubEventTime()
1959 {
1960   // Calculates time of the next subevent
1961   fEventTime = 0.;
1962   if (fEventsTime) {
1963     TArrayF &array = *fEventsTime;
1964     fEventTime = array[fCurSubEvent++];
1965   }
1966   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1967   return;
1968 }
1969
1970 Bool_t AliGenPythia6Diff::IsInBarrel(Float_t eta) const
1971 {
1972   // Is particle in Central Barrel acceptance?
1973   // etamin=-etamax
1974   if( eta < fTriggerEta  )
1975     return kTRUE;
1976   else
1977     return kFALSE;
1978 }
1979
1980 Bool_t AliGenPythia6Diff::IsInEMCAL(Float_t phi, Float_t eta) const
1981 {
1982   // Is particle in EMCAL acceptance?
1983   // phi in degrees, etamin=-etamax
1984   if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi &&
1985      eta < fEMCALEta  )
1986     return kTRUE;
1987   else
1988     return kFALSE;
1989 }
1990
1991 Bool_t AliGenPythia6Diff::IsInPHOS(Float_t phi, Float_t eta, Int_t iparticle)
1992 {
1993   // Is particle in PHOS acceptance?
1994   // Acceptance slightly larger considered.
1995   // phi in degrees, etamin=-etamax
1996   // iparticle is the index of the particle to be checked, for PHOS rotation case
1997
1998   if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi &&
1999      eta < fPHOSEta  )
2000     return kTRUE;
2001   else
2002   {
2003     if( fCheckPHOSeta && eta < fPHOSEta) fPHOSRotateCandidate = iparticle;
2004
2005     return kFALSE;
2006   }
2007 }
2008
2009 void AliGenPythia6Diff::RotatePhi(Bool_t& okdd)
2010 {
2011   //Rotate event in phi to enhance events in PHOS acceptance
2012
2013   if(fPHOSRotateCandidate < 0) return ;
2014
2015   //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
2016   Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
2017   Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
2018   Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
2019
2020   //calculate deltaphi
2021   TParticle* ph = (TParticle *) fParticles.At(fPHOSRotateCandidate);
2022   Double_t phphi = ph->Phi();
2023   Double_t deltaphi = phiPHOS - phphi;
2024
2025
2026
2027   //loop for all particles and produce the phi rotation
2028   Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
2029   Double_t oldphi, newphi;
2030   Double_t newVx, newVy, r, vZ, time;
2031   Double_t newPx, newPy, pt, pz, e;
2032   for(Int_t i=0; i< np; i++) {
2033     TParticle* iparticle = (TParticle *) fParticles.At(i);
2034     oldphi = iparticle->Phi();
2035     newphi = oldphi + deltaphi;
2036     if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
2037     if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
2038
2039     r = iparticle->R();
2040     newVx = r * TMath::Cos(newphi);
2041     newVy = r * TMath::Sin(newphi);
2042     vZ   = iparticle->Vz(); // don't transform
2043     time = iparticle->T(); // don't transform
2044
2045     pt = iparticle->Pt();
2046     newPx = pt * TMath::Cos(newphi);
2047     newPy = pt * TMath::Sin(newphi);
2048     pz = iparticle->Pz(); // don't transform
2049     e  = iparticle->Energy(); // don't transform
2050
2051     // apply rotation
2052     iparticle->SetProductionVertex(newVx, newVy, vZ, time);
2053     iparticle->SetMomentum(newPx, newPy, pz, e);
2054
2055   } //end particle loop
2056
2057   // now let's check that we put correctly the candidate photon in PHOS
2058   Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
2059   Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
2060   if(IsInPHOS(phi,eta,-1))
2061     okdd = kTRUE;
2062
2063   // reset the value for next event
2064   fPHOSRotateCandidate = -1;
2065
2066 }
2067
2068
2069 Bool_t AliGenPythia6Diff::CheckDiffraction()
2070 {
2071   // use this method only with Perugia-0 tune!
2072
2073   //  printf("AAA\n");
2074
2075   Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
2076
2077   Int_t iPart1=-1;
2078   Int_t iPart2=-1;
2079
2080   Double_t y1 = 1e10;
2081   Double_t y2 = -1e10;
2082
2083   const Int_t kNstable=20;
2084   const Int_t pdgStable[20] = {
2085     22,             // Photon
2086     11,             // Electron
2087     12,             // Electron Neutrino
2088     13,             // Muon
2089     14,             // Muon Neutrino
2090     15,             // Tau
2091     16,             // Tau Neutrino
2092     211,            // Pion
2093     321,            // Kaon
2094     311,            // K0
2095     130,            // K0s
2096     310,            // K0l
2097     2212,           // Proton
2098     2112,           // Neutron
2099     3122,           // Lambda_0
2100     3112,           // Sigma Minus
2101     3222,           // Sigma Plus
2102     3312,           // Xsi Minus
2103     3322,           // Xsi0
2104     3334            // Omega
2105   };
2106
2107   for (Int_t i = 0; i < np; i++) {
2108     TParticle *  part = (TParticle *) fParticles.At(i);
2109
2110     Int_t statusCode = part->GetStatusCode();
2111
2112     // Initial state particle
2113     if (statusCode != 1)
2114       continue;
2115
2116     Int_t pdg = TMath::Abs(part->GetPdgCode());
2117     Bool_t isStable = kFALSE;
2118     for (Int_t i1 = 0; i1 < kNstable; i1++) {
2119       if (pdg == pdgStable[i1]) {
2120         isStable = kTRUE;
2121         break;
2122       }
2123     }
2124     if(!isStable)
2125       continue;
2126
2127     Double_t y = part->Y();
2128
2129     if (y < y1)
2130           {
2131             y1 = y;
2132             iPart1 = i;
2133           }
2134     if (y > y2)
2135     {
2136       y2 = y;
2137       iPart2 = i;
2138     }
2139   }
2140
2141   if(iPart1<0 || iPart2<0) return kFALSE;
2142
2143   y1=TMath::Abs(y1);
2144   y2=TMath::Abs(y2);
2145
2146   TParticle *  part1 = (TParticle *) fParticles.At(iPart1);
2147   TParticle *  part2 = (TParticle *) fParticles.At(iPart2);
2148
2149   Int_t pdg1 = part1->GetPdgCode();
2150   Int_t pdg2 = part2->GetPdgCode();
2151
2152
2153   Int_t iPart = -1;
2154   if (pdg1 == 2212 && pdg2 == 2212)
2155   {
2156     if(y1 > y2)
2157       iPart = iPart1;
2158     else if(y1 < y2)
2159       iPart = iPart2;
2160     else {
2161       iPart = iPart1;
2162       if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
2163     }
2164   }
2165   else if (pdg1 == 2212)
2166     iPart = iPart1;
2167   else if (pdg2 == 2212)
2168     iPart = iPart2;
2169
2170
2171
2172
2173
2174   Double_t M=-1.;
2175   if(iPart>0) {
2176     TParticle *  part = (TParticle *) fParticles.At(iPart);
2177     Double_t E= part->Energy();
2178     Double_t P= part->P();
2179     Double_t M2 = (fEnergyCMS-E-P)*(fEnergyCMS-E+P);
2180     if(M2<0)  return kFALSE;
2181     M= TMath::Sqrt(M2);
2182   }
2183
2184   Double_t Mmin, Mmax, wSD, wDD, wND;
2185   if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
2186
2187   if(M>-1 && M<Mmin) return kFALSE;
2188   if(M>Mmax) M=-1;
2189
2190   Int_t procType=fPythia->GetMSTI(1);
2191   Int_t proc0=2;
2192   if(procType== 94) proc0=1;
2193   if(procType== 92 || procType== 93) proc0=0;
2194
2195   Int_t proc=2;
2196   if(M>0) proc=0;
2197   else if(proc0==1) proc=1;
2198
2199   if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
2200   if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
2201
2202
2203   //     if(proc==1 || proc==2) return kFALSE;
2204
2205   if(proc!=0) {
2206     if(proc0!=0) fProcDiff = procType;
2207     else       fProcDiff = 95;
2208     return kTRUE;
2209   }
2210
2211   if(wSD<0)  AliError("wSD<0 ! \n");
2212
2213   if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> wSD) return kFALSE;
2214
2215   //    printf("iPart = %d\n", iPart);
2216
2217   if(iPart==iPart1) fProcDiff=93;
2218   else if(iPart==iPart2) fProcDiff=92;
2219   else {
2220     printf("EROOR:  iPart!=iPart1 && iPart!=iPart2\n");
2221
2222   }
2223
2224   return kTRUE;
2225 }
2226
2227
2228
2229 Bool_t AliGenPythia6Diff::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax,
2230                                                 Double_t &wSD, Double_t &wDD, Double_t &wND)
2231 {
2232
2233   // 900 GeV
2234   if(TMath::Abs(fEnergyCMS-900)<1 ){
2235
2236     const Int_t nbin=400;
2237     Double_t bin[]={
2238       1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
2239       4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
2240       7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
2241       10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
2242       13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
2243       15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
2244       18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
2245       21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
2246       24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
2247       27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
2248       30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
2249       33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
2250       36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
2251       39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
2252       42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
2253       45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
2254       48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
2255       51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
2256       54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
2257       57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
2258       60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
2259       63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
2260       66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
2261       69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
2262       72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
2263       75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
2264       78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
2265       81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
2266       84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
2267       87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
2268       90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
2269       93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
2270       96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
2271       99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2272       102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2273       105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2274       108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2275       111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2276       114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2277       117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2278       120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2279       123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2280       126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2281       129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2282       132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2283       135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2284       138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2285       141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2286       144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2287       147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2288       150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2289       153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2290       156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2291       159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2292       162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2293       165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2294       168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2295       171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2296       174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2297       177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2298       180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2299       183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2300       186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2301       189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2302       192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2303       195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2304       198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2305     Double_t w[]={
2306       1.000000, 0.643457, 0.645609, 0.648347, 0.604563, 0.605002,
2307       0.602819, 0.611473, 0.576412, 0.562354, 0.550216, 0.529285,
2308       0.534558, 0.534364, 0.530358, 0.518475, 0.489253, 0.469754,
2309       0.469825, 0.450513, 0.455849, 0.435312, 0.437210, 0.456686,
2310       0.413577, 0.427093, 0.426894, 0.418834, 0.409475, 0.388483,
2311       0.412208, 0.388912, 0.389611, 0.382572, 0.389220, 0.370964,
2312       0.380463, 0.370873, 0.363701, 0.369363, 0.357361, 0.365759,
2313       0.348566, 0.337062, 0.348190, 0.332330, 0.359001, 0.335836,
2314       0.339154, 0.335599, 0.336035, 0.335204, 0.353440, 0.337836,
2315       0.333874, 0.307120, 0.294963, 0.324978, 0.313359, 0.317688,
2316       0.323758, 0.319304, 0.335317, 0.301765, 0.317257, 0.356331,
2317       0.323119, 0.297732, 0.303188, 0.326102, 0.316467, 0.294728,
2318       0.308135, 0.288857, 0.325692, 0.312493, 0.291100, 0.325921,
2319       0.313317, 0.295980, 0.308481, 0.328380, 0.313081, 0.296763,
2320       0.295431, 0.317325, 0.320462, 0.286918, 0.316035, 0.335208,
2321       0.283193, 0.333945, 0.292534, 0.294164, 0.330680, 0.296992,
2322       0.285509, 0.317260, 0.298269, 0.311299, 0.312129, 0.286822,
2323       0.287442, 0.319139, 0.283314, 0.318454, 0.297727, 0.301597,
2324       0.282483, 0.294792, 0.305569, 0.290957, 0.297817, 0.282908,
2325       0.272401, 0.305584, 0.300220, 0.297020, 0.298781, 0.278008,
2326       0.277727, 0.323777, 0.287419, 0.342074, 0.287259, 0.303658,
2327       0.302668, 0.279622, 0.280586, 0.313630, 0.276068, 0.257051,
2328       0.309996, 0.265534, 0.297138, 0.281738, 0.294610, 0.292882,
2329       0.286860, 0.312686, 0.293244, 0.293744, 0.271375, 0.278734,
2330       0.280308, 0.304739, 0.287907, 0.285261, 0.311180, 0.313476,
2331       0.289660, 0.289282, 0.319505, 0.271285, 0.272008, 0.289245,
2332       0.281038, 0.285284, 0.295836, 0.281416, 0.283501, 0.295417,
2333       0.304372, 0.297764, 0.291378, 0.321530, 0.315604, 0.329507,
2334       0.282609, 0.275576, 0.283721, 0.311714, 0.283760, 0.273188,
2335       0.312193, 0.264347, 0.281532, 0.301226, 0.281718, 0.336408,
2336       0.283157, 0.332010, 0.289974, 0.290256, 0.301569, 0.332228,
2337       0.288282, 0.326339, 0.313653, 0.300361, 0.289470, 0.264830,
2338       0.298659, 0.272359, 0.278878, 0.306001, 0.328168, 0.294991,
2339       0.327737, 0.278056, 0.302435, 0.284183, 0.279270, 0.307279,
2340       0.307917, 0.315196, 0.283803, 0.313333, 0.315730, 0.304818,
2341       0.307171, 0.295223, 0.333741, 0.346911, 0.310143, 0.336686,
2342       0.275459, 0.334781, 0.295405, 0.275816, 0.301678, 0.327242,
2343       0.320717, 0.309230, 0.292145, 0.294489, 0.305088, 0.300969,
2344       0.277438, 0.326159, 0.297065, 0.301177, 0.303843, 0.275382,
2345       0.304019, 0.284166, 0.289610, 0.331611, 0.317131, 0.310880,
2346       0.360456, 0.294052, 0.342694, 0.327166, 0.336797, 0.298040,
2347       0.295767, 0.260053, 0.325544, 0.335310, 0.320182, 0.301072,
2348       0.313117, 0.283407, 0.299206, 0.293525, 0.305067, 0.255978,
2349       0.327055, 0.316382, 0.317700, 0.278993, 0.283120, 0.314000,
2350       0.274396, 0.291208, 0.348813, 0.319603, 0.313076, 0.289155,
2351       0.343988, 0.311426, 0.322896, 0.328726, 0.337062, 0.389093,
2352       0.284122, 0.312184, 0.304008, 0.319170, 0.320778, 0.288390,
2353       0.337272, 0.356273, 0.343310, 0.312209, 0.330709, 0.297977,
2354       0.346146, 0.369162, 0.324385, 0.339831, 0.337037, 0.318739,
2355       0.343157, 0.277720, 0.368407, 0.321330, 0.338997, 0.314220,
2356       0.328861, 0.321824, 0.328013, 0.356925, 0.359144, 0.296314,
2357       0.345415, 0.396711, 0.347032, 0.294928, 0.343799, 0.322331,
2358       0.328656, 0.326098, 0.337338, 0.337038, 0.300179, 0.351391,
2359       0.324337, 0.330896, 0.302842, 0.310522, 0.337052, 0.359989,
2360       0.383250, 0.359355, 0.315382, 0.333113, 0.342598, 0.355348,
2361       0.320751, 0.320475, 0.351762, 0.351475, 0.338358, 0.326153,
2362       0.302507, 0.340048, 0.318685, 0.381646, 0.339320, 0.299453,
2363       0.426599, 0.393515, 0.353929, 0.328435, 0.413976, 0.292558,
2364       0.379340, 0.358344, 0.409259, 0.313821, 0.336675, 0.324521,
2365       0.408382, 0.346273, 0.312939, 0.362453, 0.343152, 0.330577,
2366       0.332831, 0.353299, 0.347745, 0.334818, 0.332234, 0.385585,
2367       0.395483, 0.395316, 0.326972, 0.349434, 0.317604, 0.328980,
2368       0.375056, 0.317290, 0.357083, 0.346165, 0.310444, 0.356873,
2369       0.359523, 0.365308, 0.365122, 0.383685, 0.370975, 0.396928,
2370       0.407654, 0.307755, 0.323033, 0.350580, 0.345231, 0.342462,
2371       0.400000, 0.318309, 0.403570, 0.322856, 0.383053, 0.422252,
2372       0.386112, 0.364314, 0.434375, 0.334629};
2373     wDD = 0.379611;
2374     wND = 0.496961;
2375     wSD = -1;
2376
2377     Mmin = bin[0];
2378     Mmax = bin[nbin];
2379     if(M<Mmin || M>Mmax) return kTRUE;
2380
2381     Int_t ibin=nbin-1;
2382     for(Int_t i=1; i<=nbin; i++)
2383       if(M<=bin[i]) {
2384         ibin=i-1;
2385         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2386         break;
2387       }
2388     wSD=w[ibin];
2389     return kTRUE;
2390   }
2391   else if(TMath::Abs(fEnergyCMS-2760)<1 ){
2392
2393     const Int_t nbin=400;
2394     Double_t bin[]={
2395       1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
2396       4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
2397       7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
2398       10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
2399       13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
2400       15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
2401       18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
2402       21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
2403       24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
2404       27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
2405       30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
2406       33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
2407       36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
2408       39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
2409       42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
2410       45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
2411       48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
2412       51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
2413       54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
2414       57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
2415       60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
2416       63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
2417       66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
2418       69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
2419       72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
2420       75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
2421       78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
2422       81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
2423       84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
2424       87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
2425       90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
2426       93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
2427       96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
2428       99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2429       102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2430       105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2431       108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2432       111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2433       114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2434       117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2435       120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2436       123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2437       126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2438       129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2439       132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2440       135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2441       138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2442       141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2443       144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2444       147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2445       150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2446       153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2447       156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2448       159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2449       162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2450       165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2451       168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2452       171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2453       174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2454       177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2455       180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2456       183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2457       186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2458       189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2459       192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2460       195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2461       198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2462     Double_t w[]={
2463       1.000000, 0.692593, 0.673384, 0.666273, 0.657285, 0.637723,
2464       0.625881, 0.643590, 0.606100, 0.589007, 0.567824, 0.578705,
2465       0.538530, 0.517937, 0.528278, 0.515918, 0.539461, 0.466186,
2466       0.489869, 0.468402, 0.465017, 0.453336, 0.460769, 0.474638,
2467       0.456347, 0.434471, 0.427478, 0.435435, 0.410934, 0.366431,
2468       0.382240, 0.379513, 0.394249, 0.386837, 0.353103, 0.382138,
2469       0.377497, 0.389479, 0.378736, 0.347933, 0.354605, 0.352077,
2470       0.324443, 0.358792, 0.339968, 0.359052, 0.330734, 0.318291,
2471       0.333703, 0.358644, 0.335819, 0.332213, 0.309051, 0.309975,
2472       0.331626, 0.304407, 0.309819, 0.312097, 0.312462, 0.320411,
2473       0.280401, 0.302311, 0.315863, 0.281479, 0.310003, 0.296911,
2474       0.313676, 0.281071, 0.294163, 0.306500, 0.283462, 0.274867,
2475       0.307149, 0.270555, 0.282264, 0.287373, 0.307849, 0.278675,
2476       0.286990, 0.278269, 0.300105, 0.286799, 0.265674, 0.275140,
2477       0.285702, 0.257352, 0.267714, 0.248204, 0.252220, 0.255678,
2478       0.282946, 0.268257, 0.282375, 0.262675, 0.275564, 0.248345,
2479       0.236259, 0.291914, 0.259936, 0.241338, 0.267389, 0.285044,
2480       0.289419, 0.253594, 0.284568, 0.231840, 0.260008, 0.268527,
2481       0.275363, 0.224115, 0.281260, 0.257913, 0.295152, 0.264399,
2482       0.232287, 0.282533, 0.223431, 0.255756, 0.244471, 0.221695,
2483       0.272450, 0.284244, 0.253682, 0.270717, 0.275403, 0.240323,
2484       0.245081, 0.241859, 0.216340, 0.244789, 0.220291, 0.238478,
2485       0.224691, 0.244058, 0.266117, 0.271478, 0.242012, 0.267321,
2486       0.248494, 0.253343, 0.255606, 0.235458, 0.241079, 0.233223,
2487       0.226813, 0.259224, 0.234239, 0.258606, 0.210892, 0.238186,
2488       0.243271, 0.222678, 0.213437, 0.273939, 0.247966, 0.232548,
2489       0.263438, 0.222089, 0.272111, 0.248818, 0.244295, 0.238368,
2490       0.236908, 0.248776, 0.232604, 0.231194, 0.227117, 0.231152,
2491       0.282140, 0.229778, 0.232631, 0.261794, 0.216633, 0.253471,
2492       0.242157, 0.227406, 0.269335, 0.230547, 0.210618, 0.285872,
2493       0.248776, 0.229875, 0.242728, 0.227388, 0.220567, 0.222062,
2494       0.235950, 0.224087, 0.228895, 0.208287, 0.235999, 0.208696,
2495       0.230367, 0.267667, 0.220484, 0.233402, 0.233815, 0.250455,
2496       0.253120, 0.219556, 0.230980, 0.236661, 0.222395, 0.226111,
2497       0.198315, 0.210555, 0.202813, 0.208594, 0.235976, 0.221490,
2498       0.243059, 0.204901, 0.216987, 0.229039, 0.231466, 0.221975,
2499       0.231220, 0.253638, 0.250448, 0.260291, 0.328345, 0.205739,
2500       0.222014, 0.251513, 0.279427, 0.270506, 0.248409, 0.222472,
2501       0.291632, 0.227796, 0.248769, 0.276896, 0.214742, 0.200139,
2502       0.230693, 0.226031, 0.268900, 0.185160, 0.245353, 0.205843,
2503       0.231155, 0.219122, 0.214811, 0.199763, 0.274179, 0.217598,
2504       0.274988, 0.237244, 0.211820, 0.225459, 0.252799, 0.235948,
2505       0.224986, 0.245385, 0.237770, 0.213373, 0.229737, 0.215487,
2506       0.234453, 0.249684, 0.239435, 0.250422, 0.257194, 0.214762,
2507       0.212266, 0.228988, 0.253798, 0.201607, 0.239946, 0.205245,
2508       0.231670, 0.212774, 0.206768, 0.231563, 0.189388, 0.227926,
2509       0.227203, 0.237754, 0.221628, 0.211138, 0.203322, 0.200985,
2510       0.231780, 0.220294, 0.232686, 0.234243, 0.218264, 0.255870,
2511       0.223213, 0.238670, 0.213713, 0.213064, 0.246700, 0.233446,
2512       0.221503, 0.206767, 0.200722, 0.226179, 0.237425, 0.239229,
2513       0.238611, 0.240419, 0.247806, 0.215923, 0.205298, 0.232778,
2514       0.272312, 0.226773, 0.258103, 0.223287, 0.269404, 0.203398,
2515       0.223782, 0.204213, 0.229664, 0.234040, 0.228419, 0.203936,
2516       0.263686, 0.199141, 0.236127, 0.214058, 0.204611, 0.224324,
2517       0.292140, 0.190735, 0.235157, 0.213018, 0.257085, 0.190554,
2518       0.203197, 0.213044, 0.237023, 0.214243, 0.193562, 0.262403,
2519       0.206256, 0.221396, 0.233588, 0.256611, 0.249731, 0.226683,
2520       0.199330, 0.251026, 0.222596, 0.201941, 0.186374, 0.221038,
2521       0.196555, 0.222560, 0.299419, 0.231979, 0.242924, 0.198310,
2522       0.217628, 0.235458, 0.278595, 0.218624, 0.277305, 0.239109,
2523       0.205600, 0.253715, 0.221173, 0.218195, 0.277647, 0.241974,
2524       0.268748, 0.268128, 0.292258, 0.249420, 0.191034, 0.219506,
2525       0.216502, 0.250677, 0.193386, 0.201310, 0.259464, 0.255351,
2526       0.269628, 0.221063, 0.294079, 0.196726, 0.233634, 0.221870,
2527       0.216236, 0.197259, 0.247433, 0.272765, 0.294079, 0.236336,
2528       0.206396, 0.238524, 0.247846, 0.269519, 0.237141, 0.230611,
2529       0.201712, 0.242225, 0.255565, 0.258738};
2530     wDD = 0.512813;
2531     wND = 0.518820;
2532     wSD = -1;
2533
2534     Mmin = bin[0];
2535     Mmax = bin[nbin];
2536     if(M<Mmin || M>Mmax) return kTRUE;
2537
2538     Int_t ibin=nbin-1;
2539     for(Int_t i=1; i<=nbin; i++)
2540       if(M<=bin[i]) {
2541         ibin=i-1;
2542         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2543         break;
2544       }
2545     wSD=w[ibin];
2546     return kTRUE;
2547   }
2548
2549
2550   else if(TMath::Abs(fEnergyCMS-7000)<1  || TMath::Abs(fEnergyCMS-8000)<1){
2551     const Int_t nbin=400;
2552     Double_t bin[]={
2553       1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
2554       4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
2555       7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
2556       10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
2557       13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
2558       15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
2559       18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
2560       21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
2561       24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
2562       27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
2563       30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
2564       33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
2565       36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
2566       39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
2567       42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
2568       45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
2569       48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
2570       51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
2571       54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
2572       57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
2573       60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
2574       63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
2575       66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
2576       69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
2577       72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
2578       75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
2579       78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
2580       81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
2581       84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
2582       87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
2583       90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
2584       93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
2585       96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
2586       99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2587       102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2588       105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2589       108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2590       111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2591       114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2592       117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2593       120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2594       123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2595       126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2596       129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2597       132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2598       135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2599       138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2600       141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2601       144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2602       147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2603       150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2604       153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2605       156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2606       159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2607       162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2608       165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2609       168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2610       171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2611       174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2612       177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2613       180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2614       183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2615       186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2616       189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2617       192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2618       195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2619       198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2620     Double_t w[]={
2621       1.000000, 0.798640, 0.770197, 0.717990, 0.699434, 0.692093,
2622       0.620940, 0.655294, 0.636496, 0.633483, 0.600518, 0.594548,
2623       0.550759, 0.550413, 0.528926, 0.525877, 0.506701, 0.504117,
2624       0.486973, 0.492930, 0.461833, 0.488695, 0.438874, 0.460274,
2625       0.451551, 0.454424, 0.449270, 0.387571, 0.427556, 0.428740,
2626       0.401682, 0.402260, 0.408961, 0.395417, 0.387426, 0.378602,
2627       0.357826, 0.359125, 0.348494, 0.363710, 0.356117, 0.340363,
2628       0.337637, 0.396084, 0.341777, 0.340551, 0.348838, 0.344014,
2629       0.340468, 0.317654, 0.355584, 0.326023, 0.373416, 0.312298,
2630       0.326522, 0.290540, 0.335557, 0.318689, 0.327544, 0.319501,
2631       0.331754, 0.312728, 0.282263, 0.274937, 0.303867, 0.307820,
2632       0.289344, 0.268934, 0.288908, 0.290018, 0.291369, 0.295242,
2633       0.289067, 0.277685, 0.267957, 0.267559, 0.320229, 0.265060,
2634       0.305931, 0.305352, 0.262064, 0.281879, 0.287780, 0.270033,
2635       0.270814, 0.276667, 0.271531, 0.275283, 0.258189, 0.287969,
2636       0.251247, 0.301527, 0.267230, 0.245860, 0.293125, 0.253421,
2637       0.272396, 0.243637, 0.236206, 0.278452, 0.246544, 0.263165,
2638       0.230484, 0.231102, 0.258281, 0.244707, 0.270111, 0.248295,
2639       0.246942, 0.245592, 0.272766, 0.254243, 0.199647, 0.262590,
2640       0.226710, 0.243836, 0.214153, 0.233535, 0.235207, 0.234245,
2641       0.247698, 0.248379, 0.241463, 0.265230, 0.223242, 0.236191,
2642       0.252700, 0.231865, 0.228352, 0.229101, 0.237385, 0.246485,
2643       0.254706, 0.245565, 0.224841, 0.257301, 0.240968, 0.202379,
2644       0.236943, 0.241683, 0.220809, 0.219014, 0.213015, 0.244470,
2645       0.221042, 0.198996, 0.225295, 0.264466, 0.200600, 0.228143,
2646       0.250138, 0.216784, 0.268002, 0.275734, 0.218144, 0.229866,
2647       0.235443, 0.208909, 0.215067, 0.189743, 0.216741, 0.242686,
2648       0.200234, 0.218882, 0.245991, 0.222815, 0.219576, 0.209773,
2649       0.205247, 0.203855, 0.202534, 0.192536, 0.223447, 0.225810,
2650       0.220583, 0.248421, 0.223424, 0.206033, 0.203130, 0.221225,
2651       0.223763, 0.239216, 0.252311, 0.206893, 0.228461, 0.233591,
2652       0.201512, 0.242382, 0.215147, 0.232578, 0.207116, 0.239213,
2653       0.196215, 0.184367, 0.235135, 0.189768, 0.274084, 0.206267,
2654       0.203283, 0.248828, 0.250214, 0.217465, 0.232080, 0.215150,
2655       0.207936, 0.176789, 0.191338, 0.188655, 0.196181, 0.223691,
2656       0.254257, 0.216874, 0.213536, 0.221399, 0.192024, 0.212534,
2657       0.218169, 0.226635, 0.201191, 0.212700, 0.211634, 0.232353,
2658       0.223636, 0.208605, 0.249132, 0.183681, 0.221842, 0.187136,
2659       0.203772, 0.249575, 0.217713, 0.205193, 0.193941, 0.219707,
2660       0.264226, 0.182105, 0.207183, 0.220845, 0.236571, 0.182390,
2661       0.212914, 0.186266, 0.195361, 0.217665, 0.204527, 0.188804,
2662       0.222832, 0.191193, 0.201073, 0.185616, 0.195011, 0.200828,
2663       0.181434, 0.233707, 0.202925, 0.211992, 0.173100, 0.205258,
2664       0.182695, 0.216520, 0.218422, 0.209490, 0.211257, 0.215801,
2665       0.220652, 0.211409, 0.195731, 0.194957, 0.198888, 0.234926,
2666       0.221377, 0.229822, 0.176700, 0.172322, 0.212265, 0.206133,
2667       0.170355, 0.253305, 0.198688, 0.240043, 0.225384, 0.174729,
2668       0.195820, 0.200093, 0.196912, 0.212308, 0.200490, 0.240415,
2669       0.159744, 0.173686, 0.223853, 0.213604, 0.193779, 0.179328,
2670       0.180873, 0.237481, 0.179640, 0.235857, 0.202847, 0.167690,
2671       0.245093, 0.215504, 0.205848, 0.184408, 0.201626, 0.209651,
2672       0.236306, 0.217803, 0.188534, 0.187861, 0.161663, 0.221718,
2673       0.152044, 0.190412, 0.173505, 0.208995, 0.174575, 0.180734,
2674       0.247704, 0.203388, 0.194184, 0.169679, 0.182703, 0.213402,
2675       0.191808, 0.178604, 0.178116, 0.198452, 0.174687, 0.169809,
2676       0.222851, 0.156811, 0.229170, 0.210359, 0.178557, 0.248570,
2677       0.208536, 0.192571, 0.178912, 0.224505, 0.170822, 0.205492,
2678       0.330973, 0.160924, 0.203724, 0.213255, 0.205827, 0.187162,
2679       0.181252, 0.191723, 0.238106, 0.182398, 0.202358, 0.212066,
2680       0.201255, 0.168159, 0.185219, 0.229176, 0.158222, 0.164092,
2681       0.215405, 0.200724, 0.234811, 0.184488, 0.213112, 0.198577,
2682       0.219622, 0.160671, 0.179349, 0.206681, 0.206091, 0.165618,
2683       0.165203, 0.174442, 0.179287, 0.187318, 0.163481, 0.217440,
2684       0.160381, 0.177025, 0.204385, 0.163676, 0.210733, 0.186519,
2685       0.230701, 0.191764, 0.185119, 0.190770, 0.242987, 0.186812,
2686       0.202906, 0.161935, 0.182426, 0.197922, 0.181475, 0.155903,
2687       0.175006, 0.223482, 0.202706, 0.218108};
2688     wDD = 0.207705;
2689     wND = 0.289628;
2690     wSD = -1;
2691
2692     Mmin = bin[0];
2693     Mmax = bin[nbin];
2694
2695     if(M<Mmin || M>Mmax) return kTRUE;
2696
2697     Int_t ibin=nbin-1;
2698     for(Int_t i=1; i<=nbin; i++)
2699       if(M<=bin[i]) {
2700         ibin=i-1;
2701         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2702         break;
2703       }
2704     wSD=w[ibin];
2705     return kTRUE;
2706   }
2707
2708   return kFALSE;
2709 }
2710
2711 Bool_t AliGenPythia6Diff::IsFromHeavyFlavor(Int_t ipart)
2712 {
2713   // Check if this is a heavy flavor decay product
2714   TParticle *  part = (TParticle *) fParticles.At(ipart);
2715   Int_t mpdg = TMath::Abs(part->GetPdgCode());
2716   Int_t mfl  = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2717   //
2718   // Light hadron
2719   if (mfl >= 4 && mfl < 6) return kTRUE;
2720
2721   Int_t imo = part->GetFirstMother()-1;
2722   TParticle* pm = part;
2723   //
2724   // Heavy flavor hadron produced by generator
2725   while (imo >  -1) {
2726     pm  =  (TParticle*)fParticles.At(imo);
2727     mpdg = TMath::Abs(pm->GetPdgCode());
2728     mfl  = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2729     if ((mfl > 3) && (mfl <6) && mpdg > 400) return kTRUE;
2730     imo = pm->GetFirstMother()-1;
2731   }
2732   return kFALSE;
2733 }
2734
2735 Bool_t AliGenPythia6Diff::CheckDetectorAcceptance(Float_t phi, Float_t eta, Int_t iparticle)
2736 {
2737   // check the eta/phi correspond to the detectors acceptance
2738   // iparticle is the index of the particle to be checked, for PHOS rotation case
2739   if     (fCheckPHOS   && IsInPHOS  (phi,eta,iparticle)) return kTRUE;
2740   else if(fCheckEMCAL  && IsInEMCAL (phi,eta)) return kTRUE;
2741   else if(fCheckBarrel && IsInBarrel(    eta)) return kTRUE;
2742   else                                         return kFALSE;
2743 }
2744
2745 Bool_t AliGenPythia6Diff::TriggerOnSelectedParticles(Int_t np)
2746 {
2747   // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
2748   // implemented primaryly for kPyJets, but extended to any kind of process.
2749
2750   //printf("Check: frag photon %d, pi0 %d, eta %d, electron %d, hadron %d, decay %d, PHOS %d, EMCAL %d, Barrel %d \n",
2751   //       fFragPhotonInCalo,fPi0InCalo, fEtaInCalo,fEleInCalo,fHadronInCalo,fDecayPhotonInCalo,fCheckPHOS,fCheckEMCAL, fCheckBarrel);
2752
2753   Bool_t ok = kFALSE;
2754   for (Int_t i=0; i< np; i++) {
2755
2756     TParticle* iparticle = (TParticle *) fParticles.At(i);
2757
2758     Int_t pdg          = iparticle->GetPdgCode();
2759     Int_t status       = iparticle->GetStatusCode();
2760     Int_t imother      = iparticle->GetFirstMother() - 1;
2761
2762     TParticle* pmother = 0x0;
2763     Int_t momStatus    = -1;
2764     Int_t momPdg       = -1;
2765     if(imother > 0 ){
2766       pmother = (TParticle *) fParticles.At(imother);
2767       momStatus    = pmother->GetStatusCode();
2768       momPdg       = pmother->GetPdgCode();
2769     }
2770
2771     ok = kFALSE;
2772
2773     //
2774     // Check the particle type: hadron (not pi0 or eta), electron, decay photon (from pi0 or eta or any), pi0 or eta
2775     //
2776     // Hadron
2777     if (fHadronInCalo && status == 1)
2778     {
2779       if(TMath::Abs(pdg) > 23 && pdg !=221 && pdg != 111) // avoid photons, electrons, muons, neutrinos and eta or pi0
2780         // (in case neutral mesons were declared stable)
2781         ok = kTRUE;
2782     }
2783     //Electron
2784     else if (fEleInCalo && status == 1 && TMath::Abs(pdg) == 11)
2785     {
2786       ok = kTRUE;
2787     }
2788     //Fragmentation photon
2789     else if (fFragPhotonInCalo && pdg == 22 && status == 1)
2790     {
2791       if(momStatus != 11) ok = kTRUE ;  // No photon from hadron decay
2792     }
2793     // Decay photon
2794     else if (fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay && pdg == 22) // pi0 status can be 1 or 11 depending on decay settings, work only for 11
2795     {
2796       if( momStatus == 11)
2797       {
2798         //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Decay photon! pdg %d, status %d, pt %2.2f, mom: pdg %d, pt %2.2f\n",
2799         //                                                   pdg,status,iparticle->Pt(),momPdg,pmother->Pt());
2800         ok = kTRUE ;  // photon from hadron decay
2801
2802         //In case only decays from pi0 or eta requested
2803         if(fPi0InCalo && momPdg!=111) ok = kFALSE;
2804         if(fEtaInCalo && momPdg!=221) ok = kFALSE;
2805       }
2806
2807     }
2808     // Pi0 or Eta particle
2809     else if ((fPi0InCalo || fEtaInCalo))
2810     {
2811       if(fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay ) continue ;
2812
2813       if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
2814       {
2815         //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2816         ok = kTRUE;
2817       }
2818       else if (fEtaInCalo && pdg == 221)
2819       {
2820         //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Eta! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2821         ok = kTRUE;
2822       }
2823
2824     }// pi0 or eta
2825
2826     //
2827     // Check that the selected particle is in the calorimeter acceptance
2828     //
2829     if(ok && iparticle->Pt() > fTriggerParticleMinPt)
2830     {
2831       //Just check if the selected particle falls in the acceptance
2832       if(!fForceNeutralMeson2PhotonDecay )
2833       {
2834         //printf("\t Check acceptance! \n");
2835         Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2836         Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2837
2838         if(CheckDetectorAcceptance(phi,eta,i))
2839         {
2840           ok =kTRUE;
2841           AliDebug(1,Form("Selected trigger pdg %d, status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2842           //printf("\t Accept \n");
2843           break;
2844         }
2845         else ok = kFALSE;
2846       }
2847       //Mesons have several decay modes, select only those decaying into 2 photons
2848       else if(fForceNeutralMeson2PhotonDecay && (fPi0InCalo || fEtaInCalo))
2849       {
2850         // In case we want the pi0/eta trigger,
2851         // check the decay mode (2 photons)
2852
2853         //printf("\t Force decay 2 gamma\n");
2854
2855         Int_t ndaughters = iparticle->GetNDaughters();
2856         if(ndaughters != 2){
2857           ok=kFALSE;
2858           continue;
2859         }
2860
2861         TParticle*d1 = (TParticle *) fParticles.At(iparticle->GetDaughter(0)-1);
2862         TParticle*d2 = (TParticle *) fParticles.At(iparticle->GetDaughter(1)-1);
2863         if(!d1 || !d2) {
2864           ok=kFALSE;
2865           continue;
2866         }
2867
2868         //iparticle->Print();
2869         //d1->Print();
2870         //d2->Print();
2871
2872         Int_t pdgD1 = d1->GetPdgCode();
2873         Int_t pdgD2 = d2->GetPdgCode();
2874         //printf("\t \t 2 daughters pdg = %d - %d\n",pdgD1,pdgD2);
2875         //printf("mother %d - %d\n",d1->GetFirstMother(),d2->GetFirstMother());
2876
2877         if(pdgD1 != 22  || pdgD2 != 22){
2878           ok = kFALSE;
2879           continue;
2880         }
2881
2882         //printf("\t accept decay\n");
2883
2884         //Trigger on the meson, not on the daughter
2885         if(!fDecayPhotonInCalo){
2886
2887           Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2888           Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2889
2890           if(CheckDetectorAcceptance(phi,eta,i))
2891           {
2892             //printf("\t Accept meson pdg %d\n",pdg);
2893             ok =kTRUE;
2894             AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2895             break;
2896           } else {
2897             ok=kFALSE;
2898             continue;
2899           }
2900         }
2901
2902         //printf("Check daughters acceptance\n");
2903
2904         //Trigger on the meson daughters
2905         //Photon 1
2906         Float_t phi = d1->Phi()*180./TMath::Pi(); //Convert to degrees
2907         Float_t eta =TMath::Abs(d1->Eta()); //in calos etamin=-etamax
2908         if(d1->Pt() > fTriggerParticleMinPt)
2909         {
2910           //printf("\t Check acceptance photon 1! \n");
2911           if(CheckDetectorAcceptance(phi,eta,i))
2912           {
2913             //printf("\t Accept Photon 1\n");
2914             ok =kTRUE;
2915             AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2916             break;
2917           }
2918           else ok = kFALSE;
2919         } // pt cut
2920         else  ok = kFALSE;
2921
2922         //Photon 2
2923         phi = d2->Phi()*180./TMath::Pi(); //Convert to degrees
2924         eta =TMath::Abs(d2->Eta()); //in calos etamin=-etamax
2925
2926         if(d2->Pt() > fTriggerParticleMinPt)
2927         {
2928           //printf("\t Check acceptance photon 2! \n");
2929           if(CheckDetectorAcceptance(phi,eta,i))
2930           {
2931             //printf("\t Accept Photon 2\n");
2932             ok =kTRUE;
2933             AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2934             break;
2935           }
2936           else ok = kFALSE;
2937         } // pt cut
2938         else ok = kFALSE;
2939       } // force 2 photon daughters in pi0/eta decays
2940       else ok = kFALSE;
2941     } else ok = kFALSE; // check acceptance
2942   } // primary loop
2943
2944   //
2945   // If requested, rotate the particles event in phi to enhance/speed PHOS selection
2946   // A particle passing all trigger conditions except phi position in PHOS, is used as reference
2947   //
2948   if(fCheckPHOSeta)
2949   {
2950     RotatePhi(ok);
2951   }
2952
2953   return ok;
2954 }
2955
2956
2957 AliGenerator* GenPythia6Diff()
2958 {
2959   AliGenPythia6Diff* pythia = new AliGenPythia6Diff(-1);
2960   pythia->SetMomentumRange(0, 999999.);
2961   pythia->SetThetaRange(0., 180.);
2962   pythia->SetYRange(-12.,12.);
2963   pythia->SetPtRange(0,1000.);
2964   pythia->SetProcess(kPyMb);
2965   pythia->SetEnergyCMS(VAR_PYTHIA6_CMS_ENERGY);
2966   //    Tune
2967   //    320     Perugia 0
2968   pythia->SetTune(320);
2969   pythia->UseNewMultipleInteractionsScenario();
2970   pythia->SetCrossingAngle(0,0.000280);
2971   pythia->SetTuneForDiff();
2972
2973   return pythia;
2974 }