Finish run with print-out of run statistics added.
[u/mrichter/AliRoot.git] / TDPMjet / AliGenDPMjet.cxx
index 313638e..4692ccf 100644 (file)
 #include <TParticleClassPDG.h>
 #include <TPDGCode.h>
 #include <TLorentzVector.h>
-
+#include <TClonesArray.h>
+#include "AliRunLoader.h"
 #include "AliGenDPMjet.h"
 #include "AliGenDPMjetEventHeader.h"
-#include "AliPythia.h"
 #include "AliRun.h"
 #include "AliDpmJetRndm.h"
+#include "AliHeader.h"
+#include "AliStack.h"
+#include "AliMC.h"
 
- ClassImp(AliGenDPMjet)
-
+ClassImp(AliGenDPMjet)
 
 //______________________________________________________________________________
 AliGenDPMjet::AliGenDPMjet()
-                 :AliGenMC()
+    :AliGenMC(), 
+     fBeamEn(2750.),
+     fEnergyCMS(5500.),
+     fMinImpactParam(0.),
+     fMaxImpactParam(5.),
+     fICentr(0),
+     fSelectAll(0),
+     fFlavor(0),
+     fTrials(0),
+     fSpectators(1),
+     fSpecn(0),
+     fSpecp(0),
+     fDPMjet(0),
+     fNoGammas(0),
+     fLHC(0),
+     fPi0Decay(0),
+     fGenImpPar(0.),
+     fProcess(kDpmMb)
 {
 // Constructor
-    fParticles = 0;
-    fDPMjet = 0;
     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
 }
 
 
 //______________________________________________________________________________
 AliGenDPMjet::AliGenDPMjet(Int_t npart)
-    :AliGenMC(npart)
+    :AliGenMC(npart),
+     fBeamEn(2750.),
+     fEnergyCMS(5500.),
+     fMinImpactParam(0.),
+     fMaxImpactParam(5.),
+     fICentr(0),
+     fSelectAll(0),
+     fFlavor(0),
+     fTrials(0),
+     fSpectators(1),
+     fSpecn(0),
+     fSpecp(0),
+     fDPMjet(0),
+     fNoGammas(0),
+     fLHC(0),
+     fPi0Decay(0),
+     fGenImpPar(0.),
+     fProcess(kDpmMb)
 {
 // Default PbPb collisions at 5. 5 TeV
 //
     fName = "DPMJET";
     fTitle= "Particle Generator using DPMJET";
-
-    SetEnergyCMS();
     SetTarget();
     SetProjectile();
-    SetCentral();
-    SetImpactParameterRange();
-    SetBoostLHC();
-    
-    fKeep       =  0;
-    fDecaysOff  =  1;
-    fEvaluate   =  0;
-    fSelectAll  =  0;
-    fFlavor     =  0;
-    fSpectators =  1;
     fVertex.Set(3);
-        
-    fParticles = new TClonesArray("TParticle",10000);    
-
-    // Set random number generator   
-    fDPMjet = 0;
-    // Instance AliPythia
-    AliPythia::Instance(); 
     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
 }
 
 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
-:AliGenMC()
+    :AliGenMC(),
+     fBeamEn(2750.),
+     fEnergyCMS(5500.),
+     fMinImpactParam(0.),
+     fMaxImpactParam(5.),
+     fICentr(0),
+     fSelectAll(0),
+     fFlavor(0),
+     fTrials(0),
+     fSpectators(1),
+     fSpecn(0),
+     fSpecp(0),
+     fDPMjet(0),
+     fNoGammas(0),
+     fLHC(0),
+     fPi0Decay(0),
+     fGenImpPar(0.),
+     fProcess(kDpmMb)
 {
+    // Dummy copy constructor
 }
 
 //______________________________________________________________________________
 AliGenDPMjet::~AliGenDPMjet()
 {
 // Destructor
-    delete fParticles;
 }
-
-
 //______________________________________________________________________________
 void AliGenDPMjet::Init()
 {
@@ -115,7 +145,7 @@ void AliGenDPMjet::Init()
     
     fDPMjet->SetfFCentr(fICentr);  
     fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam); 
-    
+    fDPMjet->SetPi0Decay(fPi0Decay);
 //
 //  Initialize DPMjet  
 //    
@@ -130,7 +160,6 @@ void AliGenDPMjet::Generate()
 
   Float_t polar[3]    =   {0,0,0};
   Float_t origin[3]   =   {0,0,0};
-  Float_t origin0[3]  =   {0,0,0};
   Float_t p[3];
   Float_t tof;
 
@@ -154,13 +183,13 @@ void AliGenDPMjet::Generate()
       fDPMjet->GenerateEvent();
       fTrials++;
 
-      fDPMjet->ImportParticles(fParticles,"All");      
+      fDPMjet->ImportParticles(&fParticles,"All");      
       if (fLHC) Boost();
 
       // Temporaneo
       fGenImpPar = fDPMjet->GetBImpac();
       
-      Int_t np = fParticles->GetEntriesFast();
+      Int_t np = fParticles.GetEntriesFast();
       printf("\n **************************************************%d\n",np);
       Int_t nc=0;
       if (np==0) continue;
@@ -173,16 +202,10 @@ void AliGenDPMjet::Generate()
          pSelected[i] = 0;
       }
       
