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