Updates
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 12 Oct 2012 19:14:48 +0000 (19:14 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 12 Oct 2012 19:14:48 +0000 (19:14 +0000)
Chiara Oppedisano

DPMJET/dpmjet3.0-5.f
TDPMjet/AliGenDPMjet.cxx
TDPMjet/AliGenDPMjet.h
TDPMjet/DPMcommon.h
TDPMjet/TDPMjet.cxx
TDPMjet/fastGenDPMjet.C

index 027d9b8..3c0b6a9 100644 (file)
@@ -5178,7 +5178,8 @@ C           ENDIF
      &                INTER1(MAXINT),INTER2(MAXINT)
 * Glauber formalism: collision properties
       COMMON /DTGLCP/ RPROJ,RTARG,BIMPAC,
-     &                NWTSAM,NWASAM,NWBSAM,NWTACC,NWAACC,NWBACC
+     &                NWTSAM,NWASAM,NWBSAM,NWTACC,NWAACC,NWBACC,
+     &                NCP,NCT
 * central particle production, impact parameter biasing
       COMMON /DTIMPA/ BIMIN,BIMAX,XSFRAC,ICENTR
 **temporary
@@ -5191,6 +5192,8 @@ C           ENDIF
       IREJ   = 0
       ICREQU = ICREQU+1
       NC     = 0
+      NCP    = 0
+      NCT    = 0
 
     1 CONTINUE
       ICSAMP = ICSAMP+1
@@ -5230,6 +5233,14 @@ C           ENDIF
          ITOLD  = IT
          JJPOLD = JJPROJ
          EPROLD = EPROJ
+        DO 8 I=1, IP
+           NCP = NCP+JSSH(I)
+*          WRITE(6,*)' PROJ.NUCL. ',I,' NCOLL = ',NCP 
+    8 CONTINUE
+        DO 9 I=1, IT
+           NCT = NCT+JTSH(I)
+*          WRITE(6,*)' TAR.NUCL. ',I,' NCOLL = ',NCT 
+    9 CONTINUE
       ENDIF
 
 * force diffractive particle production in h-K interactions
index 3fd2c40..78a3dbf 100644 (file)
@@ -45,9 +45,9 @@ ClassImp(AliGenDPMjet)
 //______________________________________________________________________________
 AliGenDPMjet::AliGenDPMjet()
     :AliGenMC(), 
-     fBeamEn(2750.),
+     fBeamEn(0.),
      fMinImpactParam(0.),
-     fMaxImpactParam(5.),
+     fMaxImpactParam(100.),
      fICentr(0),
      fSelectAll(0),
      fFlavor(0),
@@ -67,7 +67,9 @@ AliGenDPMjet::AliGenDPMjet()
      fTriggerMultiplicityPtMin(0),
      fkTuneForDiff(0),
      fProcDiff(0),
-     fFragmentation(kFALSE)
+     fFragmentation(kFALSE),
+     fHeader(AliGenDPMjetEventHeader("DPMJET"))
+
 {
 // Constructor
     fEnergyCMS = 5500.;
@@ -78,9 +80,9 @@ AliGenDPMjet::AliGenDPMjet()
 //______________________________________________________________________________
 AliGenDPMjet::AliGenDPMjet(Int_t npart)
     :AliGenMC(npart),
-     fBeamEn(2750.),
+     fBeamEn(0.),
      fMinImpactParam(0.),
-     fMaxImpactParam(5.),
+     fMaxImpactParam(100.),
      fICentr(0),
      fSelectAll(0),
      fFlavor(0),
@@ -100,7 +102,9 @@ AliGenDPMjet::AliGenDPMjet(Int_t npart)
      fTriggerMultiplicityPtMin(0),
      fkTuneForDiff(0),
      fProcDiff(0),
-     fFragmentation(kFALSE)
+     fFragmentation(kFALSE),
+     fHeader(AliGenDPMjetEventHeader("DPMJET"))
+
 
 {
 // Default PbPb collisions at 5. 5 TeV
@@ -116,9 +120,9 @@ AliGenDPMjet::AliGenDPMjet(Int_t npart)
 
 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
     :AliGenMC(),
-     fBeamEn(2750.),
+     fBeamEn(0.),
      fMinImpactParam(0.),
-     fMaxImpactParam(5.),
+     fMaxImpactParam(100.),
      fICentr(0),
      fSelectAll(0),
      fFlavor(0),
@@ -138,7 +142,9 @@ AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
      fTriggerMultiplicityPtMin(0),
      fkTuneForDiff(0),
      fProcDiff(0),
-     fFragmentation(kFALSE)
+     fFragmentation(kFALSE),
+     fHeader(0x0)
+
 
 {
     // Dummy copy constructor
@@ -155,26 +161,25 @@ void AliGenDPMjet::Init()
 {
 // Initialization
     
+    if(fEnergyCMS>0. && fBeamEn<0.1) fBeamEn = fEnergyCMS/2;
     SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget, 
                      fBeamEn,fEnergyCMS));
 
     fDPMjet=(TDPMjet*) fMCEvGen;
-    //
+    if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
+    
     // **** Flag to force central production
     // fICentr=1. central production forced 
     // fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam       
     // fICentr<-99 -> fraction of x-sec. = XSFRAC                
     // fICentr=-1. -> evaporation/fzc suppressed                 
-    // fICentr<-1. -> evaporation/fzc allowed            
-    if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
-    
+    // fICentr<-1. -> evaporation/fzc allowed                
     fDPMjet->SetfFCentr(fICentr);  
+    
     fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam); 
     fDPMjet->SetPi0Decay(fPi0Decay);
     fDPMjet->SetDecayAll(fDecayAll);
     fDPMjet->SetFragmentProd(fFragmentation);
-
-    AliGenMC::Init();
     
 //
 //  Initialize DPMjet  
@@ -462,40 +467,40 @@ Bool_t AliGenDPMjet::Stable(TParticle*  particle)
 //______________________________________________________________________________
 void AliGenDPMjet::MakeHeader()
 {
-//  printf("MakeHeader %13.3f \n", fDPMjet->GetBImpac());
 // Builds the event header, to be called after each event
-    AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
-    ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
-    ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
-    ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
-    ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetProjParticipants(), 
-                                                        fDPMjet->GetTargParticipants());
-
-    if(fProcDiff>0){
-    ((AliGenDPMjetEventHeader*) header)->SetProcessType(fProcDiff);
-    }
-    else 
-      ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
+    fHeader.SetNProduced(fDPMjet->GetNumStablePc());
+    fHeader.SetImpactParameter(fDPMjet->GetBImpac());
+    fHeader.SetTotalEnergy(fDPMjet->GetTotEnergy());
+    fHeader.SetParticipants(fDPMjet->GetProjParticipants(), 
+                           fDPMjet->GetTargParticipants());
+
+    fHeader.SetCollisions(DTGLCP.ncp, DTGLCP.nct, 
+       fDPMjet->GetProjWounded(),fDPMjet->GetTargWounded());
+               
+    if(fProcDiff>0)  fHeader.SetProcessType(fProcDiff);
+    else fHeader.SetProcessType(fDPMjet->GetProcessCode());
 
     // Bookkeeping for kinematic bias
-    ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
+    fHeader.SetTrials(fTrials);
     // Event Vertex
-    header->SetPrimaryVertex(fVertex);
-    header->SetInteractionTime(fTime);
-    ((AliGenDPMjetEventHeader*) header)->SetNDiffractive(POEVT1.nsd1, POEVT1.nsd2, POEVT1.ndd);
-    gAlice->SetGenEventHeader(header);    
-    AddHeader(header);
+    fHeader.SetPrimaryVertex(fVertex);
+    fHeader.SetInteractionTime(fTime);
+    fHeader.SetNDiffractive(POEVT1.nsd1, POEVT1.nsd2, POEVT1.ndd);
+//    gAlice->SetGenEventHeader(fHeader);    
+    AddHeader(&fHeader);
+    fCollisionGeometry = &fHeader;
 }
 
-void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
+//______________________________________________________________________________
+/*void AliGenDPMjet::AddHeader(AliGenEventHeader* fHeader)
 {
-    // Add header to container or runloader
+    // Add fHeader to container or runloader
     if (fContainer) {
-        fContainer->AddHeader(header);
+        fContainer->AddHeader(fHeader);
     } else {
-        AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
+        AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(fHeader);
     }
-}
+}*/
 
 
 //______________________________________________________________________________
index 380eec2..0b5a786 100644 (file)
@@ -8,6 +8,7 @@
 
 #include "AliGenMC.h"
 #include "TDPMjet.h"
+#include "AliGenDPMjetEventHeader.h"
 #include <TString.h>
 #include <TArrayI.h>
 
@@ -34,8 +35,8 @@ class AliGenDPMjet : public AliGenMC
     virtual void    Generate();
     virtual void    Init();
     virtual void    FinishRun();
-    virtual void    SetEnergyCMS(Float_t energy = 14000.) {fEnergyCMS = energy; fBeamEn = energy / 2.;}
-    virtual void    SetpBeamEnergy(Float_t benergy = 14000.) {fBeamEn = benergy;}
+    virtual void    SetEnergyCMS(Float_t energy = 14000.) {fEnergyCMS = energy;}
+    virtual void    SetProjectileBeamEnergy(Float_t benergy = 7000.) {fBeamEn = benergy;}
     virtual void    SetImpactParameterRange(Float_t bmin=0., Float_t bmax=1.)
                        {fMinImpactParam=bmin; fMaxImpactParam=bmax;}
     virtual void    SetProcess(DpmProcess_t iproc) {fProcess = iproc;}
@@ -59,15 +60,16 @@ class AliGenDPMjet : public AliGenMC
       fTriggerMultiplicityPtMin = ptmin;}
 
     AliGenDPMjet &  operator=(const AliGenDPMjet & rhs);
