Code for the simulations at generator level, used for model comparison in the first...
authormfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 21 Apr 2010 13:26:22 +0000 (13:26 +0000)
committermfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 21 Apr 2010 13:26:22 +0000 (13:26 +0000)
Initial commit.

PWG0/PWG0baseLinkDef.h
PWG0/genLevelSimulation/AliAnalysisTaskdNdetaMC.cxx [new file with mode: 0644]
PWG0/genLevelSimulation/AliAnalysisTaskdNdetaMC.h [new file with mode: 0644]
PWG0/genLevelSimulation/MergeCollectionFromGrid.C [new file with mode: 0644]
PWG0/genLevelSimulation/fastGen.C [new file with mode: 0644]
PWG0/genLevelSimulation/runTask.C [new file with mode: 0644]
PWG0/genLevelSimulation/runTaskProof.C [new file with mode: 0644]
PWG0/genLevelSimulation/rungen.C [new file with mode: 0644]
PWG0/libPWG0base.pkg

index 993c4c7..104f75e 100644 (file)
@@ -22,4 +22,6 @@
 
 #pragma link C++ class AliMultiplicityCorrection+;
 
+#pragma link C++ class AliAnalysisTaskdNdetaMC+;
+
 #endif
diff --git a/PWG0/genLevelSimulation/AliAnalysisTaskdNdetaMC.cxx b/PWG0/genLevelSimulation/AliAnalysisTaskdNdetaMC.cxx
new file mode 100644 (file)
index 0000000..7a597f4
--- /dev/null
@@ -0,0 +1,511 @@
+
+//----------------------------------------------------------------
+//      Implementation of Class AliAnalysisTaskdNdetaMC
+//
+// Task used to analize simulations at generation level (i.e. only
+// needs galice.root and Kinematics.root).
+// 
+// The tasks produces multiplicity, eta and pt histograms for
+// different event classes (listed in the enum in the header file).
+// The multiplicity histogram are further produced for different eta
+// ranges, defined in the constructor.
+//
+// Author: Michele Floris, CERN
+//
+//----------------------------------------------------------------
+
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TCanvas.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliVEvent.h"
+#include "AliESDEvent.h"
+#include "AliMCEvent.h"
+#include "AliAODEvent.h"
+
+#include "AliAnalysisTaskdNdetaMC.h"
+#include "TGraphErrors.h"
+#include "AliLog.h"
+#include "AliStack.h"
+#include "AliGenEventHeader.h"
+#include "AliGenPythiaEventHeader.h"
+
+#include <iostream>
+#include "AliPDG.h"
+#include "AliGenDPMjetEventHeader.h"
+#include "TFile.h"
+
+using namespace std;
+
+ClassImp(AliAnalysisTaskdNdetaMC)
+
+Float_t  AliAnalysisTaskdNdetaMC::fEtaMax = 0.5;
+
+AliAnalysisTaskdNdetaMC::AliAnalysisTaskdNdetaMC() 
+  : AliAnalysisTaskSE(), fNchDens(0), fMyOut(0), 
+    fHistIev(0), fHistNParticlesAtMidRapidity(0), fSkipNormalization(0)
+{
+
+  // default constructor
+  for(Int_t ihist = 0; ihist < kNHist; ihist++){
+    fHistEta[ihist]  = 0;
+    fHistPt[ihist]  = 0;
+    for(Int_t ihist2 = 0; ihist2 < kNMultHist; ihist2++){
+      fHistMult[ihist][ihist2] = 0;    
+    }
+  }  
+
+  fEtaBins[kMult05] = 0.5;
+  fEtaBins[kMult10] = 1.0;
+  fEtaBins[kMult14] = 1.3;
+
+  fEtaMax=0.5;
+
+}
+
+//________________________________________________________________________
+AliAnalysisTaskdNdetaMC::AliAnalysisTaskdNdetaMC(const char *name) 
+  : AliAnalysisTaskSE(name), fNchDens(0), fMyOut(0),
+    fHistIev(0), fHistNParticlesAtMidRapidity(0), fSkipNormalization(0)
+{
+  // constructor
+
+  for(Int_t ihist = 0; ihist < kNHist; ihist++){
+    fHistEta[ihist]  = 0;
+    fHistPt[ihist]  = 0;
+    for(Int_t ihist2 = 0; ihist2 < kNMultHist; ihist2++){
+      fHistMult[ihist][ihist2] = 0;    
+    }
+  }
+  fEtaBins[kMult05] = 0.5;
+  fEtaBins[kMult10] = 1.0;
+  fEtaBins[kMult14] = 1.3;
+
+
+  fEtaMax=0.5;
+
+  AliPDG::AddParticlesToPdgDataBase();
+
+  DefineOutput(1, TList::Class());
+}
+
+AliAnalysisTaskdNdetaMC::~AliAnalysisTaskdNdetaMC() {
+
+  // destructor
+
+  if(fMyOut) {
+    // fMyOut owns the histos
+    delete fMyOut;
+    fMyOut = 0;
+  }
+
+}
+
+
+AliAnalysisTaskdNdetaMC::AliAnalysisTaskdNdetaMC(const char *name, const char * fname) 
+  : AliAnalysisTaskSE(name), fNchDens(0), fMyOut(0),
+    fHistIev(0), fHistNParticlesAtMidRapidity(0), fSkipNormalization(0)
+{
+  // This constructor open list from a dndeta file (useful to finalize after merging)
+  for(Int_t ihist = 0; ihist < kNHist; ihist++){
+    fHistEta[ihist]  = 0;
+    fHistPt[ihist]  = 0;
+    for(Int_t ihist2 = 0; ihist2 < kNMultHist; ihist2++){
+      fHistMult[ihist][ihist2] = 0;    
+    }
+  }
+  fEtaBins[kMult05] = 0.5;
+  fEtaBins[kMult10] = 1.0;
+  fEtaBins[kMult14] = 1.3;
+
+
+  fEtaMax=0.5;
+
+  AliPDG::AddParticlesToPdgDataBase();
+
+  TFile * f =new TFile(fname);
+  if (!f) AliFatal(Form("Cannot open file %s!",fname));
+  fMyOut  = (TList*) f->Get("coutput");
+  if (!fMyOut) AliFatal(Form("Cannot get output from file %s!",fname));
+
+  DefineOutput(1, TList::Class());
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskdNdetaMC::UserCreateOutputObjects()
+{
+  // Create histograms
+  // Called once
+
+  fMyOut = new TList();
+
+  fHistEta[kHistINEL] = BookHetaHist("fHistdNdetaMCInel", "MC dN/deta distribution Inel");
+  fHistEta[kHistNSD]  = BookHetaHist("fHistdNdetaMCNSD",  "MC dN/deta distribution NSD");
+  fHistEta[kHistSiD]  = BookHetaHist("fHistdNdetaMCSiD",  "MC dN/deta distribution SiD");
+  fHistEta[kHistND]   = BookHetaHist("fHistdNdetaMCND",   "MC dN/deta distribution Non-Diffractive");
+
+  fHistEta[kHistHL]   = BookHetaHist("fHistdNdetaMCHL",     "MC dN/deta distribution, at least 1 particle |eta| < 1");
+
+
+  fHistPt[kHistINEL] = BookHptHist("fHistdNdptMCInel", "MC dN/dpt distribution Inel (|#eta| < 0.8)");
+  fHistPt[kHistNSD]  = BookHptHist("fHistdNdptMCNSD",  "MC dN/dpt distribution NSD (|#eta| < 0.8)");
+  fHistPt[kHistSiD]  = BookHptHist("fHistdNdptMCSiD",  "MC dN/dpt distribution SiD (|#eta| < 0.8)");
+  fHistPt[kHistND]   = BookHptHist("fHistdNdptMCND",   "MC dN/dpt distribution Non-Diffractive (|#eta| < 0.8)");
+  fHistPt[kHistHL]   = BookHptHist("fHistdNdptMCHL",   "MC dN/dpt distribution at least 1 particle |eta| < 1 (|#eta| < 0.8)");
+
+
+  const char * labelType[] = {"INEL", "NSD", "SiD", "ND", "HL"};
+  for(Int_t ihist = 0; ihist < kNHist; ihist++){ // type
+    for(Int_t ihist2 = 0; ihist2 < kNMultHist; ihist2++) { // eta range
+      fHistMult[ihist][ihist2] = BookMultHisto(Form("fHistMult_%s_%1.1f",labelType[ihist],fEtaBins[ihist2]),
+                                              Form("(dN/dN_{ch})_{|#eta_{max}|<%1.1f} (%s)",fEtaBins[ihist2],labelType[ihist]));    
+    }
+  }
+
+  
+  
+  fNchDens = new TGraphErrors();
+  fNchDens -> SetName  ("fNchDens");
+  fNchDens -> SetTitle ("Charged tracks density at mid-rapidity (|#eta| < fEtaMax)");
+  fNchDens->SetMarkerStyle(kFullCircle);
+  fMyOut->Add(fNchDens);
+
+  fHistNParticlesAtMidRapidity = new TH1I("fHistNParticlesAtMidRapidity","Number of particles at midrapidity", kNHist, -0.5, kNHist-0.5); 
+  fMyOut->Add(fHistNParticlesAtMidRapidity);
+
+  fHistIev = new TH1I("fHistIev","Number of particles at midrapidity", kNHist, -0.5, kNHist-0.5); 
+  fMyOut->Add(fHistIev);
+
+
+  fMyOut->SetOwner();
+
+  // Suppress annoying printout
+  AliLog::SetGlobalLogLevel(AliLog::kError);
+
+}
+
+TH1F* AliAnalysisTaskdNdetaMC::BookHetaHist(const char * name, const char * title) {
+
+  // Book Eta histos
+  TH1F * h = new TH1F(name, title, 200, -10., 10.);
+  h->GetXaxis()->SetTitle("#eta");
+  h->GetYaxis()->SetTitle("dN/d#eta");
+  h->SetMarkerStyle(kFullCircle);
+  h->Sumw2();
+  fMyOut->Add(h);
+  return h;
+
+}
+
+TH1F* AliAnalysisTaskdNdetaMC::BookHptHist(const char * name, const char * title) {
+
+  // Book Pt histos
+  TH1F * h = new TH1F(name, title, 200, 0., 20.);
+  h->GetXaxis()->SetTitle("p_{T}");
+  h->GetYaxis()->SetTitle("dN/dp_{T}");
+  h->SetMarkerStyle(kFullCircle);
+  h->Sumw2();
+  fMyOut->Add(h);
+  return h;
+
+}
+
+TH1F* AliAnalysisTaskdNdetaMC::BookMultHisto(const char * name, const char * title) {
+
+  // Book multiplicity histos
+  Int_t maxmult = 100;
+  TH1F * h = new TH1F(name, title, maxmult, -0.5, maxmult-0.5);
+  h->GetXaxis()->SetTitle("N_{ch}");
+  h->GetYaxis()->SetTitle("dN/dN_{ch}");
+  h->SetMarkerStyle(kFullCircle);
+  h->Sumw2();
+  fMyOut->Add(h);
+  return h;
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskdNdetaMC::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+
+  // also a AliEvent...
+  //  AliVEvent* mcEvent = MCEvent();
+  AliMCEvent* mcEvent = MCEvent();
+  AliGenPythiaEventHeader * headPy  = 0;
+  AliGenDPMjetEventHeader * headPho = 0;
+  AliGenEventHeader * htmp = mcEvent->GenEventHeader();
+  if(!htmp) {
+    AliError("Cannot Get MC Header!!");
+    return;
+  }
+  if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
+    headPy =  (AliGenPythiaEventHeader*) htmp;
+  } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
+    headPho = (AliGenDPMjetEventHeader*) htmp;
+  } else {
+    cout << "Unknown header" << endl;
+    
+  }
+
+  if (!mcEvent) {
+     Printf("ERROR: Could not retrieve MC event");
+     return;
+  }
+
+  //  Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
+  //  Check if the evend is single diffractive
+
+  Bool_t isSD = kFALSE;
+  Bool_t isND = kFALSE;
+  if(headPy)   {
+    //    cout << "Process: " << headPy->ProcessType() << endl;
+    if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
+      isSD = kTRUE; // is single difractive
+    }
+    if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {     
+      isND = kTRUE; // is non-diffractive
+    }
+
+  } else if (headPho) {
+    if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
+      isSD = kTRUE;
+    }       
+    if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6  && headPho->ProcessType() != 7 ) {
+      isND = kTRUE;
+    }       
+  }
+
+  // HL definition: is there at least one particle in |eta|<1?
+  Bool_t isThereOneCentralPart = kFALSE;
+  for (Int_t iTrack = 0; iTrack < mcEvent->GetNumberOfTracks(); iTrack++) {
+    AliMCParticle *track = (AliMCParticle*)mcEvent->GetTrack(iTrack);
+    if (!track) {
+      Printf("ERROR: Could not receive track %d", iTrack);
+      continue;
+    }
+    Bool_t isPrimary = mcEvent->Stack()->IsPhysicalPrimary(iTrack);
+    if (isPrimary && track->Charge() != 0){
+      if (track->Eta() > -1 && track->Eta() < 1) {
+       isThereOneCentralPart = kTRUE;
+       break;
+      }
+    }
+  }
+
+  fHistIev->Fill(kHistINEL);
+  if (!isSD)                 fHistIev->Fill(kHistNSD);
+  if (isSD)                  fHistIev->Fill(kHistSiD);
+  if (isND)                  fHistIev->Fill(kHistND);
+  if (isThereOneCentralPart) fHistIev->Fill(kHistHL);
+  if(!(Int_t(fHistIev->GetBinContent(fHistIev->FindBin(kHistINEL)))%500)) 
+    cout << "Event " << Int_t(fHistIev->GetBinContent(fHistIev->FindBin(kHistINEL))) << endl;
+  
+
+
+  // Track loop
+  Int_t multiplicity[kNHist][kNMultHist] = {{0}};  
+  for (Int_t iTrack = 0; iTrack < mcEvent->GetNumberOfTracks(); iTrack++) {
+    AliMCParticle *track = (AliMCParticle*)mcEvent->GetTrack(iTrack);
+    if (!track) {
+      Printf("ERROR: Could not receive track %d", iTrack);
+      continue;
+    }
+    Bool_t isPrimary = mcEvent->Stack()->IsPhysicalPrimary(iTrack);
+    if (isPrimary && track->Charge() != 0){
+      Bool_t isEtaLess08 = (track->Eta() > -0.8 && track->Eta() < 0.8);
+
+      fHistEta[kHistINEL]->Fill(track->Eta());            
+      if (isEtaLess08) fHistPt[kHistINEL] ->Fill(track->Pt());            
+      if(track->Eta() > -fEtaMax && track->Eta() < fEtaMax) fHistNParticlesAtMidRapidity->Fill(kHistINEL);
+      
+      if(!isSD) {
+       fHistEta[kHistNSD]->Fill(track->Eta());      
+       if (isEtaLess08) fHistPt [kHistNSD]->Fill(track->Pt());      
+       if(track->Eta() > -fEtaMax && track->Eta() < fEtaMax) fHistNParticlesAtMidRapidity->Fill(kHistNSD);
+      }
+      if(isSD) {
+       fHistEta[kHistSiD]->Fill(track->Eta());      
+       if (isEtaLess08) fHistPt [kHistSiD]->Fill(track->Pt());      
+       if(track->Eta() > -fEtaMax && track->Eta() < fEtaMax) fHistNParticlesAtMidRapidity->Fill(kHistSiD);
+      }
+      if (isND) {
+       fHistEta[kHistND]->Fill(track->Eta());      
+       if (isEtaLess08) fHistPt [kHistND]->Fill(track->Pt());      
+       if(track->Eta() > -fEtaMax && track->Eta() < fEtaMax) fHistNParticlesAtMidRapidity->Fill(kHistND);
+      }
+      if (isThereOneCentralPart) {
+       fHistEta[kHistHL]->Fill(track->Eta());      
+       if (isEtaLess08) fHistPt [kHistHL]->Fill(track->Pt());      
+       if(track->Eta() > -fEtaMax && track->Eta() < fEtaMax) fHistNParticlesAtMidRapidity->Fill(kHistHL);
+      }
+
+      // fill array of multiplicity for different classes of events
+      // and in different eta ranges
+      for(Int_t ihist = 0; ihist < kNMultHist; ihist++){
+       if(track->Eta() > -fEtaBins[ihist] && track->Eta() < fEtaBins[ihist]) {
+         multiplicity[kHistINEL][ihist]++;
+         if(!isSD)                  multiplicity[kHistNSD][ihist]++;
+         if(isSD)                   multiplicity[kHistSiD][ihist]++;
+         if(isND)                   multiplicity[kHistND] [ihist]++;
+         if(isThereOneCentralPart)  multiplicity[kHistHL] [ihist]++;
+  
+       }
+      }
+    }
+  } //track loop 
+
+
+
+  // Fill multiplicity histos
+  for(Int_t ihist = 0; ihist < kNMultHist; ihist++){
+    fHistMult[kHistINEL][ihist] ->Fill(multiplicity[kHistINEL][ihist]);
+    if(!isSD)                  fHistMult[kHistNSD][ihist]->Fill(multiplicity[kHistNSD][ihist]);
+    if(isSD)                   fHistMult[kHistSiD][ihist]->Fill(multiplicity[kHistSiD][ihist]);
+    if(isND)                   fHistMult[kHistND] [ihist]->Fill(multiplicity[kHistND] [ihist]);
+    if(isThereOneCentralPart)  fHistMult[kHistHL] [ihist]->Fill(multiplicity[kHistHL] [ihist]);    
+  }
+
+  // Post output data.
+  PostData(1, fMyOut);
+}      
+
+//________________________________________________________________________
+void AliAnalysisTaskdNdetaMC::Terminate(Option_t *) 
+{
+  // Draw result to the screen
+  // Called once at the end of the query   
+
+  //  fHistPt = dynamic_cast<TH1F*> (GetOutputData(1));
+  fMyOut  = dynamic_cast<TList*> (GetOutputData(1));
+  
+  Finalize();
+
+}
+
+
+
+void AliAnalysisTaskdNdetaMC::Finalize() {
+
+  // Scale histos and computes dNdeta
+
+  fHistEta[kHistINEL] = (TH1F*) fMyOut->FindObject("fHistdNdetaMCInel");
+  fHistEta[kHistNSD]  = (TH1F*) fMyOut->FindObject("fHistdNdetaMCNSD" );
+  fHistEta[kHistSiD]  = (TH1F*) fMyOut->FindObject("fHistdNdetaMCSiD" );
+  fHistEta[kHistND]   = (TH1F*) fMyOut->FindObject("fHistdNdetaMCND"  );
+  fHistEta[kHistHL]   = (TH1F*) fMyOut->FindObject("fHistdNdetaMCHL"  );
+  fHistPt[kHistINEL] = (TH1F*) fMyOut->FindObject("fHistdNdptMCInel");
+  fHistPt[kHistNSD]  = (TH1F*) fMyOut->FindObject("fHistdNdptMCNSD" );
+  fHistPt[kHistSiD]  = (TH1F*) fMyOut->FindObject("fHistdNdptMCSiD" );
+  fHistPt[kHistND]   = (TH1F*) fMyOut->FindObject("fHistdNdptMCND"  );
+  fHistPt[kHistHL]   = (TH1F*) fMyOut->FindObject("fHistdNdptMCHL"  );
+  fNchDens            = (TGraphErrors*) fMyOut->FindObject("fNchDens");
+
+  const char * labelType[] = {"INEL", "NSD", "SiD", "ND", "HL"};
+  for(Int_t ihist = 0; ihist < kNHist; ihist++){
+    for(Int_t ihist2 = 0; ihist2 < kNMultHist; ihist2++){
+      fHistMult[ihist][ihist2] = (TH1F*) fMyOut->FindObject(Form("fHistMult_%s_%1.1f",labelType[ihist],fEtaBins[ihist2]));      
+      if (!fHistMult[ihist][ihist2]) cout << "Cannot get histo " << Form("fHistMult_%s_%1.1f",labelType[ihist],fEtaBins[ihist2]) << endl;
+    }
+  }
+
+  fHistIev  = (TH1I*) fMyOut->FindObject("fHistIev"  );
+  fHistNParticlesAtMidRapidity  = (TH1I*) fMyOut->FindObject("fHistNParticlesAtMidRapidity"  );
+    
+
+//   fHistIev->Draw();
+
+  cout << "Eta Max: " << fEtaMax << endl;
+
+
+  for(Int_t ihist = 0; ihist < kNHist; ihist++){
+
+    Int_t iev = (Int_t) fHistIev->GetBinContent(fHistIev->FindBin(ihist));
+    Int_t npm = (Int_t) fHistNParticlesAtMidRapidity->GetBinContent(fHistNParticlesAtMidRapidity->FindBin(ihist));
+
+    // Renormalize dNdeta to the number of events
+    //    cout << fHistEta[ihist] << " " << iev << endl;
+    
+    // compute density at midrapidity (Delta_y = 1);
+    Int_t firstbin =  fHistEta[ihist]->FindBin(-fEtaMax);
+    Int_t lastbin =  fHistEta[ihist]->FindBin(fEtaMax);
+
+    if (!fSkipNormalization) {
+      fHistEta[ihist]->Scale(1./iev, "width");       
+      fHistPt[ihist] ->Scale(1./iev, "width");       
+      for(Int_t ihist2 = 0; ihist2 < kNMultHist; ihist2++){
+       fHistMult[ihist][ihist2]->Scale(1./iev);
+      }
+    }
+
+
+    Double_t meaneta = 0;
+    Double_t meanerr = 0;
+    Double_t sumweight  = 0;
+    for(Int_t ibin = firstbin; ibin <=lastbin ; ibin++){
+      Double_t x    = fHistEta[ihist]->GetBinContent(ibin);
+      Double_t xerr = fHistEta[ihist]->GetBinError(ibin);
+      
+      
+      Double_t xerr2 = xerr*xerr;
+      if(xerr2){
+       //      cout << "xe2 " << xerr2 << endl;
+       Double_t weight = 1. / xerr2;
+       sumweight += weight;
+       meaneta += weight * x;
+      }
+      
+    }
+    
+    if(sumweight){
+      meaneta /= sumweight;
+      meanerr = TMath::Sqrt(1./ sumweight);
+    }
+    else {
+      meaneta = 0;
+      meanerr = 0;
+    }
+
+
+    Double_t range = 2*fEtaMax;
+
+//     meaneta /= fHistEta[ihist]->GetBinWidth(firstbin);
+//     meanerr /= fHistEta[ihist]->GetBinWidth(firstbin);
+    cout << "Histo: " << fHistEta[ihist]->GetName() << endl;
+    cout << " Evts: " << iev << endl;
+    
+    cout << " Density at midrapidity:       " << meaneta << "+-" << meanerr << endl;
+
+    // Direct computation
+    Double_t errdir = TMath::Sqrt(npm) / iev;
+    Double_t etadir = Double_t(npm) / iev;
+    
+    cout << " Density at midrapidity: (dir) " << etadir/range << "+-" << errdir/range << endl;
+  
+    fNchDens->SetPoint     (ihist, ihist+1, etadir);
+    fNchDens->SetPointError(ihist, 0, errdir);
+
+    cout << " Density at midrapidity: (TH1, Eta<0.5) " << fHistMult[ihist][0]->GetMean() << "+-" << fHistMult[ihist][0]->GetMeanError() << endl;
+
+  }
+  
+  // Draw fHistEtaAll
+  TCanvas *c1 = new TCanvas("AliAnalysisTaskdNdetaMC","dNdetaMC",10,10,800,400);
+  //  c1->cd(1)->SetLogy();
+  c1->Divide(2,1);
+  c1->cd(1);
+  fHistEta[0]->DrawCopy("");
+  fHistEta[1]->SetLineColor(kRed);
+  fHistEta[1]->SetMarkerColor(kRed);
+  fHistEta[1]->SetMarkerStyle(kOpenSquare);
+  fHistEta[1]->DrawCopy("same");
+  c1->cd(2);
+  fNchDens->Draw("AP");
+
+}
diff --git a/PWG0/genLevelSimulation/AliAnalysisTaskdNdetaMC.h b/PWG0/genLevelSimulation/AliAnalysisTaskdNdetaMC.h
new file mode 100644 (file)
index 0000000..43caa80
--- /dev/null
@@ -0,0 +1,54 @@
+#ifndef AliAnalysisTaskdNdetaMC_H
+#define AliAnalysisTaskdNdetaMC_H
+
+//
+// Task used to analize simulations at generation level (i.e. only
+// needs galice.root and Kinematics.root).
+// 
+
+class TH1F;
+class TH1I;
+class TGraphErrors;
+
+enum {kHistINEL,kHistNSD,kHistND,kHistSiD,kHistHL,kNHist};
+enum {kMult05,kMult10,kMult14,kNMultHist};// 
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskdNdetaMC : public AliAnalysisTaskSE {
+ public:
+  AliAnalysisTaskdNdetaMC();
+  AliAnalysisTaskdNdetaMC(const char *name );
+  AliAnalysisTaskdNdetaMC(const char *name, const char *file );
+  virtual ~AliAnalysisTaskdNdetaMC();
+  
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+  TH1F* BookHetaHist(const char * name, const char * title);
+  TH1F* BookHptHist(const char * name, const char * title);
+  TH1F* BookMultHisto(const char * name, const char * title) ;
+  void SetEtaMax(Float_t eta) {fEtaMax = eta;}
+
+  void SkipNormalization(Bool_t flag = kTRUE) { fSkipNormalization = flag; }
+  void Finalize();
+  TList * GetList() const { return fMyOut;} 
+ private:
+  TH1F         *fHistEta[kNHist]; //Eta spectrum 
+  TH1F         *fHistPt[kNHist]; //Eta spectrum  , |eta| < 0.8
+  TGraphErrors *fNchDens; // <dN/deta>
+  TList * fMyOut; // list of output histos
+  TH1I * fHistIev; // number of events per class
+  TH1I * fHistNParticlesAtMidRapidity;  // number of particles at midrapidity per class
+  static Float_t fEtaMax; // max eta
+  Bool_t fSkipNormalization; // Use this when you are running the job on the grid, so that you can normalize dNdeta after merging
+
+  Float_t  fEtaBins[kNMultHist];    // array of eta_max values
+  TH1F * fHistMult[kNHist][kNMultHist];   // array of multiplicity histos in the different eta ranges values, for the different event classes
+
+  AliAnalysisTaskdNdetaMC(const AliAnalysisTaskdNdetaMC&); // not implemented
+  AliAnalysisTaskdNdetaMC& operator=(const AliAnalysisTaskdNdetaMC&); // not implemented
+  
+  ClassDef(AliAnalysisTaskdNdetaMC, 2); 
+};
+
+#endif
diff --git a/PWG0/genLevelSimulation/MergeCollectionFromGrid.C b/PWG0/genLevelSimulation/MergeCollectionFromGrid.C
new file mode 100644 (file)
index 0000000..89f92f5
--- /dev/null
@@ -0,0 +1,58 @@
+#include <iostream>
+#include "TSystem.h"
+#include "TFileMerger.h"
+#include "TGrid.h"
+#include "TAlienCollection.h"
+#include "TFile.h"
+#include "TH1F.h"
+#include "TList.h"
+#include "TH1I.h"
+#include "TMath.h"
+#include "TGraphErrors.h"
+#include "AliAnalysisTaskdNdetaMC.h"
+#include "TROOT.h"
+
+using namespace std;
+
+
+
+void MergeCollectionFromGrid(const char * incollection = "test.xml", const char * outputfile= "dndeta_merged.root")
+{
+  // for running with root only
+  gSystem->Load("libTree.so");
+  gSystem->Load("libGeom.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libSTEERBase.so");
+  gSystem->Load("libESD.so");
+  gSystem->Load("libAOD.so"); 
+
+  // load analysis framework
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+
+  // load dndeta task
+  gROOT->LoadMacro("AliAnalysisTaskdNdetaMC.cxx+");
+
+
+  TFileMerger * fileMerger  = new TFileMerger(0); // dont merge local files
+
+  TGrid::Connect("alien://");
+  TGridCollection * coll = TAlienCollection::Open (incollection);
+  Int_t  ifile=0;
+  while(coll->Next()){
+    fileMerger->AddFile(TString("alien://")+coll->GetLFN());
+    ifile++;
+    //    if(ifile>2) break;
+  }
+  fileMerger->OutputFile("tmp.root");
+  fileMerger->Merge();
+
+  // Reopen the merged file, normalize histos and save them back.
+  // SOME DUPLICATED CODE... SOME CLEAN UP WOULD BE GOOD
+  AliAnalysisTaskdNdetaMC * localTask  = new AliAnalysisTaskdNdetaMC("merger", "tmp.root");
+  localTask->Finalize();
+  
+  TFile * fout = new TFile (outputfile, "recreate");
+  localTask->GetList()->Write();
+  fout->Close();
+}
diff --git a/PWG0/genLevelSimulation/fastGen.C b/PWG0/genLevelSimulation/fastGen.C
new file mode 100644 (file)
index 0000000..d820e38
--- /dev/null
@@ -0,0 +1,164 @@
+typedef enum {kPhojet = -1,kPyTuneCDFA=100,kPyTuneAtlasCSC=306, kPyTuneCMS6D6T=109, kPyTunePerugia0=320 } Tune_t;
+
+AliGenerator*  CreateGenerator(Tune_t tune = kPyTuneCDFA , Float_t energy);
+
+void fastGen(Tune_t tune = kPyTuneCDFA , Float_t energy, Int_t nev = 1, TString process)
+{
+  // Add all particles to the PDG database
+  AliPDG::AddParticlesToPdgDataBase();
+
+  // set the random seed
+  TDatime date;
+  UInt_t seed    = date.Get()+gSystem->GetPid();
+  gRandom->SetSeed(seed);
+  cout<<"Seed for random number generation= "<<seed<<endl; 
+
+
+  //  Runloader  
+  AliRunLoader* rl = AliRunLoader::Open("galice.root", "FASTRUN","recreate");
+    
+  rl->SetCompressionLevel(2);
+  rl->SetNumberOfEventsPerFile(nev);
+  rl->LoadKinematics("RECREATE");
+  rl->MakeTree("E");
+  gAlice->SetRunLoader(rl);
+
+  //  Create stack
+  rl->MakeStack();
+  AliStack* stack      = rl->Stack();
+  //  Header
+  AliHeader* header = rl->GetHeader();
+  //
+  //  Create and Initialize Generator
+  AliGenerator *gener = CreateGenerator(tune,energy);
+  gener->Init();
+  // if nsd switch off single diffraction
+  if ( process == "NSD"){
+    if(tune != kPhojet) {
+      AliPythia::Instance()->  SetMSUB(92,0);             // single diffraction AB-->XB
+      AliPythia::Instance()-> SetMSUB(93,0);             // single diffraction AB-->AX
+    }
+    else {
+      cout << "NSD not yet implemented in the phojet case" << endl;
+      exit(1);
+    }
+  }
+  gener->SetStack(stack);
+    
+  //
+  //                        Event Loop
+  //
+  Int_t iev;
+     
+  for (iev = 0; iev < nev; iev++) {
+
+    if(!(iev%500)) printf("\n \n Event number %d \n \n", iev);
+       
+    //  Initialize event
+    header->Reset(0,iev);
+    rl->SetEventNumber(iev);
+    stack->Reset();
+    rl->MakeTree("K");
+    // stack->ConnectTree();
+    
+    //  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 mpart  = MPart->GetPdgCode();
+    //             printf("Particle %d\n", mpart);
+    //         }
+       
+    //  Finish event
+    header->SetNprimary(stack->GetNprimary());
+    header->SetNtrack(stack->GetNtrack());  
+    //      I/O
+    // 
+    stack->FinishEvent();
+    header->SetStack(stack);
+    rl->TreeE()->Fill();
+    rl->WriteKinematics("OVERWRITE");
+
+  } // event loop
+    //
+    //                         Termination
+    //  Generator
+  gener->FinishRun();
+  //  Write file
+  rl->WriteHeader("OVERWRITE");
+  gener->Write();
+  rl->Write();
+    
+}
+
+
+AliGenerator*  CreateGenerator(Tune_t tune, Float_t energy)
+{
+
+  
+
+  if (tune == -1) {
+
+    // phojet
+    AliGenDPMjet* gener = new AliGenDPMjet(1);
+    gener->SetProcess(kDpmMb);
+    gener->SetProjectile("P", 1, 1);
+    gener->SetTarget("P", 1, 1);
+
+    gener->SetEnergyCMS(energy);
+    return gener;
+  }
+
+  if (tune != kPyTuneAtlasCSC && tune != kPyTuneCDFA && tune != kPyTuneCMS6D6T && tune != kPyTunePerugia0) {
+    
+    Printf("Unknown pythia tune, quitting");
+    exit(1);
+
+  } else {
+
+    AliGenPythia * gener =  new AliGenPythia(1);
+    //
+    //  set pythia tune
+    gener->SetTune(tune);
+
+    //   structure function  
+    if(tune == kPyTuneAtlasCSC) {
+      cout << "Setting structure function" << endl;      
+      gener->SetStrucFunc(kCTEQ61);
+    }
+    if(tune == kPyTuneCMS6D6T) {
+      cout << "Setting structure function" << endl;      
+      gener->SetStrucFunc(kCTEQ6l);
+    }
+    if(tune == kPyTunePerugia0) {
+      cout << "Setting new parton shower" << endl;      
+      gener->UseNewMultipleInteractionsScenario();
+    }
+    //   charm, beauty, charm_unforced, beauty_unforced, jpsi, jpsi_chi, mb
+    gener->SetProcess(kPyMb);
+    //   Centre of mass energy 
+    gener->SetEnergyCMS(energy);
+
+    // Set target/projectile // Is this needded?
+    gener->SetProjectile("P", 1, 1);
+    gener->SetTarget("P", 1, 1);
+
+    return gener;
+  }
+}
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/PWG0/genLevelSimulation/runTask.C b/PWG0/genLevelSimulation/runTask.C
new file mode 100644 (file)
index 0000000..76bedcd
--- /dev/null
@@ -0,0 +1,77 @@
+void runTask(Float_t etamax=0.5,const char * incollection = 0, const char * outfile = "dndeta.root", Bool_t skipNorm = kFALSE)
+{
+  // for running with root only
+  gSystem->Load("libTree.so");
+  gSystem->Load("libGeom.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libSTEERBase.so");
+  gSystem->Load("libESD.so");
+  gSystem->Load("libAOD.so"); 
+
+  // load analysis framework
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+
+
+  TChain * chain = new TChain ("TE");
+  if (incollection == 0) {
+    chain->Add("galice.root");
+  }
+  else if (TString(incollection).Contains("xml")){
+    TGrid::Connect("alien://");
+    TAlienCollection * coll = TAlienCollection::Open (incollection);
+    while(coll->Next()){
+      chain->Add(TString("alien://")+coll->GetLFN());
+    }
+  } else {
+    ifstream file_collect(incollection);
+    TString line;
+    while (line.ReadLine(file_collect) ) {
+      chain->Add(line.Data());
+    }
+  }
+  chain->GetListOfFiles()->Print();
+
+  // for includes use either global setting in $HOME/.rootrc
+  // ACLiC.IncludePaths: -I$(ALICE_ROOT)/include
+  // or in each macro
+  gSystem->AddIncludePath("-I$ALICE_ROOT/include");
+
+  // Create the analysis manager
+  AliAnalysisManager *mgr = new AliAnalysisManager("dNdeta");
+
+  AliVEventHandler* esdH = new AliESDInputHandler;
+  mgr->SetInputEventHandler(esdH);
+
+  // Create tasks
+  gROOT->LoadMacro("AliAnalysisTaskdNdetaMC.cxx++g");
+
+  AliAnalysisTask *task1 = new AliAnalysisTaskdNdetaMC("TaskdNdeta");
+  ((AliAnalysisTaskdNdetaMC*)task1)->SetEtaMax(etamax);
+  if (skipNorm) ((AliAnalysisTaskdNdetaMC*)task1)->SkipNormalization();
+  // Enable MC event handler
+
+  AliMCEventHandler* handler = new AliMCEventHandler;
+  handler->SetReadTR(kFALSE);
+  mgr->SetMCtruthEventHandler(handler);
+
+  // Add tasks
+  mgr->AddTask(task1);
+
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("coutput", TList::Class(),    AliAnalysisManager::kOutputContainer, outfile);
+
+  // Connect input/output
+  mgr->ConnectInput(task1, 0, cinput);
+  mgr->ConnectOutput(task1, 1, coutput1);
+
+  // Enable debug printouts
+  mgr->SetDebugLevel(0);
+
+  if (!mgr->InitAnalysis()) return;
+
+  mgr->PrintStatus();
+
+  mgr->StartAnalysis("local", chain);
+}
diff --git a/PWG0/genLevelSimulation/runTaskProof.C b/PWG0/genLevelSimulation/runTaskProof.C
new file mode 100644 (file)
index 0000000..d58e338
--- /dev/null
@@ -0,0 +1,57 @@
+void runTaskProof(const char * dataset, const char * datasetpath="/COMMON/COMMON/", const char * outdir="NchDistributions") {
+  gEnv->SetValue("XSec.GSI.DelegProxy","2");
+  TProof::Open("alicecaf");
+  
+  //  gSystem->AddIncludePath("-I${ALICE_ROOT}/include/ -I${ALICE_ROOT}/PWG0/ -I${ALICE_ROOT}/PWG0/dNdEta/");
+  gSystem->AddIncludePath("-I${ALICE_ROOT}/include/");
+  gProof->UploadPackage("$ALICE_ROOT/STEERBase");
+  gProof->EnablePackage("$ALICE_ROOT/STEERBase");
+  gProof->UploadPackage("$ALICE_ROOT/ESD");
+  gProof->EnablePackage("$ALICE_ROOT/ESD");
+  gProof->UploadPackage("$ALICE_ROOT/AOD");
+  gProof->EnablePackage("$ALICE_ROOT/AOD");
+  gProof->UploadPackage("$ALICE_ROOT/ANALYSIS");
+  gProof->EnablePackage("$ALICE_ROOT/ANALYSIS");
+  gProof->UploadPackage("$ALICE_ROOT/ANALYSISalice");
+  gProof->EnablePackage("$ALICE_ROOT/ANALYSISalice");
+
+    // Make the analysis manager
+  AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
+  //AliVEventHandler *esdH = new AliESDInputHandler();
+  AliESDInputHandler* esdH = new AliESDInputHandler();
+  //esdH->SetActiveBranches("ESDfriend");
+  mgr->SetInputEventHandler(esdH);
+       
+  AliMCEventHandler *mc = new AliMCEventHandler();
+  mc->SetReadTR(kFALSE);
+  mgr->SetMCtruthEventHandler(mc);
+
+  // assign simple task
+  gProof->Load("AliAnalysisTaskdNdetaMC.cxx+g");
+  AliAnalysisTask *task = new AliAnalysisTaskdNdetaMC("TaskdNdeta");
+  mgr->AddTask(task);
+
+  char outfname[1000];
+  sprintf(outfname, "dndeta_%s.root",dataset);
+  
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
+  
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("clist1",TList::Class(),
+                                                           AliAnalysisManager::kOutputContainer,outfname);
+
+  mgr->ConnectInput(task,0,cinput1);
+  mgr->ConnectOutput(task,1,coutput1);
+       
+  if (!mgr->InitAnalysis()) return;
+       
+  mgr->PrintStatus();
+  mgr->StartAnalysis("proof",Form("%s%s#esdTree",datasetpath,dataset));
+  //  mgr->StartAnalysis("proof","/COMMON/COMMON/LHC09b14_7TeV_0.5T_Phojet#esdTree");
+
+  if (!strcmp(outdir,"")){
+    gSystem->Exec(Form("mv %s %s", outfname, outdir));
+  }
+
+}
diff --git a/PWG0/genLevelSimulation/rungen.C b/PWG0/genLevelSimulation/rungen.C
new file mode 100644 (file)
index 0000000..a7a250a
--- /dev/null
@@ -0,0 +1,34 @@
+typedef enum {kPhojet = -1, kPyTuneCDFA=100,kPyTuneAtlasCSC=306, kPyTuneCMS6D6T=109, kPyTunePerugia0=320 } Tune_t;
+
+void rungen(Tune_t tune = kPyTuneCDFA, Float_t energy, Int_t nev=1, TString process, Int_t random_index = 0){
+  // Simulation and reconstruction
+  TStopwatch timer;
+  timer.Start();
+  gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I$ALICE_ROOT");
+  gSystem->Load("liblhapdf.so");      // Parton density functions
+  if (tune == kPhojet) {
+    cout << "Loading phojet" << endl;
+    
+    // => phojet
+    gSystem->Load("libpythia6.so");     // Pythia
+    gSystem->Load("libdpmjet.so");     // 
+    gSystem->Load("libTDPMjet.so");     // 
+  } 
+  else if (tune == kPyTunePerugia0) {
+    gSystem->Load("libEGPythia6.so");   // TGenerator interface 
+    gSystem->Load("libpythia6.4.21.so");     // Pythia
+    gSystem->Load("libAliPythia6.so");  // ALICE specific implementations
+  }
+
+  else {
+    gSystem->Load("libEGPythia6.so");   // TGenerator interface 
+    gSystem->Load("libqpythia.so");     // Pythia
+    gSystem->Load("libAliPythia6.so");  // ALICE specific implementations
+   }
+  gROOT->LoadMacro("fastGen.C");
+  fastGen(tune, energy, nev, process);
+
+
+  timer.Stop();
+  timer.Print();
+}
index 2570a3d..5663683 100644 (file)
@@ -10,13 +10,14 @@ SRCS = dNdEta/dNdEtaAnalysis.cxx \
       AliCorrection.cxx \
       AliPWG0Helper.cxx \
       AliUnfolding.cxx \
-      multiplicity/AliMultiplicityCorrection.cxx
+      multiplicity/AliMultiplicityCorrection.cxx \
+      genLevelSimulation/AliAnalysisTaskdNdetaMC.cxx
 
 HDRS = $(SRCS:.cxx=.h)
 
 DHDR= PWG0baseLinkDef.h
 
-EINCLUDE= PWG0/dNdEta PWG0/multiplicity
+EINCLUDE= PWG0/dNdEta PWG0/multiplicity PWG0/genLevelSimulation
 
 ifeq (win32gcc,$(ALICE_TARGET))
 PACKSOFLAGS:= $(SOFLAGS) -L$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET) -lSTEERBase \
@@ -24,4 +25,4 @@ PACKSOFLAGS:= $(SOFLAGS) -L$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET) -lSTEERBase \
                          -L$(ROOTLIBDIR) -lEG
 endif
 
-EXPORT=AliUnfolding.h
\ No newline at end of file
+EXPORT=AliUnfolding.h