-//      Get event vertex
-      
-      fVertex[0] = origin0[0];
-      fVertex[1] = origin0[1]; 
-      fVertex[2] = origin0[2];
-      
 //      First select parent particles
 
       for (i = 0; i<np; i++) {
-         TParticle *iparticle = (TParticle *) fParticles->At(i);
+         TParticle *iparticle = (TParticle *) fParticles.At(i);
 
 // Is this a parent particle ?
 
@@ -214,7 +237,7 @@ void AliGenDPMjet::Generate()
 
 
       for (i=0; i<np; i++) {
-         TParticle *iparticle = (TParticle *) fParticles->At(i);
+         TParticle *iparticle = (TParticle *) fParticles.At(i);
 
 // Is this a final state particle ?
 
@@ -248,7 +271,7 @@ void AliGenDPMjet::Generate()
 // Write particles to stack
 
       for (i = 0; i<np; i++) {
-         TParticle *  iparticle = (TParticle *) fParticles->At(i);
+         TParticle *  iparticle = (TParticle *) fParticles.At(i);
          Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
          if (pSelected[i]) {
              
@@ -268,8 +291,8 @@ void AliGenDPMjet::Generate()
              TParticle* mother = 0;
              if (hasMother) {
                  imo = iparticle->GetFirstMother();
-                 mother = (TParticle *) fParticles->At(imo);
-                 imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
+                 mother = (TParticle *) fParticles.At(imo);
+                 imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
              } // if has mother   
 
              Bool_t tFlag = (fTrackIt && (ks == 1));
@@ -292,81 +315,6 @@ void AliGenDPMjet::Generate()
   SetHighWaterMark(nt);
 }
 
-
-//______________________________________________________________________________
-void AliGenDPMjet::KeepFullEvent()
-{
-    fKeep=1;
-}
-
-
-//______________________________________________________________________________
-/*void AliGenDPMjet::EvaluateCrossSections()
-{
-//     Glauber Calculation of geometrical x-section
-//
-    Float_t xTot       = 0.;          // barn
-    Float_t xTotHard   = 0.;          // barn 
-    Float_t xPart      = 0.;          // barn
-    Float_t xPartHard  = 0.;          // barn 
-    Float_t sigmaHard  = 0.1;         // mbarn
-    Float_t bMin       = 0.;
-    Float_t bMax       = fDPMjet->GetProjRadius()+fDPMjet->GetTargRadius();
-    const Float_t kdib = 0.2;
-    Int_t   kMax       = Int_t((bMax-bMin)/kdib)+1;
-
-
-    printf("\n Projectile Radius (fm): %f \n",fDPMjet->GetProjRadius());
-    printf("\n Target     Radius (fm): %f \n",fDPMjet->GetTargRadius());    
-    Int_t i;
-    Float_t oldvalue= 0.;
-
-    Float_t* b   = new Float_t[kMax];
-    Float_t* si1 = new Float_t[kMax];    
-    Float_t* si2 = new Float_t[kMax];    
-    
-    for (i = 0; i < kMax; i++)
-    {
-       Float_t xb  = bMin+i*kdib;
-       Float_t ov;
-       ov=fDPMjet->Profile(xb);
-       // ATT!->Manca la x-sec anel. nucleone-nucleone
-       Float_t gb  =  2.*0.01*fDPMjet->TMath::Pi()*kdib*xb*(1.-TMath::Exp(-fDPMjet->GetXSFrac()*ov));
-       Float_t gbh =  2.*0.01*fDPMjet->TMath::Pi()*kdib*xb*sigmaHard*ov;
-       xTot+=gb;
-       xTotHard += gbh;
-       if (xb > fMinImpactParam && xb < fMaxImpactParam)
-       {
-           xPart += gb;
-           xPartHard += gbh;
-       }
-       
-       if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
-       oldvalue = xTot;
-       printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
-       printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
-       if (i>0) {
-           si1[i] = gb/kdib;
-           si2[i] = gbh/gb;
-           b[i]  = xb;
-       }
-    }
-
-    printf("\n Total cross section (barn): %f \n",xTot);
-    printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
-    printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
-    printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
-
-//  Store result as a graph
-    b[0] = 0;
-    si1[0] = 0;
-    si2[0]=si2[1];
-    
-    fDsigmaDb  = new TGraph(i, b, si1);
-    fDnDb      = new TGraph(i, b, si2);
-}*/
-
-
 //______________________________________________________________________________
 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
 {
@@ -383,7 +331,7 @@ Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
        imin = iparticle->GetFirstDaughter();
        imax = iparticle->GetLastDaughter();       
        for (i = imin; i <= imax; i++){
-           TParticle *  jparticle = (TParticle *) fParticles->At(i);   
+           TParticle *  jparticle = (TParticle *) fParticles.At(i);    
            Int_t ip = jparticle->GetPdgCode();
            if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
                selected=kTRUE; break;
@@ -444,14 +392,24 @@ void AliGenDPMjet::MakeHeader()
     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(), 
                                                         fDPMjet->GetfIt());
-
+ ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
 // Bookkeeping for kinematic bias
     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
 // Event Vertex
     header->SetPrimaryVertex(fVertex);
     gAlice->SetGenEventHeader(header);    
+    AddHeader(header);
 }
 
+void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
+{
+    // Add header to container or runloader
+    if (fContainer) {
+        fContainer->AddHeader(header);
+    } else {
+        AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
+    }
+}
 
 
 //______________________________________________________________________________
@@ -462,5 +420,12 @@ AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
 }
 
 
+void AliGenDPMjet::FinishRun()
+{
+    // Print run statistics
+    fDPMjet->Dt_Dtuout();
+}
+
+    
 
 //______________________________________________________________________________