-    void     AddHeader(AliGenEventHeader* header);
+    //void     AddHeader(AliGenEventHeader* header);
 
    void SetTuneForDiff(Bool_t a=kTRUE) {fkTuneForDiff=a;}
 
    virtual void  SetFragmentProd(Bool_t val) {fFragmentation = val;}
+   virtual Bool_t ProvidesCollisionGeometry() const {return kTRUE;}
 
  protected:
-    Bool_t SelectFlavor(Int_t pid);
-    void   MakeHeader();
+   Bool_t SelectFlavor(Int_t pid);
+   void   MakeHeader();
 
  protected:
     Float_t       fBeamEn;        // beam energy
@@ -96,6 +98,7 @@ class AliGenDPMjet : public AliGenMC
     Int_t  fProcDiff;
     
     Bool_t fFragmentation; // Allows evaporation and fragments production
+    AliGenDPMjetEventHeader fHeader; // MC header
 
  private:
     // adjust the weight from kinematic cuts
@@ -109,7 +112,7 @@ class AliGenDPMjet : public AliGenMC
    Bool_t GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax, 
                                               Double_t &wSD, Double_t &wDD, Double_t &wND);
 
-    ClassDef(AliGenDPMjet,5) // AliGenerator interface to DPMJET
+    ClassDef(AliGenDPMjet,6) // AliGenerator interface to DPMJET
 };
 #endif
 
