PWGLF/FORWARD/analysis2/sim/FastSim.C

   1 #ifndef FASTSIM_H
   2 #define FASTSIM_H
   3 #include <TSelector.h>
   4 #include <TQObject.h>
   5 #ifndef __CINT__
   6 # include "AliGenerator.h"
   7 # include "AliRunLoader.h"
   8 # include "AliStack.h"
   9 # include "AliHeader.h"
  10 # include "AliGenEventHeader.h"
  11 # include "AliRun.h"
  12 # include "AliCollisionGeometry.h"
  13 # include "AliGenPythiaEventHeader.h"
  14 # include "AliGenDPMjetEventHeader.h"
  15 # include "AliGenGeVSimEventHeader.h"
  16 # include "AliGenHerwigEventHeader.h"
  17 # include <TROOT.h>
  18 # include <TString.h>
  19 # include <TMath.h>
  20 # include <TParticle.h>
  21 # include <TH1.h>
  22 # include <TTree.h>
  23 # include <TClonesArray.h>
  24 # include <TList.h>
  25 # include <TProof.h>
  26 # include <TParticlePDG.h>
  27 # include <TStopwatch.h>
  28 # include <TFile.h>
  29 # include <TProofOutputFile.h>
  30 # include <TCanvas.h>
  31 # include <TTimer.h>
  32 # include <fstream>
  33 #else
  34 class AliGenerator;
  35 class AliRunLoader;
  36 class AliStack;
  37 class AliHeader;
  38 class AliGenEventHeader;
  39 class TH1;
  40 class TTree;
  41 class TClonesArray;
  42 class TBrowser;
  43 class TList;
  44 class TFile;
  45 class TProofOutputFile;
  46 class TCanvas;
  47 class TVirtualPad;
  48 class TTimer;
  49 #endif
  50
  51 //====================================================================
  52 /**
  53  * Monitor output objects
  54  */
  55 struct FastMonitor : public TObject, public TQObject
  56 {
  57   /**
  58    * Constructor
  59    *
  60    *
  61    * @return
  62    */
  63   FastMonitor(TSelector* s=0)
  64     : fName("FastMonitor"),
  65       fCanvas(0),
  66       fSelector(s)
  67   {
  68     if (gROOT->IsBatch()) {
  69       Warning("FastMonitor", "Batch processing, no monitoring");
  70       return;
  71     }
  72
  73     if (gProof) {
  74       fName = gProof->GetSessionTag();
  75       gDirectory->Add(this);
  76       Bool_t ret = gProof->Connect("Feedback(TList *objs)", "FastMonitor", this,
  77                                    "Feedback(TList *objs)");
  78       if (!ret) {
  79         Warning("FastMonitor", "Failed to connect to Proof");
  80         return;
  81       }
  82     }
  83     else if (!s) return;
  84
  85     fCanvas = new TCanvas(fName, Form("Monitor %s", fName.Data()), 1000, 800);
  86     fCanvas->SetFillColor(0);
  87     fCanvas->SetFillStyle(0);
  88     fCanvas->SetTopMargin(0.01);
  89     fCanvas->SetRightMargin(0.01);
  90
  91     fCanvas->Divide(2,2);
  92     RegisterDraw(1, "type",   "", 0);
  93     RegisterDraw(2, "b",      "", 0);
  94     RegisterDraw(3, "cent",   "", 0);
  95     RegisterDraw(4, "dNdeta", "", 0x8);
  96   }
  97   /**
  98    * Register a draw of a an object
  99    *
 100    * @param i      Pad number
 101    * @param name   Name of object
 102    * @param option Drawing option
 103    * @param flags  Flags
 104    *
 105    *  - 0x1   Log(x)
 106    *  - 0x2   Log(y)
 107    *  - 0x4   Log(z)
 108    *  - 0x8   Scale to events and bin width
 109    */
 110   void RegisterDraw(Int_t i,
 111                     const char* name,
 112                     const char* option,
 113                     UShort_t    flags=0)
 114   {
 115     TVirtualPad* p = fCanvas->GetPad(i);
 116     if (!p) {
 117       Warning("RegisterDraw", "Not enough sub-pads (%d)", i);
 118       return;
 119     }
 120     p->SetFillColor(0);
 121     p->SetFillStyle(0);
 122     p->SetTopMargin(0.01);
 123     p->SetRightMargin(0.01);
 124     p->SetName(Form("p_%s", name));
 125     p->SetTitle(option);
 126     if (flags & 0x1) p->SetLogx();
 127     if (flags & 0x2) p->SetLogy();
 128     if (flags & 0x4) p->SetLogz();
 129     if (flags & 0x8) p->SetBit(BIT(15));
 130   }
 131   /**
 132    * Desctructor
 133    */
 134   virtual ~FastMonitor()
 135   {
 136     if (!gProof) return;
 137     gProof->Disconnect("Feedback(TList *objs)",this,
 138                        "Feedback(TList* objs)");
 139   }
 140   /**
 141    * Set name of this object
 142    *
 143    * @param name Name
 144    */
 145   void SetName(const char* name) { fName = name; }
 146   /**
 147    * Get the name of this object
 148    *
 149    * @return Name
 150    */
 151   const char* GetName() const { return fName.Data(); }
 152   /**
 153    * Find pad corresponding to an object
 154    *
 155    * @param name Name of object
 156    *
 157    * @return Pointer to pad or null
 158    */
 159   TVirtualPad* FindPad(const TString& name)
 160   {
 161     TVirtualPad* p = 0;
 162     Int_t        i = 1;
 163     TString      t = Form("p_%s", name.Data());
 164     while ((p = fCanvas->GetPad(i))) {
 165       if (t.EqualTo(p->GetName())) return p;
 166       i++;
 167     }
 168     return 0;
 169   }
 170   /**
 171    * Called when we get notified of
 172    *
 173    * @param objs List of monitored objects
 174    */
 175   void Feedback(TList* objs)
 176   {
 177     // Info("FeedBack", "List is %p", objs);
 178     // if (objs) objs->ls();
 179     if (!fCanvas) return;
 180
 181     TList* l = static_cast<TList*>(objs->FindObject("histograms"));
 182     if (!l) {
 183       Warning("Feedback", "No histograms");
 184       return;
 185     }
 186     Int_t nEvents = 1;
 187     TObject* oIpz = l->FindObject("ipZ");
 188     if (oIpz && oIpz->IsA()->InheritsFrom(TH1::Class()))
 189       nEvents = static_cast<TH1*>(oIpz)->GetEntries();
 190     else
 191       Warning("Feedback", "Histogram ipZ not found");
 192
 193     TIter next(l);
 194     TObject* o = 0;
 195     while ((o = next())) {
 196       TVirtualPad* p = FindPad(o->GetName());
 197       if (!p)
 198         // Info("FeedBack", "no pad for %s", o->GetName());
 199         continue;
 200
 201       p->cd();
 202       if (o->IsA()->InheritsFrom(TH1::Class())) {
 203         TH1* h = static_cast<TH1*>(o);
 204         TH1* c = h->DrawCopy(p->GetTitle());
 205         c->SetDirectory(0);
 206         c->SetBit(TObject::kCanDelete);
 207         if (p->TestBit(BIT(15))) {
 208           Info("Feedback", "Scaling %s by 1./%d and width",
 209                c->GetName(), nEvents);
 210           c->Scale(1./nEvents, "width");
 211         }
 212       }
 213       else {
 214         TObject* c = o->DrawClone(p->GetTitle());
 215         c->SetBit(TObject::kCanDelete);
 216       }
 217       p->Modified();
 218     }
 219     fCanvas->Modified();
 220     fCanvas->Update();
 221     fCanvas->cd();
 222   }
 223   /**
 224    * Function to handle connect signals
 225    *
 226    */
 227   void Handle()
 228   {
 229     HandleTimer(0);
 230   }
 231   /**
 232    * Function to handle timer events
 233    */
 234   Bool_t HandleTimer(TTimer*)
 235   {
 236     Info("HandleTimer", "Selector=%p", fSelector);
 237     if (!fSelector) return false;
 238     Feedback(fSelector->GetOutputList());
 239     return true;
 240   }
 241   /** Our name */
 242   TString fName;
 243   /** Our canvas */
 244   TCanvas* fCanvas;
 245   /** Possibly link to selector */
 246   TSelector* fSelector;
 247   ClassDef(FastMonitor,1);
 248 };
 249
 250
 251 //====================================================================
 252 /**
 253  * Run a event generator simulation
 254  */
 255 struct FastSim : public TSelector
 256 {
 257   /**
 258    * Constructor
 259    *
 260    * @param eg     Event generator
 261    * @param runNo  Run number to simulate
 262    * @param bMin   Lease impact parameter
 263    * @param bMax   Largest impact parameter
 264    */
 265   FastSim(const char* eg="",
 266           ULong_t runNo=0,
 267           Double_t bMin=0,
 268           Double_t bMax=20,
 269           Long64_t nEvents=0)
 270     : TSelector(),
 271       fEGName(eg),
 272       fRunNo(runNo),
 273       fBMin(bMin),
 274       fBMax(bMax),
 275       fGRP(0),
 276       fNEvents(nEvents),
 277       fGenerator(0),
 278       fRunLoader(0),
 279       fStack(0),
 280       fHeader(0),
 281       fTree(0),
 282       fParticles(0),
 283       fList(0),
 284       fHEta(0),
 285       fHIpz(0),
 286       fHType(0),
 287       fHCent(0),
 288       fHB(0),
 289       fHPhiR(0),
 290       fHTime(0),
 291       fProofFile(0),
 292       fFile(0),
 293       fFileName("")
 294   {}
 295   const char* FileName() const
 296   {
 297     static TString fn;
 298     if (fn.IsNull()) {
 299       if (!fFileName.IsNull())  fn = fFileName;
 300       else {
 301         const char* egName = (fGenerator ?
 302                               fGenerator->GetName() :
 303                               fEGName.Data());
 304         fn = Form("%s_%09d", egName, fRunNo);
 305         if (fNEvents > 0) {
 306           if (fNEvents >= 1000000)
 307             fn.Append(Form("_%lldM", fNEvents/1000000));
 308           else if (fNEvents >= 1000)
 309             fn.Append(Form("_%lldk", fNEvents/1000));
 310           else
 311             fn.Append(Form("_%lld", fNEvents));
 312         }
 313         fn.Append(".root");
 314         fFileName = fn;
 315       }
 316     }
 317     return fn.Data();
 318     /*
 319       if (fFileName.IsNull())
 320       fFileName = Form("%s_%09d.root", fEGName.Data(), fRunNo);
 321       return fFileName.Data();*/
 322   }
 323   const char* GetName() const { return "FastSim"; }
 324   const char* GetTitle() const { return "ALICE Event Generator simulation"; }
 325   /**
 326    * Create our outputs
 327    *
 328    *
 329    * @return true on success
 330    */
 331   Bool_t SetupOutput()
 332   {
 333     Info("SetupOutput", "First the file");
 334     Bool_t isProof = false;
 335     if (fInput && fInput->FindObject("PROOF_Ordinal"))
 336       isProof = true;
 337     if (isProof) {
 338       Info("SetupOutput", "Making Proof File");
 339       fProofFile = new TProofOutputFile(FileName(), "M");
 340       // TProofOutputFile::kMerge,
 341       // TProofOutputFile::kRemote);
 342       fFile = fProofFile->OpenFile("RECREATE");
 343     }
 344     else
 345       fFile = TFile::Open(FileName(), "RECREATE");
 346
 347     Info("SetupOutput", "Making our tree");
 348     fTree      = new TTree("T", "T");
 349     fParticles = new TClonesArray("TParticle");
 350     fTree->Branch("header", &fShortHead,
 351                   "run/i:event:npart:nbin:type:ipx/D:ipy:ipz:b:c:phir");
 352     fTree->Branch("particles", &fParticles);
 353     fTree->AutoSave();
 354     fTree->SetDirectory(fFile);
 355     fTree->SetAlias("primary", "(particles.fBits&(1<<14))");
 356     fTree->SetAlias("weak",    "(particles.fBits&(1<<15))");
 357     fTree->SetAlias("charged", "(particles.fBits&(1<<16))");
 358     fTree->SetAlias("pt",      "(sqrt(pow(particles.fPx,2)+"
 359                     /*       */"pow(particles.fPy,2)))");
 360     fTree->SetAlias("eta",     "(pt<1e-10?1024:"
 361                     "-log(tan(atan2(particles.Pt(),particles.fPz)/2)))");
 362     fTree->SetAlias("good",    "(primary&&charged&&abs(eta)<1000)");
 363     fTree->SetAlias("sd",      "(header.fType & 0x1)");
 364     fTree->SetAlias("dd",      "(header.fType & 0x2)");
 365     fTree->SetAlias("pion",    "(abs(particles.fPdgCode)==211)");
 366     fTree->SetAlias("kaon",    "(abs(particles.fPdgCode)==321)");
 367     fTree->SetAlias("proton",  "(abs(particles.fPdgCode)==2212)");
 368     fTree->SetAlias("electron","(abs(particles.fPdgCode)==11)");
 369     fTree->SetAlias("other",   "(!pion&&!kaon&&!proton&&!electron)");
 370     fTree->SetAlias("beta",    "(particles.P()/particle.Energy())");
 371     fTree->SetAlias("gamma",   "(1./sqrt(1-beta*beta))");
 372
 373     Info("SetupOutput", "Making histograms");
 374     Double_t maxEta = 10;
 375     Double_t dEta   = 10./200;
 376     fHEta = new TH1D("dNdeta", "Charged particle pseudo-rapidity density",
 377                      Int_t(2*maxEta/dEta+.5), -maxEta, +maxEta);
 378     fHEta->Sumw2();
 379     fHEta->SetXTitle("#it{#eta}");
 380     fHEta->SetYTitle("1/N d#it{N}_{ch}/d#it{#eta}");
 381     fHEta->SetMarkerColor(kRed+2);
 382     fHEta->SetMarkerStyle(20);
 383     fHEta->SetDirectory(0);
 384
 385     fHIpz = new TH1D("ipZ", "Z-coordinate of interaction point",
 386                      10, -10, 10);
 387     fHIpz->SetMarkerColor(kGreen+2);
 388     fHIpz->SetFillColor(kGreen+2);
 389     fHIpz->SetFillStyle(3001);
 390     fHIpz->SetXTitle("IP_{#it{z}} [cm]");
 391     fHIpz->SetDirectory(0);
 392
 393     fHType = new TH1D("type", "Diffractive", 3, .5, 3.5);
 394     fHType->SetMarkerColor(kOrange+2);
 395     fHType->SetFillColor(kOrange+2);
 396     fHType->SetFillStyle(3001);
 397     fHType->SetDirectory(0);
 398     fHType->GetXaxis()->SetBinLabel(1, "Non");
 399     fHType->GetXaxis()->SetBinLabel(2, "Single");
 400     fHType->GetXaxis()->SetBinLabel(3, "Double");
 401
 402     fHCent = new TH1D("cent", "Centrality", 20, 0, 100);
 403     fHCent->SetMarkerColor(kPink+2);
 404     fHCent->SetFillColor(kPink+2);
 405     fHCent->SetFillStyle(3001);
 406     fHCent->SetDirectory(0);
 407     fHCent->SetXTitle("Centrality [%]");
 408
 409     fHB = new TH1D("b", "Impact parameter", 20, 0, 20);
 410     fHB->SetMarkerColor(kYellow+2);
 411     fHB->SetFillColor(kYellow+2);
 412     fHB->SetFillStyle(3001);
 413     fHB->SetXTitle("#it{b} [fm]");
 414     fHB->SetDirectory(0);
 415
 416     fHPhiR = new TH1D("phiR", "Event plane angle", 360, 0, 360);
 417     fHPhiR->SetMarkerColor(kMagenta+2);
 418     fHPhiR->SetFillColor(kMagenta+2);
 419     fHPhiR->SetFillStyle(3001);
 420     fHPhiR->SetXTitle("#it{#Phi}_{R} [degrees]");
 421     fHPhiR->SetDirectory(0);
 422
 423     fHTime = new TH1D("timing", "Timing of processing", 5,0.5,5.5);
 424     fHTime->SetMarkerColor(kBlue+2);
 425     fHTime->SetFillColor(kBlue+2);
 426     fHTime->SetFillStyle(3001);
 427     fHTime->GetXaxis()->SetBinLabel(1,"Reset");
 428     fHTime->GetXaxis()->SetBinLabel(2,"Generate");
 429     fHTime->GetXaxis()->SetBinLabel(3,"Header");
 430     fHTime->GetXaxis()->SetBinLabel(4,"Particles");
 431     fHTime->GetXaxis()->SetBinLabel(5,"Filling");
 432     fHTime->SetDirectory(0);
 433
 434     fList = new TList;
 435     fList->SetName("histograms");
 436     fList->SetOwner(true);
 437     fList->Add(fHEta);
 438     fList->Add(fHIpz);
 439     fList->Add(fHType);
 440     fList->Add(fHCent);
 441     fList->Add(fHB);
 442     fList->Add(fHPhiR);
 443     fList->Add(fHTime);
 444
 445     Info("SetupOutput", "Adding list ot outputs");
 446     fOutput->Add(fList);
 447
 448     return true;
 449   }
 450   Bool_t SetupGen()
 451   {
 452     Printf(" === Setup ==============================");
 453     Printf("  Run #:                          %6d", fRunNo);
 454     Printf("  EG:     %30s", fEGName.Data());
 455     Printf("  B range:             %5.1ffm - %5.1ffm", fBMin, fBMax);
 456     Printf(" ========================================");
 457     Printf("Macro path: %s", gROOT->GetMacroPath());
 458
 459     // --- Check if we shoud get the GRP line ------------------------
 460     if (!fGRP && fInput) {
 461       fGRP = fInput->FindObject("GRP");
 462       std::ofstream* pout = new std::ofstream("grp.dat");
 463       if (pout) {
 464         Info("SetupGen", "Writing GRP line '%s' to \"grp.dat\"",
 465              fGRP->GetTitle());
 466         std::ostream& out = *pout;
 467         out << fGRP->GetTitle() << std::endl;
 468         pout->close();
 469       }
 470     }
 471
 472     // --- Load our settings -----------------------------------------
 473     Info("SetupGen", "Loading scripts");
 474     gROOT->Macro(Form("GRP.C(%d)", fRunNo));
 475     gROOT->Macro("BaseConfig.C");
 476     gROOT->Macro("EGConfig.C");
 477
 478     gROOT->ProcessLine(Form("VirtualEGCfg::LoadGen(\"%s\")",fEGName.Data()));
 479
 480     // --- Make our generator ----------------------------------------
 481     Info("SetupGen", "Creating generator");
 482     TString egMk = Form("egCfg->MakeGenerator(\"%s\",%f,%f)",
 483                         fEGName.Data(), fBMin, fBMax);
 484     Long64_t egPtr = gROOT->ProcessLine(egMk);
 485     if (egPtr == 0) {
 486       Error("Setup", "Failed to make generator");
 487       return false;
 488     }
 489     fGenerator = reinterpret_cast<AliGenerator*>(egPtr);
 490
 491
 492     if (fFileName.IsNull()) FileName();
 493     Info("SetupRun", "File name is '%s'", fFileName.Data());
 494
 495     return true;
 496   }
 497   /**
 498    * Setup the generator etc. of the job
 499    *
 500    * @param nev Maximum number of events per file
 501    *
 502    * @return true on success
 503    */
 504   Bool_t SetupRun()
 505   {
 506     // --- gAlice (bare ROOT) ----------------------------------------
 507     if (!gAlice)
 508       new AliRun("gAlice", "The ALICE Off-line framework");
 509
 510     Long64_t nev = (fNEvents <= 0 ? 0xFFFFFFFF : fNEvents);
 511     // --- Run-loader, stack, etc  -----------------------------------
 512     Info("SetupRun", "Set-up run Loader");
 513     fRunLoader = AliRunLoader::Open("galice.root", "FASTRUN", "RECREATE");
 514     fRunLoader->SetCompressionLevel(2);
 515     fRunLoader->SetNumberOfEventsPerFile(nev);
 516     fRunLoader->LoadKinematics("RECREATE");
 517     fRunLoader->MakeTree("E");
 518     gAlice->SetRunLoader(fRunLoader);
 519     fRunLoader->MakeStack();
 520     fStack  = fRunLoader->Stack();
 521     fHeader = fRunLoader->GetHeader();
 522
 523     // --- Initialize generator --------------------------------------
 524     Info("SetupRun", "Initializing generator");
 525     fGenerator->Init();
 526     fGenerator->SetStack(fStack);
 527
 528     return true;
 529   }
 530   /**
 531    * Read the previously created grp.dat file
 532    *
 533    */
 534   Bool_t ReadGRPLine()
 535   {
 536     std::ifstream* pin = new std::ifstream("grp.dat");
 537     if (!pin) {
 538       Warning("ReadGRPLine", "Failed to open \"grp.dat\"");
 539       return false;
 540     }
 541     std::istream&  in  = *pin;
 542     TString line;
 543     TString env;
 544     do {
 545       line.ReadLine(in);
 546       if (line.IsNull()) continue;
 547       if (line.BeginsWith("#")) continue;
 548       env = line;
 549       break;
 550     } while (!in.eof());
 551     pin->close();
 552
 553     if (env.IsNull()) {
 554       Warning("ReadGRPLine", "Got no line from \"grp.dat\"");
 555       return false;
 556     }
 557
 558     fGRP = new TNamed("GRP",env.Data());
 559     return true;
 560   }
 561
 562   /**
 563    * Set up job
 564    *
 565    */
 566   void Init(TTree*)
 567   {
 568   }
 569   /**
 570    * Set up job
 571    *
 572    */
 573   void Begin(TTree*)
 574   {
 575     // Make a monitor
 576     Info("Begin", "gProof=%p Nomonitor=%p",
 577          gProof, (gProof ? gProof->GetParameter("NOMONITOR") : 0));
 578
 579     if (gProof && !gProof->GetParameter("NOMONITOR")) {
 580       new FastMonitor;
 581       gProof->AddFeedback("histograms");
 582       Info("Begin", "Adding monitoring");
 583     }
 584     gROOT->Macro(Form("GRP.C(%d)", fRunNo));
 585     if (ReadGRPLine()) {
 586       if(gProof) {
 587         gProof->AddInput(fGRP);
 588       }
 589     }
 590   }
 591   /**
 592    * Set-up this sub-job
 593    *
 594    */
 595   void SlaveBegin(TTree*)
 596   {
 597     SetupGen();
 598     SetupOutput();
 599     SetupRun();
 600   }
 601   /* Reset internal caches etc.
 602    *
 603    * @param iEv Event number
 604    *
 605    * @return true on success
 606    */
 607   Bool_t PreEvent(Long64_t iEv)
 608   {
 609     // --- Reset header ----------------------------------------------
 610     fShortHead.fRunNo   = fRunNo;
 611     fShortHead.fEventId = iEv;
 612     fShortHead.fIpX     = 1024;
 613     fShortHead.fIpY     = 1024;
 614     fShortHead.fIpZ     = 1024;
 615     fShortHead.fNpart   = -1;
 616     fShortHead.fNbin    = -1;
 617     fShortHead.fPhiR    = -1;
 618     fShortHead.fB       = -1;
 619     fShortHead.fC       = -1;
 620     fParticles->Clear();
 621     // --- Reset header, etc.  ---------------------------------------
 622     fHeader->Reset(fRunNo, iEv);
 623     fRunLoader->SetEventNumber(iEv);
 624     fStack->Reset();
 625     fRunLoader->MakeTree("K");
 626
 627     return true;
 628   }
 629   /**
 630    * Process the event header
 631    *
 632    * @return true if the event should be diagnosed
 633    */
 634   Bool_t ProcessHeader()
 635   {
 636     // --- Copy to short header --------------------------------------
 637     fShortHead.fRunNo   = fHeader->GetRun();
 638     fShortHead.fEventId = fHeader->GetEvent();
 639     TArrayF ip;
 640     fHeader->GenEventHeader()->PrimaryVertex(ip);
 641     fShortHead.fIpX     = ip[0];
 642     fShortHead.fIpY     = ip[1];
 643     fShortHead.fIpZ     = ip[2];
 644
 645     // --- Check header type -----------------------------------------
 646     AliGenEventHeader* genHeader = fHeader->GenEventHeader();
 647     AliCollisionGeometry* geometry =
 648       dynamic_cast<AliCollisionGeometry*>(genHeader);
 649     AliGenPythiaEventHeader* pythia    =
 650       dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 651     AliGenDPMjetEventHeader* dpm       =
 652       dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
 653     AliGenGeVSimEventHeader* gev       =
 654       dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
 655     AliGenHerwigEventHeader* herwig    =
 656       dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
 657     if (geometry) {
 658       fShortHead.fB     = geometry->ImpactParameter();
 659       fShortHead.fNpart = (geometry->ProjectileParticipants() +
 660                            geometry->TargetParticipants());
 661       fShortHead.fNbin  = geometry->NN();
 662       fShortHead.fPhiR  = geometry->ReactionPlaneAngle();
 663     }
 664     // --- Determine diffraction flags -------------------------------
 665     Bool_t sd = false;
 666     Bool_t dd = false;
 667     if (pythia) {
 668       Int_t type = pythia->ProcessType();
 669       if (type < 100) { // pythia6
 670         switch (type) {
 671         case 92: case 93: sd = true; break;
 672         case 94:          dd = true; break;
 673         }
 674       }
 675       else {
 676         switch (type) { // Pythia8
 677         case 103: case 104: sd = true; break;
 678         case 105:           dd = true; break;
 679         }
 680       }
 681       fShortHead.fB     = pythia->GetImpactParameter();
 682       fShortHead.fNpart = 2;
 683       fShortHead.fNbin  = 1;
 684     }
 685     if (dpm) {
 686       Int_t type = dpm->ProcessType();
 687       switch (type) {
 688       case 5: case 6: sd = true;
 689
 690       case 7:         dd = true;
 691       }
 692     }
 693     if (gev) fShortHead.fPhiR = gev->GetEventPlane();
 694     if (herwig) {
 695       Int_t type = herwig->ProcessType();
 696       switch (type) {
 697       case 5: case 6: sd = true; break;
 698       }
 699       fShortHead.fNpart = 2;
 700       fShortHead.fNbin  = 1;
 701     }
 702     fShortHead.fType = (sd ? 0x1 : 0) | (dd ? 0x2 : 0);
 703
 704     // --- Check centrality -----------------------------------------
 705     // PbPb only
 706     // Updated 4th of November 2014 from
 707     // cern.ch/twiki/bin/view/ALICE/CentStudies#Tables_with_centrality_bins_AN1
 708     Float_t  np = 0;
 709     UInt_t   nc = 0;
 710     Double_t c  = -1;
 711     Double_t b  = fShortHead.fB;
 712     if (b >= 0) {
 713       if      (0.00 >= b  && b < 1.57)  { c=0.5;  np=403.8; nc=1861; }
 714       else if (1.57 >= b  && b < 2.22)  { c=1.5;  np=393.6; nc=1766; }
 715       else if (2.22 >= b  && b < 2.71)  { c=2.5;  np=382.9; nc=1678; }
 716       else if (2.71 >= b  && b < 3.13)  { c=3.5;  np=372;   nc=1597; }
 717       else if (3.13 >= b  && b < 3.50)  { c=4.5;  np=361.1; nc=1520; }
 718       else if (3.50 >= b  && b < 4.94)  { c=7.5;  np=329.4; nc=1316; }
 719       else if (4.94 >= b  && b < 6.05)  { c=12.5; np=281.2; nc=1032; }
 720       else if (6.05 >= b  && b < 6.98)  { c=17.5; np=239;   nc=809.8; }
 721       else if (6.98 >= b  && b < 7.81)  { c=22.5; np=202.1; nc=629.6; }
 722       else if (7.81 >= b  && b < 8.55)  { c=27.5; np=169.5; nc=483.7; }
 723       else if (8.55 >= b  && b < 9.23)  { c=32.5; np=141;   nc=366.7; }
 724       else if (9.23 >= b  && b < 9.88)  { c=37.5; np=116;   nc=273.4; }
 725       else if (9.88 >= b  && b < 10.47) { c=42.5; np=94.11; nc=199.4; }
 726       else if (10.47 >= b && b < 11.04) { c=47.5; np=75.3;  nc=143.1; }
 727       else if (11.04 >= b && b < 11.58) { c=52.5; np=59.24; nc=100.1; }
 728       else if (11.58 >= b && b < 12.09) { c=57.5; np=45.58; nc=68.46; }
 729       else if (12.09 >= b && b < 12.58) { c=62.5; np=34.33; nc=45.79; }
 730       else if (12.58 >= b && b < 13.05) { c=67.5; np=25.21; nc=29.92; }
 731       else if (13.05 >= b && b < 13.52) { c=72.5; np=17.96; nc=19.08; }
 732       else if (13.52 >= b && b < 13.97) { c=77.5; np=12.58; nc=12.07; }
 733       else if (13.97 >= b && b < 14.43) { c=82.5; np=8.812; nc=7.682; }
 734       else if (14.43 >= b && b < 14.96) { c=87.5; np=6.158; nc=4.904; }
 735       else if (14.96 >= b && b < 15.67) { c=92.5; np=4.376; nc=3.181; }
 736       else if (15.67 >= b && b < 20.00) { c=97.5; np=3.064; nc=1.994; }
 737       fShortHead.fC = c;
 738       // Be careful to round off
 739       if (fShortHead.fNpart <= 0) fShortHead.fNpart = Int_t(np+.5);
 740       if (fShortHead.fNbin  <= 0) fShortHead.fNbin  = Int_t(nc+.5)/2;
 741     }
 742
 743     // --- Check if within vertex cut -------------------------------
 744     Bool_t selected = (fShortHead.fIpZ <= fHIpz->GetXaxis()->GetXmax() &&
 745                        fShortHead.fIpZ >= fHIpz->GetXaxis()->GetXmin());
 746
 747     // --- Only update histograms if within IPz cut ------------------
 748     if (selected) {
 749       fHPhiR->Fill(fShortHead.fPhiR*TMath::RadToDeg());
 750       fHB->Fill(fShortHead.fB);
 751       fHIpz->Fill(fShortHead.fIpZ);
 752       fHType->Fill(dd ? 3 : sd ? 2 : 1);
 753       fHCent->Fill(c);
 754     }
 755     return selected;
 756   }
 757   /**
 758    * Process all particles
 759    *
 760    * @param selected True if particle information should be diagnosed
 761    *
 762    * @return true on success
 763    */
 764   Bool_t ProcessParticles(Bool_t selected)
 765   {
 766     Int_t nPart = fStack->GetNprimary();
 767     for (Int_t iPart = 0; iPart < nPart; iPart++) {
 768       TParticle*    particle  = fStack->Particle(iPart);
 769       TParticlePDG* pdg       = particle->GetPDG();
 770       Bool_t        primary   = fStack->IsPhysicalPrimary(iPart);
 771       Bool_t        weakDecay = fStack->IsSecondaryFromWeakDecay(iPart);
 772       Bool_t        charged   = (pdg &&  TMath::Abs(pdg->Charge()) > 0);
 773       if (primary)   particle->SetBit(BIT(14));
 774       if (weakDecay) particle->SetBit(BIT(15));
 775       if (charged)   particle->SetBit(BIT(16));
 776
 777       new ((*fParticles)[iPart]) TParticle(*particle);
 778
 779       if (!selected || !charged || !primary) continue;
 780       Double_t pT    = particle->Pt();
 781       if (pT < 1e-10) continue; /// Along beam axis
 782       Double_t pZ    = particle->Pz();
 783       Double_t theta = TMath::ATan2(pT, pZ);
 784       Double_t eta   = -TMath::Log(TMath::Tan(theta/2));
 785       fHEta->Fill(eta);
 786     }
 787     return true;
 788   }
 789   /**
 790    * Do final event processing (fill output)
 791    *
 792    */
 793   void PostEvent()
 794   {
 795     fHeader->SetNprimary(fStack->GetNprimary());
 796     fHeader->SetNtrack(fStack->GetNtrack());
 797
 798     fTree->Fill();
 799
 800     fStack->FinishEvent();
 801     fHeader->SetStack(fStack);
 802
 803     fRunLoader->TreeE()->Fill();
 804     fRunLoader->WriteKinematics("OVERWRITE");
 805   }
 806   /**
 807    * Process one event
 808    *
 809    * @param iEv Event number
 810    *
 811    * @return true on success, false otherwize
 812    */
 813   Bool_t Process(Long64_t iEv)
 814   {
 815     // --- The stopwatch ---------------------------------------------
 816     TStopwatch timer;
 817     timer.Start();
 818     PreEvent(iEv);
 819     fHTime->Fill(1, timer.RealTime());
 820
 821     // --- Generate event --------------------------------------------
 822     timer.Start();
 823     fGenerator->Generate();
 824     fHTime->Fill(2, timer.RealTime());
 825
 826     // --- Process the header ----------------------------------------
 827     timer.Start();
 828     Bool_t selected = ProcessHeader();
 829     fHTime->Fill(3, timer.RealTime());
 830
 831     // --- Loop over particles ---------------------------------------
 832     timer.Start();
 833     ProcessParticles(selected);
 834     fHTime->Fill(4, timer.RealTime());
 835
 836     // --- Do final stuff --------------------------------------------
 837     timer.Start();
 838     PostEvent();
 839     fHTime->Fill(5, timer.RealTime());
 840
 841     return true;
 842   }
 843   /**
 844    * Finalize this sub-job
 845    *
 846    */
 847   void SlaveTerminate()
 848   {
 849     fGenerator->FinishRun();
 850     fRunLoader->WriteHeader("OVERWRITE");
 851     fGenerator->Write();
 852     fRunLoader->Write();
 853
 854     if (fFile) {
 855       if (fProofFile) {
 856         fOutput->Add(fProofFile);
 857         fOutput->Add(new TH1F("filename", fFileName.Data(),1,0,1));
 858       }
 859       // Flush out tree
 860       fFile->cd();
 861       fTree->Write(0, TObject::kOverwrite);
 862       fFile->Close();
 863       fFile->Delete();
 864       fFile = 0;
 865     }
 866   }
 867   /**
 868    * Final processing of the data
 869    *
 870    */
 871   void Terminate()
 872   {
 873     if (gProof) gProof->ClearFeedback();
 874
 875     if (!fList)
 876       fList = static_cast<TList*>(fOutput->FindObject("histograms"));
 877     if (!fList) {
 878       Error("Terminate", "No output list");
 879       return;
 880     }
 881
 882     if (!fProofFile) {
 883       TObject* fn = fOutput->FindObject("filename");
 884       if (fn) fFileName  = fn->GetTitle();
 885       fProofFile =
 886         static_cast<TProofOutputFile*>(fOutput->FindObject(FileName()));
 887     }
 888     if (fProofFile)
 889       fFile = fProofFile->OpenFile("UPDATE");
 890     if (!fFile)
 891       fFile = TFile::Open(FileName(),"UPDATE");
 892
 893
 894     fHEta  = static_cast<TH1*>(fList->FindObject("dNdeta"));
 895     fHIpz  = static_cast<TH1*>(fList->FindObject("ipZ"));
 896     fHType = static_cast<TH1*>(fList->FindObject("type"));
 897     fHCent = static_cast<TH1*>(fList->FindObject("cent"));
 898     fHB    = static_cast<TH1*>(fList->FindObject("b"));
 899     fHPhiR = static_cast<TH1*>(fList->FindObject("phiR"));
 900     fHTime = static_cast<TH1*>(fList->FindObject("timing"));
 901
 902     if (!(fHEta && fHIpz && fHType && fHB && fHPhiR && fHTime)) {
 903       Warning("Terminate", "Missing histograms (%p,%p,%p,%p,%p,%p)",
 904               fHEta, fHIpz, fHType, fHB, fHPhiR, fHTime);
 905       return;
 906     }
 907
 908     Int_t nTotal = fHIpz->GetEntries();
 909     fHEta ->Scale(1./nTotal, "width");
 910     fHB   ->Scale(1./nTotal, "width");
 911     fHPhiR->Scale(1./nTotal, "width");
 912     fHTime->Scale(1./nTotal, "width");
 913
 914     if (!fFile){
 915       Warning("Terminate", "No file to write to");
 916       return;
 917     }
 918
 919     fHEta ->Write();
 920     fHIpz ->Write();
 921     fHType->Write();
 922     fHCent->Write();
 923     fHB   ->Write();
 924     fHPhiR->Write();
 925     fHTime->Write();
 926
 927     fTree = static_cast<TTree*>(fFile->Get("T"));
 928     if (!fTree)  Warning("Terminate", "No tree");
 929
 930     fFile->Close();
 931   }
 932   /**
 933    * Interface version used
 934    *
 935    * @return 1
 936    */
 937   Int_t Version() const { return 1; }
 938
 939   /**
 940    * @{
 941    * @name Parameters
 942    */
 943   TString  fEGName;               // Name of event generator
 944   Int_t    fRunNo;                // Run to simulate
 945   Double_t fBMin;                 // Least impact parameter
 946   Double_t fBMax;                 // Largest impact parameter
 947   TObject* fGRP;                  //! GRP in one line
 948   Long64_t fNEvents;              //  Number of requested events
 949   /* @} */
 950   /**
 951    * @{
 952    * @name ALICE EG interface
 953    */
 954   AliGenerator* fGenerator;       //! Event generator
 955   AliRunLoader* fRunLoader;       //! Loader of trees
 956   AliStack*     fStack;           //! Stack of particles
 957   AliHeader*    fHeader;          //! Header handler
 958   /* @} */
 959   /**
 960    * @{
 961    * @name Custom output
 962    */
 963   TTree*        fTree;            //! Custom tree
 964   TClonesArray* fParticles;       //! List of particles
 965   /**
 966    * @{
 967    * @name Diagnostics
 968    */
 969   TList* fList;                   //! List of outputs
 970   TH1*   fHEta;                   //! dN/deta
 971   TH1*   fHIpz;                   //! IPz histogram
 972   TH1*   fHType;                  //! Event type histogram
 973   TH1*   fHCent;                  //! Event type histogram
 974   TH1*   fHB;                     //! B histogram
 975   TH1*   fHPhiR;                  //! Reaction plane
 976   TH1*   fHTime;                  //! Timing
 977   /* @} */
 978   /**
 979    * @{
 980    * @name Output files
 981    */
 982   TProofOutputFile* fProofFile;   //! Proof output file
 983   TFile*            fFile;        //! Output file
 984   mutable TString   fFileName;    //! Output file name
 985   /* @} */
 986
 987   // Hide from CINT
 988 #ifndef __CINT__
 989   struct ShortHeader {
 990     UInt_t   fRunNo;
 991     UInt_t   fEventId;
 992     UInt_t   fNpart;
 993     UInt_t   fNbin;
 994     UInt_t   fType;
 995     Double_t fIpX;
 996     Double_t fIpY;
 997     Double_t fIpZ;
 998     Double_t fB;
 999     Double_t fC;
1000     Double_t fPhiR;
1001   } fShortHead;
1002 #endif
1003   /**
1004    * Run this selector as a normal process
1005    *
1006    * @param nev        Number of events
1007    * @param run        Run number to anchor in
1008    * @param gen        Generator
1009    * @param bMin       Least impact parameter [fm]
1010    * @param bMax       Largest impact parameter [fm]
1011    * @param monitor    Monitor frequency [s]
1012    *
1013    * @return true on succes
1014    */
1015   static Bool_t  Run(Long64_t    nev,
1016                      UInt_t      run,
1017                      const char* gen,
1018                      Double_t    bMin,
1019                      Double_t    bMax,
1020                      Int_t       monitor)
1021   {
1022     TStopwatch stopwatch;
1023     stopwatch.Start();
1024
1025     FastSim* sim = new FastSim(gen,run,bMin,bMax,nev);
1026     sim->Begin(0);
1027     sim->SlaveBegin(0);
1028
1029     TTimer* timer = 0;
1030     if (monitor > 0) {
1031       // timer = new TTimer(new FastMonitor(sim), monitor*1000,true);
1032       timer = new TTimer(1000);
1033       timer->Connect("Timeout()","FastMonitor",
1034                      new FastMonitor(sim), "Handle()");
1035       ::Info("Run", "Turning on monitoring");
1036       timer->Start(-1,false);
1037     }
1038
1039     for (Long64_t i=0; i <nev; i++) {
1040       Printf("=== Event # %6lld/%6lld ==========================",
1041              i+1, nev);
1042       sim->Process(i);
1043       if (timer && (i > 0) && (i % 500 == 0)) {
1044         if (timer->CheckTimer(gSystem->Now()))
1045           Printf("Fired timer");
1046       }
1047     }
1048     if (timer) timer->TurnOff();
1049     sim->SlaveTerminate();
1050     sim->Terminate();
1051
1052     stopwatch.Print();
1053     return true;
1054   }
1055   static void ProofLoadLibs()
1056   {
1057     if (!gProof) return;
1058
1059     // Remember to copy changes to RunFast.C
1060     TList clsLib;
1061     clsLib.Add(new TNamed("TVirtualMC",              "libVMC"));
1062     clsLib.Add(new TNamed("TLorentzVector",          "libPhysics"));
1063     clsLib.Add(new TNamed("TLinearFitter",           "libMinuit"));
1064     clsLib.Add(new TNamed("TTree",                   "libTree"));
1065     clsLib.Add(new TNamed("TProof",                  "libProof"));
1066     clsLib.Add(new TNamed("TGFrame",                 "libGui"));
1067     clsLib.Add(new TNamed("TSAXParser",              "libXMLParser"));
1068     clsLib.Add(new TNamed("AliVEvent",               "libSTEERBase"));
1069     clsLib.Add(new TNamed("AliESDEvent",             "libESD"));
1070     clsLib.Add(new TNamed("AliAODEvent",             "libAOD"));
1071     clsLib.Add(new TNamed("AliAnalysisManager",      "libANALYSIS"));
1072     clsLib.Add(new TNamed("AliCDBManager",           "libCDB"));
1073     clsLib.Add(new TNamed("AliRawVEvent",            "libRAWDatabase"));
1074     clsLib.Add(new TNamed("AliHit",                  "libSTEER"));
1075     clsLib.Add(new TNamed("AliGenMC",                "libEVGEN"));
1076     clsLib.Add(new TNamed("AliFastEvent",            "libFASTSIM"));
1077
1078     TIter next(&clsLib);
1079     TObject* obj = 0;
1080     while ((obj = next())) {
1081       gProof->Exec(Form("gROOT->LoadClass(\"%s\",\"%s\");",
1082                         obj->GetName(), obj->GetTitle()));
1083     }
1084   }
1085   /**
1086    * Run this selector in PROOF(Lite)
1087    *
1088    * @param url        Proof URL
1089    * @param nev        Number of events
1090    * @param run        Run number to anchor in
1091    * @param gen        Generator
1092    * @param bMin       Least impact parameter [fm]
1093    * @param bMax       Largest impact parameter [fm]
1094    * @param monitor    Monitor frequency [s]
1095    * @param opt        Compilation options
1096    *
1097    * @return true on succes
1098    */
1099   static Bool_t Proof(const char*  url,
1100                       Long64_t     nev,
1101                       UInt_t       run,
1102                       const char*  gen,
1103                       Double_t     bMin,
1104                       Double_t     bMax,
1105                       Int_t        monitor=-1,
1106                       const char*  opt="")
1107   {
1108     Printf("# events:  %lld", nev);
1109     Printf("Run #:     %u",   run);
1110     Printf("Generator: %s",   gen);
1111     Printf("b range:   %5.1f-%5.1f", bMin, bMax);
1112     Printf("monitor:   %ds", monitor);
1113
1114     TStopwatch timer;
1115     timer.Start();
1116
1117     TProof::Reset(url);
1118     TProof::Open(url);
1119     gProof->ClearCache();
1120
1121     TString ali = gSystem->ExpandPathName("$(ALICE_ROOT)");
1122     // TString fwd = gSystem->ExpandPathName("$ANA_SRC");
1123     TString fwd = ali + "/PWGLF/FORWARD/analysis2";
1124
1125     gProof->AddIncludePath(Form("%s/include", ali.Data()));
1126     ProofLoadLibs();
1127     gProof->Load(Form("%s/sim/GRP.C",fwd.Data()), true);
1128     gProof->Load(Form("%s/sim/BaseConfig.C",fwd.Data()), true);
1129     gProof->Load(Form("%s/sim/EGConfig.C",fwd.Data()), true);
1130
1131     // gROOT->ProcessLine("gProof->SetLogLevel(5);");
1132     gProof->Load(Form("%s/sim/FastSim.C+%s", fwd.Data(), opt),true);
1133
1134     if (monitor <= 0) gProof->SetParameter("NOMONITOR", true/*ignored*/);
1135     else              gProof->SetParameter("PROOF_FeedbackPeriod",
1136                                            monitor*1000/*ms*/);
1137
1138     FastSim* sim = new FastSim(gen,run,bMin,bMax,nev);
1139     gProof->Process(sim, nev, "");
1140
1141     timer.Print();
1142     return true; // status >= 0;
1143   }
1144   ClassDef(FastSim,1);
1145 };
1146
1147 #endif
1148 //
1149 // EOF
1150 //