index aaff398..ede9b27 100644 (file)
@@ -200,6 +200,8 @@ typedef struct {
    Int_t       nwtacc;
    Int_t       nwtaac;
    Int_t       nwtbac;
+   Int_t        ncp;
+   Int_t        nct;
 } DtglcpCommon;
 
 #define DTGLCP COMMON_BLOCK(DTGLCP,dtglcp)
index ab1e95b..536e15b 100644 (file)
@@ -141,7 +141,7 @@ TDPMjet::TDPMjet(DpmProcess_t  iproc, Int_t Ip=208, Int_t Ipz=82, Int_t It=208,
       fPi0Decay(0),
       fDecayAll(0),
       fProcess(iproc),
-      fFragmentation(-1)
+      fFragmentation(kFALSE)
 {  
     printf("TDPMJet Constructor %d %d %d %d \n", Ip, Ipz, It, Itz);
 }
@@ -268,6 +268,8 @@ void TDPMjet::Initialize()
 //
 //  Write standard DPMJET input cards 
 //
+    if(fFragmentation) printf("\tTDPMJet fragmentation/evaporation applied\n");
+
     FILE* out = fopen("dpmjet.inp","w");
 //  Projectile and Target definition 
     if (fIp == 1 && fIpz ==1) {
@@ -287,7 +289,7 @@ void TDPMjet::Initialize()
     }
 
 //  Beam energy and crossing-angle
-    fprintf(out, "CMENERGY      %10.1f\n",fCMEn);      
+    fprintf(out, "CMENERGY      %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fCMEn, 0., 0., 0., 0., 0.);      
     if(fIt == 1 && fIp ==1){ 
       fprintf(out, "BEAM      %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fEpn, fEpn, 0., 0., 0., 0.); //p-p
     }
@@ -295,11 +297,11 @@ void TDPMjet::Initialize()
       if(fIp>1 && fIt>1) fprintf(out, "BEAM      %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fEpn, fEpn, 0., 0., 0., 0.);//A-A
       else if(fIp==1 && fIt>1){ // proton towards A side (directed z>0)
         fprintf(out, "BEAM      %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", fEpn,fEpn*fItz/fIt, 0., 0., 0., 0.);//pA
-       printf("\n  TDPMjet::Initialize() -> p-A: p beam energy =  %10.1f, CMS energy = %10.1f\n\n",fEpn,fCMEn/2);
+       printf("\n  TDPMjet::Initialize() -> p-A: projectile (p) energy =  %10.1f, CMS energy = %10.1f\n\n",fEpn,fCMEn);
       }
       else if(fIt==1 && fIp>1){ // proton towards C side (directed z<0)
-        fprintf(out, "BEAM      %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", fEpn*fIpz/fIp, fEpn, 0., 0., 0., 0.);//A-p
-       printf("\n  TDPMjet::Initialize() -> A-p: p beam energy =  %10.1f, CMS energy = %10.1f\n\n",fEpn,fCMEn/2);
+        fprintf(out, "BEAM      %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", fEpn, fEpn*fIp/fIpz, 0., 0., 0., 0.);//A-p
+       printf("\n  TDPMjet::Initialize() -> A-p: projectile (A) energy =  %10.1f, CMS energy = %10.1f\n\n",fEpn,fCMEn);
       }
     }
 //  Centrality
index 01b2b39..919de38 100644 (file)
@@ -1,14 +1,18 @@
 AliGenerator*  CreateGenerator();
 
-void fastGenDPMjet(Int_t nev = 1, char* filename = "galice.root")
+void fastGenDPMjet(Int_t nev = 1, char* filename = "dpmjet.root")
 {
 //  Runloader
+#if defined(__CINT__)
     gSystem->Load("liblhapdf"); 
     gSystem->Load("libpythia6");     
     gSystem->Load("libdpmjet.so");
     gSystem->Load("libTDPMjet.so");
+#endif
+
+    AliPDG::AddParticlesToPdgDataBase();
     
-    AliRunLoader* rl = AliRunLoader::Open("galice.root","FASTRUN","recreate");
+    AliRunLoader* rl = AliRunLoader::Open(filename,"FASTRUN","recreate");
     AliPythiaRndm::SetPythiaRandom(new TRandom());
 
     rl->SetCompressionLevel(2);
@@ -17,49 +21,42 @@ void fastGenDPMjet(Int_t nev = 1, char* filename = "galice.root")
     rl->MakeTree("E");
     gAlice->SetRunLoader(rl);
 
-//  Create stack
+    //  Create stack
     rl->MakeStack();
-    AliStack* stack      = rl->Stack();
+    AliStack* stack   = rl->Stack();
  
-//  Header
+    //  Header
     AliHeader* header = rl->GetHeader();
-//
-//  Create and Initialize Generator
-    AliGenerator *gener = CreateGenerator();
+    //AliGenDPMjetEventHEader* dpmHeader;
+
+    //  Create and Initialize Generator
+    AliGenDPMjet *gener = CreateGenerator();
+    //AliCollisionGeometry *coll = gener->CollisionGeometry();
     gener->Init();
     gener->SetStack(stack);
     
-//
-//                        Event Loop
-//
+    // Event Loop
     Int_t iev;
      
     for (iev = 0; iev < nev; iev++) {
 
-       printf("\n \n Event number %d \n \n", iev);
+       if(iev%500==0) printf("\n  Event number %d  \n", iev);
        
-//  Initialize event
+       //  Initialize event
        header->Reset(0,iev);
        rl->SetEventNumber(iev);
        stack->Reset();
        rl->MakeTree("K");
-//     stack->ConnectTree();
     
-//  Generate event
+       //  Generate event
        gener->Generate();
-//  Analysis
-       Int_t npart = stack->GetNprimary();
-       printf("Analyse %d Particles\n", npart);
-       for (Int_t part=0; part<npart; part++) {
-           TParticle *MPart = stack->Particle(part);
-           Int_t pdg  = MPart->GetPdgCode();
-           Int_t pis  = MPart->GetStatusCode();
-//         printf("Particle %5d %5d\n", part, pdg, pis);
-       }
        
-//  Finish event
+       Int_t npart = stack->GetNprimary();
+        
+       //  Finish event
        header->SetNprimary(stack->GetNprimary());
        header->SetNtrack(stack->GetNtrack());  
+
 //      I/O
 //     
        stack->FinishEvent();
@@ -68,11 +65,9 @@ void fastGenDPMjet(Int_t nev = 1, char* filename = "galice.root")
        rl->WriteKinematics("OVERWRITE");
 
     } // event loop
-//
-//                         Termination
-//  Generator
+
     gener->FinishRun();
-//  Write file
+
     rl->WriteHeader("OVERWRITE");
     gener->Write();
     rl->Write();
@@ -82,17 +77,36 @@ void fastGenDPMjet(Int_t nev = 1, char* filename = "galice.root")
 
 AliGenerator*  CreateGenerator()
 { 
-    //  kDpmMb
-    //  kDpmMbNonDiffr
-    //  kDpmDiffr
-    //  kDpmSingleDiffr
-    //  kDpmDoubleDiffr
-    gener = new AliGenDPMjet(1);
-    gener->SetProcess(kDpmMb);
-    gener->SetProjectile("P", 1, 1);
-    gener->SetTarget("P", 1, 1);
-    gener->SetEnergyCMS(14000.);
-    return gener;
+  AliGenDPMjet *gener = new AliGenDPMjet(-1);
+  
+  //  p-A
+  /*Float_t pEnergy = 4000.;
+  gener->SetProjectile("P", 1, 1);
+  gener->SetTarget("A", 208, 82);
+  //
+  gener->SetEnergyCMS(TMath::Sqrt(82./208.) * 2* pEnergy);
+  gener->SetProjectileBeamEnergy(pEnergy); */ 
+  
+  //  A-p
+  Float_t pEnergy = 1577.;
+  gener->SetProjectile("A", 208, 82);
+  gener->SetTarget("P", 1, 1);
+  //
+  gener->SetEnergyCMS(TMath::Sqrt(208./82.) * 2* pEnergy);
+  gener->SetProjectileBeamEnergy(pEnergy);
+
+  // Pb-Pb
+  /*  gener->SetProjectile("A", 208, 82);
+  gener->SetTarget("A", 208, 82);
+  gener->SetEnergyCMS(2760.);
+  //
+  gener->SetProcess(kDpmMb);
+  gener->SetImpactParameterRange(0., 100.);
+  //gener->SetFragmentProd(kTRUE);*/
+
+  gener->SetTrackingFlag(0);
+  
+  return gener;
 }