AliPhysics  68dfc25 (68dfc25)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FastSim.C
Go to the documentation of this file.
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 "AliGenHijingEventHeader.h"
18 # include "AliGenCocktailEventHeader.h"
19 # include "AliPDG.h"
20 # include <TROOT.h>
21 # include <TString.h>
22 # include <TMath.h>
23 # include <TParticle.h>
24 # include <TH1.h>
25 # include <TTree.h>
26 # include <TClonesArray.h>
27 # include <TList.h>
28 # include <TProof.h>
29 # include <TChain.h>
30 # include <TParticlePDG.h>
31 # include <TStopwatch.h>
32 # include <TFile.h>
33 # include <TProofOutputFile.h>
34 # include <TCanvas.h>
35 # include <TTimer.h>
36 # include <TRandom.h>
37 # include <TUrl.h>
38 # include <TMacro.h>
39 # include <TSystemDirectory.h>
40 # include <TFileCollection.h>
41 # include <TPRegexp.h>
42 # include <fstream>
43 # include <TLeaf.h>
44 # include <algorithm>
45 # include "FastShortHeader.C"
46 # include "FastCentEstimators.C"
47 # include "FastMonitor.C"
48 #else
49 class AliGenerator;
50 class AliRunLoader;
51 class AliStack;
52 class AliHeader;
53 class AliGenEventHeader;
54 class TH1;
55 class TTree;
56 class TClonesArray;
57 class TBrowser;
58 class TList;
59 class TFile;
60 class TProofOutputFile;
61 class TCanvas;
62 class TVirtualPad;
63 class TTimer;
64 class TUrl;
65 class TAxis;
66 class TParticle;
67 class TMacro;
68 class TLeaf;
69 class BCentEstimator;
70 class FastMonitor;
71 #endif
72 
73 
74 
76 typedef struct {
77  Double_t rproj;
78  Double_t rtarg;
79  Double_t bimpac;
80  Int_t nwtsam;
81  Int_t nwasam;
82  Int_t nwbsam;
83  Int_t nwtacc;
84  Int_t nwtaac;
85  Int_t nwtbac;
86  Int_t ncp;
87  Int_t nct;
88 } DtglcpCommon;
89 DtglcpCommon* _dtglcp = 0;
90 
91 //====================================================================
92 struct FastSimMonitor : public FastMonitor
93 {
94  FastSimMonitor(TSelector* s=0)
95  : FastMonitor(s, "FastMonitor")
96  {
97  Register("list/histograms/b", "hist", 0, false);
98  Register("list/histograms/cent", "hist", 0, false);
99  Register("list/histograms/timing", "hist", 0x18, false);
100  Register("list/histograms/dNdeta", "", 0x18, false);
101  Register("list/histograms/dNdy", "", 0x18, false);
102  Register("list/histograms/trigger", "e", 0x10, false);
103  Register("list/estimators/rawV0MP", "e", 0x02, false);
104  if (gProof) gProof->AddFeedback("list");
105  }
106 };
107 
108 //====================================================================
112 struct FastSim : public TSelector
113 {
124  FastSim(const char* eg="",
125  ULong_t runNo=0,
126  Double_t bMin=0,
127  Double_t bMax=20,
128  Long64_t nEvents=0,
129  Int_t monitor=0)
130  : TSelector(),
131  fEGName(eg),
132  fRunNo(runNo),
133  fBMin(bMin),
134  fBMax(bMax),
135  fGRP(0),
136  fOverrides(0),
137  fNEvents(nEvents),
138  fIsTgtA(false),
139  fIsProjA(false),
140  fGenerator(0),
141  fRunLoader(0),
142  fStack(0),
143  fHeader(0),
144  fTree(0),
145  fParticles(0),
146  fList(0),
147  fHEta(0),
148  fHY(0),
149  fHIpz(0),
150  fHType(0),
151  fHCent(0),
152  fHB(0),
153  fHPhiR(0),
154  fHTime(0),
155  fHTrig(0),
156  fBEstimator(0),
157  fCentEstimators(0),
158  fProofFile(0),
159  fAliceFile(0),
160  fKineFile(0),
161  fFile(0),
162  fFileName(""),
163  fVerbose(true),
164  fMonitor(monitor)
165  {}
166  const char* FileName() const
167  {
168  static TString fn;
169  if (fn.IsNull()) {
170  if (!fFileName.IsNull()) fn = fFileName;
171  else {
172  TString egName = (fGenerator ? fGenerator->GetName() : "");
173  if (egName.IsNull()) egName = fEGName;
174 
175  fn = Form("%s_%09d", egName.Data(), fRunNo);
176  if (fGenerator) {
177  TString tgt, proj;
178  Int_t tgtA, tgtZ, projA, projZ;
179  fGenerator->GetTarget(tgt, tgtA, tgtZ);
180  fGenerator->GetProjectile(proj, projA, projZ);
181  fn.Append(Form("_%s%s", tgt.Data(), proj.Data()));
182  fn.Append(Form("_%05d", Int_t(fGenerator->GetEnergyCMS())));
183  }
184  else {
185  Long_t en = gROOT->ProcessLine("grp->energy");
186  // Bool_t a1 = gROOT->ProcessLine("(GRPData*)%p->beam1.IsA()",fGRP);
187  // Bool_t a2 = gROOT->ProcessLine("(GRPData*)%p->beam2.IsA()",fGRP);
188  fn.Append(Form("_%s%s", (fIsTgtA ? "A" : "p"),
189  (fIsProjA ? "A" : "p")));
190  fn.Append(Form("_%05ld", (fGRP ? en : 0)));
191  }
192 
193  if (fNEvents > 0) {
194  if (fNEvents >= 1000000)
195  fn.Append(Form("_%lldM", fNEvents/1000000));
196  else if (fNEvents >= 1000)
197  fn.Append(Form("_%lldk", fNEvents/1000));
198  else
199  fn.Append(Form("_%lld", fNEvents));
200  }
201  fn.Append(".root");
202  fFileName = fn;
203  }
204  }
205  return fn.Data();
206  }
212  const char* GetName() const { return "FastSim"; }
218  const char* GetTitle() const { return "ALICE Event Generator simulation"; }
224  virtual UShort_t GetSNN() const
225  {
226  Float_t fsNN = 0;
227  if (fGenerator) fsNN = fGenerator->GetEnergyCMS();
228  else fsNN = gROOT->ProcessLine("grp->energy");
229 
230  UShort_t sNN = (TMath::Abs(fsNN - 2760) < 10 ? 2760 :
231  TMath::Abs(fsNN - 5023) < 10 ? 5023 :
232  TMath::Abs(fsNN - 2360) < 10 ? 2360 :
233  TMath::Abs(fsNN - 900) < 10 ? 900 :
234  TMath::Abs(fsNN - 7000) < 10 ? 7000 :
235  TMath::Abs(fsNN - 8000) < 10 ? 8000 :
236  TMath::Abs(fsNN - 13000) < 10 ? 13000 :
237  0);
238  return sNN;
239  }
245  virtual const char* GetEGTitle() const
246  {
247 
248  if (fGenerator) return fGenerator->GetTitle();
249  return "";
250  }
256  Bool_t SetupOutput()
257  {
258  Printf("=== Output =============================\n"
259  " File name: %s\n"
260  "========================================", FileName());
261 
262  if (fVerbose) Info("SetupOutput", "First the file");
263  Bool_t isProof = false;
264  if (fInput && fInput->FindObject("PROOF_Ordinal"))
265  isProof = true;
266  if (isProof) {
267  Info("SetupOutput", "Making Proof File");
268  fProofFile = new TProofOutputFile(FileName(), "M");
269  // TProofOutputFile::kMerge,
270  // TProofOutputFile::kRemote);
271  fFile = fProofFile->OpenFile("RECREATE");
272  }
273  else
274  fFile = TFile::Open(FileName(), "RECREATE");
275 
276  UShort_t sNN = GetSNN();
277  TString tit = GetEGTitle();
278  if (fVerbose) Info("SetupOutput", "Making our tree: %s", tit.Data());
279  fTree = new TTree("T", tit.Data());
280  fParticles = new TClonesArray("TParticle");
281  fTree->Branch("header", &fShortHead,
282  "run/i:event:ntgt:nproj:nbin:type:"
283  "ipx/D:ipy:ipz:b:c:phir:"
284  "nsnp/i:nsnt:nspp:nspt:mask");
285  fTree->Branch("particles", &fParticles);
286  fTree->SetAutoSave(500); // Save every on every 100 events
287  fTree->SetDirectory(fFile);
288  fTree->SetAlias("primary", "(particles.fBits&(1<<14))");
289  fTree->SetAlias("weak", "(particles.fBits&(1<<15))");
290  fTree->SetAlias("charged", "(particles.fBits&(1<<16))");
291  fTree->SetAlias("pt", "(sqrt(pow(particles.fPx,2)+"
292  /* */"pow(particles.fPy,2)))");
293  fTree->SetAlias("eta", "(particles.Pt()<1e-100?"
294  "(particles.fPz>0?100:-100):"
295  "-log(tan(atan2(particles.Pt(),particles.fPz)/2)))");
296  fTree->SetAlias("good", "(primary&&charged&&abs(eta)<1000)");
297  fTree->SetAlias("sd", "(header.fType & 0x1)");
298  fTree->SetAlias("dd", "(header.fType & 0x2)");
299  fTree->SetAlias("pion", "(abs(particles.fPdgCode)==211)");
300  fTree->SetAlias("kaon", "(abs(particles.fPdgCode)==321)");
301  fTree->SetAlias("proton", "(abs(particles.fPdgCode)==2212)");
302  fTree->SetAlias("electron","(abs(particles.fPdgCode)==11)");
303  fTree->SetAlias("other", "(!pion&&!kaon&&!proton&&!electron)");
304  fTree->SetAlias("beta", "(particles.P()/particle.Energy())");
305  fTree->SetAlias("gamma", "(1./sqrt(1-beta*beta))");
306  fTree->SetAlias("npart", "(header.ntgt+header.nproj)");
307  fTree->SetAlias("v0a", "(header.mask&0x1)");
308  fTree->SetAlias("v0c", "(header.mask&0x2)");
309  fTree->SetAlias("ada", "(header.mask&0x4)");
310  fTree->SetAlias("adc", "(header.mask&0x8)");
311  fTree->SetAlias("eta1", "(header.mask&0x10)");
312 
313  // Let's add the title of the generator to the output. We make a
314  // histogram because that can be merged.
315  Info("SetupOutput", "Making generator histogram: %s",tit.Data());
316  TH1* egTitle = new TH1I("eg", tit.Data(), 1, 0, 1);
317  egTitle->SetDirectory(0);
318  egTitle->Fill(0.5);
319  egTitle->SetXTitle("N_{\\hbox{workers}}");
320  egTitle->SetYTitle("Count");
321  egTitle->SetFillColor(kYellow+4);
322  egTitle->SetLineColor(kYellow+4);
323  egTitle->SetMarkerColor(kYellow+4);
324  egTitle->SetFillStyle(1001);
325  egTitle->SetMarkerStyle(1);
326  egTitle->SetLineStyle(1);
327  egTitle->SetLineWidth(2);
328  egTitle->SetMarkerSize(1);
329 
330 
331  if (fVerbose) Info("SetupOutput", "Making histograms");
332  Double_t maxEta = 10;
333  Double_t dEta = 10./200;
334  fHEta = new TH1D("dNdeta", "Charged particle pseudo-rapidity density",
335  Int_t(2*maxEta/dEta+.5), -maxEta, +maxEta);
336  fHEta->Sumw2();
337  fHEta->SetXTitle("\\eta");
338  fHEta->SetYTitle("\\hbox{d}N_{\\hbox{ch}}/\\hbox{d}\\eta");
339  fHEta->SetMarkerColor(kRed+2);
340  fHEta->SetMarkerStyle(20);
341  fHEta->SetDirectory(0);
342 
343  fHY = new TH1D("dNdy", "Charged particle rapidity density",
344  Int_t(2*maxEta/dEta+.5), -maxEta, +maxEta);
345  fHY->Sumw2();
346  fHY->SetXTitle("\\mathit{y}");
347  fHY->SetYTitle("\\hbox{d}N_{\\hbox{ch}}/\\hbox{d}y");
348  fHY->SetMarkerColor(kRed+2);
349  fHY->SetMarkerStyle(20);
350  fHY->SetDirectory(0);
351 
352  fHIpz = new TH1D("ipZ", "Z-coordinate of interaction point",
353  10, -10, 10);
354  fHIpz->SetMarkerColor(kGreen+2);
355  fHIpz->SetFillColor(kGreen+2);
356  fHIpz->SetFillStyle(3001);
357  fHIpz->SetXTitle("IP_{#it{z}} [cm]");
358  fHIpz->SetDirectory(0);
359 
360  fHType = new TH1D("type", "Diffractive", 3, .5, 3.5);
361  fHType->SetMarkerColor(kOrange+2);
362  fHType->SetFillColor(kOrange+2);
363  fHType->SetFillStyle(3001);
364  fHType->SetDirectory(0);
365  fHType->GetXaxis()->SetBinLabel(1, "Non");
366  fHType->GetXaxis()->SetBinLabel(2, "Single");
367  fHType->GetXaxis()->SetBinLabel(3, "Double");
368 
369  fHCent = new TH1D("cent", "Centrality", 20, 0, 100);
370  fHCent->SetMarkerColor(kPink+2);
371  fHCent->SetFillColor(kPink+2);
372  fHCent->SetFillStyle(3001);
373  fHCent->SetDirectory(0);
374  fHCent->SetYTitle("Events");
375  fHCent->SetXTitle("Centrality [%]");
376 
377  fHB = new TH1D("b", "Impact parameter", 20, 0, 20);
378  fHB->SetMarkerColor(kYellow+2);
379  fHB->SetFillColor(kYellow+2);
380  fHB->SetFillStyle(3001);
381  fHB->SetYTitle("Events");
382  fHB->SetXTitle("#it{b} [fm]");
383  fHB->SetDirectory(0);
384 
385  fHPhiR = new TH1D("phiR", "Event plane angle", 360, 0, 360);
386  fHPhiR->SetMarkerColor(kMagenta+2);
387  fHPhiR->SetFillColor(kMagenta+2);
388  fHPhiR->SetFillStyle(3001);
389  fHPhiR->SetXTitle("#it{#Phi}_{R} [degrees]");
390  fHPhiR->SetDirectory(0);
391 
392  fHTime = new TH1D("timing", "Timing of processing", 5,0.5,5.5);
393  fHTime->SetMarkerColor(kBlue+2);
394  fHTime->SetFillColor(kBlue+2);
395  fHTime->SetFillStyle(3001);
396  fHTime->SetYTitle("seconds / event");
397  fHTime->GetXaxis()->SetBinLabel(1,"Reset");
398  fHTime->GetXaxis()->SetBinLabel(2,"Generate");
399  fHTime->GetXaxis()->SetBinLabel(3,"Header");
400  fHTime->GetXaxis()->SetBinLabel(4,"Particles");
401  fHTime->GetXaxis()->SetBinLabel(5,"Filling");
402  fHTime->SetDirectory(0);
403 
404  fHTrig = new TH1D("trigger", "Trigger bits set", 6, -0.5, 5.5);
405  fHTrig->SetMarkerColor(kMagenta+2);
406  fHTrig->SetFillColor(kMagenta+2);
407  fHTrig->SetFillStyle(3001);
408  fHTrig->SetYTitle("Events");
409  fHTrig->GetXaxis()->SetBinLabel(1, "All");
410  fHTrig->GetXaxis()->SetBinLabel(2, "V0A");
411  fHTrig->GetXaxis()->SetBinLabel(3, "V0C");
412  fHTrig->GetXaxis()->SetBinLabel(4, "ADA");
413  fHTrig->GetXaxis()->SetBinLabel(5, "ADC");
414  fHTrig->GetXaxis()->SetBinLabel(6, "N_{ch}|_{|#eta|<1}#ge1");
415  fHTrig->SetDirectory(0);
416 
417  fList = new TList;
418  fList->SetName("list");
419  fList->Add(egTitle);
420 
421  TList* histos = new TList;
422  histos->SetName("histograms");
423  histos->SetOwner(true);
424  histos->Add(fHEta);
425  histos->Add(fHY);
426  histos->Add(fHIpz);
427  histos->Add(fHType);
428  histos->Add(fHCent);
429  histos->Add(fHB);
430  histos->Add(fHPhiR);
431  histos->Add(fHTime);
432  histos->Add(fHTrig);
433  fList->Add(histos);
434 
435  TList* estimators = new TList;
436  estimators->SetName("estimators");
437  estimators->SetOwner(true);
438  fList->Add(estimators);
439 
440 
441  TIter next(fCentEstimators);
442  FastCentEstimator* estimator = 0;
443  while ((estimator = static_cast<FastCentEstimator*>(next()))) {
444  estimator->Setup(estimators, fTree,sNN,fIsTgtA,fIsProjA);
445  estimator->SetVerbose(fVerbose);
446  // estimator->Print("nah");
447  }
448 
449  if (fVerbose) Info("SetupOutput", "Adding list ot outputs");
450  fOutput->Add(fList);
451  // fOutput->ls();
452 
453  return true;
454  }
460  virtual void SetupSeed()
461  {
462  std::ifstream in("/dev/urandom");
463  UInt_t seed = 0;
464  in.read(reinterpret_cast<char*>(&seed), sizeof(seed));
465  in.close();
466  Printf("=== Random =============================\n"
467  " Random number seed: %d\n"
468  "========================================", seed);
469  gRandom->SetSeed(seed);
470  }
471  virtual Bool_t SetupGRP()
472  {
473  Printf(" === Setup ==============================");
474  Printf(" Run #: %6d", fRunNo);
475  Printf(" EG: %30s", fEGName.Data());
476  Printf(" B range: %5.1ffm - %5.1ffm", fBMin, fBMax);
477  Printf(" ========================================");
478  Printf("Macro path: %s", gROOT->GetMacroPath());
479 
480  // --- Check if we shoud get the GRP line ------------------------
481  if (!fGRP && fInput) {
482  fGRP = fInput->FindObject("GRP");
483  std::ofstream* pout = new std::ofstream("grp.dat");
484  if (pout) {
485  if (fVerbose)
486  Info("SetupGRP", "Writing GRP line '%s' to \"grp.dat\"",
487  fGRP->GetTitle());
488  std::ostream& out = *pout;
489  out << fGRP->GetTitle() << std::endl;
490  pout->close();
491  }
492  }
493  if (fVerbose)
494  Info("SetupGRP", "Overrides: %p Input: %p", fOverrides, fInput);
495  if (!fOverrides && fInput) {
496  fOverrides = static_cast<TList*>(fInput->FindObject("overrides"));
497  if (!fOverrides && fVerbose)
498  Info("SetupGRP", "No GRP overrides found in input:");
499  }
500 
501  // --- Load our settings -----------------------------------------
502  if (fVerbose) Info("SetupGRP", "Loading scripts");
503  // Check if we have the global "grp" already
504  if (gROOT->ProcessLine("grp") == 0)
505  gROOT->Macro(Form("GRP.C(%d)", fRunNo));
506  if (fVerbose) Info("SetupGRP", "Perhaps override");
507  OverrideGRP();
508  gROOT->ProcessLine("grp->Print()");
509  return true;
510  }
516  virtual Bool_t SetupGen()
517  {
518  if (fVerbose) Info("SetupGen", "Load base config");
519  gROOT->Macro("BaseConfig.C");
520  if (fVerbose) Info("SetupGen", "Load EG config");
521  gROOT->Macro("EGConfig.C");
522 
523  gROOT->ProcessLine(Form("VirtualEGCfg::LoadGen(\"%s\")",fEGName.Data()));
524 
525  // --- Make our generator ----------------------------------------
526  // Info("SetupGen", "Creating generator");
527  TString egMk = Form("egCfg->MakeGenerator(\"%s\",%f,%f)",
528  fEGName.Data(), fBMin, fBMax);
529  Long64_t egPtr = gROOT->ProcessLine(egMk);
530  if (egPtr == 0) {
531  Error("SetupGen", "Failed to make generator");
532  return false;
533  }
534  fGenerator = reinterpret_cast<AliGenerator*>(egPtr);
535  TString tgt, proj;
536  Int_t tgtA=0, tgtZ=0, projA=0, projZ=0;
537  fGenerator->GetTarget(tgt, tgtA, tgtZ);
538  fGenerator->GetProjectile(proj, projA, projZ);
539  fIsTgtA = !(tgtA == tgtZ && tgtA == 1);
540  fIsProjA = !(projA == projZ && projZ == 1);
541  if (fVerbose)
542  Info("SetupGen", "tgt=%s (%3d,%2d) proj=%s (%3d,%2d) CMS=%fGeV",
543  tgt.Data(), tgtA, tgtZ, proj.Data(), projA, projZ,
544  fGenerator->GetEnergyCMS());
545  Info("SetupGen", "Generator: %s", fGenerator->GetTitle());
546  if (fFileName.IsNull()) FileName();
547 
548  return true;
549  }
555  virtual Bool_t SetupRun()
556  {
557  // --- gAlice (bare ROOT) ----------------------------------------
558  if (!gAlice)
559  new AliRun("gAlice", "The ALICE Off-line framework");
560 
561  Long64_t nev = (fNEvents <= 0 ? 0xFFFFFFFF : fNEvents);
562  Printf("=== Run ================================\n"
563  " Number of events: %lld\n"
564  "========================================", nev);
565  TObject* ord = (fInput ? fInput->FindObject("PROOF_Ordinal") : 0);
566  UShort_t saveMode = 0;
567  TString post = "";
568  TString dir = "";
569  if (ord) {
570  TObject* save = fInput->FindObject("PROOF_SaveGALICE");
571  if (save && fVerbose) {
572  Info("SetupRun", "Got save option:");
573  save->Print();
574  }
575  TString optSave(save ? save->GetTitle() : "split");
576  optSave.ToLower();
577  if (optSave.EqualTo("none")) saveMode = 0;
578  else if (optSave.EqualTo("merge")) saveMode = 1;
579  else if (optSave.EqualTo("split")) saveMode = 2;
580  if (fProofFile && saveMode > 0)
581  dir = fProofFile->GetDir(true);
582  if (saveMode > 1)
583  post = Form("_%s", ord->GetTitle());
584  }
585  TString galiceName(Form("%sgalice.root",dir.Data()));
586  TString kineName(Form("%sKinematics.root",dir.Data()));
587 
588  // --- Run-loader, stack, etc -----------------------------------
589  // Info("SetupRun", "Set-up run Loader");
590  fRunLoader = AliRunLoader::Open(galiceName, "FASTRUN", "RECREATE");
591  fRunLoader->SetKineFileName(kineName);
592  fRunLoader->SetCompressionLevel(2);
593  fRunLoader->SetNumberOfEventsPerFile(nev);
594  fRunLoader->LoadKinematics("RECREATE");
595  fRunLoader->MakeTree("E");
596  gAlice->SetRunLoader(fRunLoader);
597  fRunLoader->MakeStack();
598  fStack = fRunLoader->Stack();
599  fHeader = fRunLoader->GetHeader();
600 
601  // --- Initialize generator --------------------------------------
602  // Info("SetupRun", "Initializing generator");
603  fGenerator->Init();
604  fGenerator->SetStack(fStack);
605 
606  if (saveMode < 1) {
607  if (ord)
608  Info("SetupRun", "Not saving galice.root and Kinematics.root");
609  return true;
610  }
611 
612  TString aliceOut = Form("galice%s.root", post.Data());
613  fAliceFile = new TProofOutputFile(aliceOut, "M");
614 
615  TString kineOut = Form("Kinematics%s.root", post.Data());
616  fKineFile = new TProofOutputFile(kineOut, "M");
617 
618  return true;
619  }
624  Bool_t ReadGRPLine()
625  {
626  std::ifstream* pin = new std::ifstream("grp.dat");
627  if (!pin) {
628  Warning("ReadGRPLine", "Failed to open \"grp.dat\"");
629  return false;
630  }
631  std::istream& in = *pin;
632  TString line;
633  TString env;
634  do {
635  line.ReadLine(in);
636  if (line.IsNull()) continue;
637  if (line.BeginsWith("#")) continue;
638  env = line;
639  break;
640  } while (!in.eof());
641  pin->close();
642 
643  if (env.IsNull()) {
644  Warning("ReadGRPLine", "Got no line from \"grp.dat\"");
645  return false;
646  }
647 
648  fGRP = new TNamed("GRP",env.Data());
649  if (fVerbose) Info("ReadGRPLine", "Read \"%s\"", env.Data());
650  return true;
651  }
656  void OverrideGRP()
657  {
658  Long_t ret = gROOT->ProcessLine("grp");
659  if (ret == 0) {
660  Warning("OverrideGRP", "GRP not set yet, cannot override");
661  return;
662  }
663  if (!fOverrides) {
664  if (fVerbose) Info("OverrideGRP", "No overrides defined");
665  return;
666  }
667  TIter next(fOverrides);
668  TObject* o = 0;
669  while ((o = next())) {
670  if (fVerbose)
671  Info("OverrideGRP", "Overriding GRP setting %s with %s",
672  o->GetName(), o->GetTitle());
673  gROOT->ProcessLine(Form("grp->%s = %s;",
674  o->GetName(), o->GetTitle()));
675  }
676  Info("OverrideGRP", "After overriding:");
677  gROOT->ProcessLine("grp->Print()");
678  }
685  void AddOverride(const TString& field, const TString& value)
686  {
687  if (!fOverrides) {
688  fOverrides = new TList;
689  fOverrides->SetName("overrides");
690  }
691  fOverrides->Add(new TNamed(field, value));
692  }
697  virtual void Init(TTree*)
698  {
699  }
704  virtual void Begin(TTree*)
705  {
706  // Make a monitor
707  // Info("Begin", "gProof=%p Nomonitor=%p",
708  // gProof, (gProof ? gProof->GetParameter("NOMONITOR") : 0));
709  if (fVerbose) Info("Begin", "Called for FastSim");
710 
711  if (fMonitor > 0 && !gROOT->IsBatch()) {
712  FastSimMonitor* m = new FastSimMonitor;
713  m->Connect(fMonitor);
714  }
715  gROOT->Macro(Form("GRP.C(%d)", fRunNo));
716  if (ReadGRPLine()) {
717  if(gProof) {
718  gProof->AddInput(fGRP);
719  if (fOverrides) gProof->AddInput(fOverrides);
720  }
721  }
722  if (fVerbose) Info("Begin", "Perhaps override");
723  OverrideGRP();
724  if (fVerbose) Info("Begin", "Defining centrality estimators");
725 
726  fBEstimator = new BCentEstimator;
727  fCentEstimators = new TList;
728  // fCentEstimators->Add(new V0CentEstimator(-1)); //V0C
729  // fCentEstimators->Add(new V0CentEstimator( 0)); //V0M
730  // fCentEstimators->Add(new V0CentEstimator(+1)); //V0A
731  fCentEstimators->Add(fBEstimator);
732  fCentEstimators->Add(new V0CentEstimator(-1,true)); //V0CP
733  fCentEstimators->Add(new V0CentEstimator( 0,true)); //V0MP
734  fCentEstimators->Add(new V0CentEstimator(+1,true)); //V0AP
735  // fCentEstimators->Add(new ZNCentEstimator(-1,zN|zE)); //ZNCE
736  // fCentEstimators->Add(new ZNCentEstimator(+1,zN|zE)); //ZNAE
737  fCentEstimators->Add(new ZNCentEstimator(-1,true,false)); //ZNCE
738  fCentEstimators->Add(new ZNCentEstimator(+1,true,false)); //ZNAE
739  fCentEstimators->Add(new ZNCentEstimator(-1,true,false,true)); //ZNCEP
740  fCentEstimators->Add(new ZNCentEstimator(+1,true,false,true)); //ZNAE
741  fCentEstimators->Add(new ZNCentEstimator(-1,true,true)); //ZNCS
742  fCentEstimators->Add(new ZNCentEstimator(+1,true,true)); //ZNAS
743  // fCentEstimators->Add(new ZNCentEstimator(-1,zE|zP)); //ZPCE
744  // fCentEstimators->Add(new ZNCentEstimator(+1,zE|zP)); //ZPAE
745  // fCentEstimators->Add(new ZNCentEstimator(-1,zS)); //ZPCS
746  // fCentEstimators->Add(new ZNCentEstimator(+1,zS)); //ZPAS
747  // fCentEstimators->Add(new RefMultEstimator(0.8));
748  // fCentEstimators->Add(new RefMultEstimator(0.5));
749 #if 0
750  TIter next(fCentEstimators);
751  FastCentEstimator* estimator = 0;
752  while ((estimator = static_cast<FastCentEstimator*>(next()))) {
753  estimator->Print("nah");
754  }
755 #endif
756  if (fVerbose) Info("Begin", "End of begin");
757  }
762  void SlaveBegin(TTree*)
763  {
764  if (fVerbose) Info("SlavesBegin", "Called for FastSim");
765  SetupSeed();
766  SetupGRP();
767  SetupGen();
768  SetupOutput();
769  SetupRun();
770  }
771  /* Reset internal caches etc.
772  *
773  * @param iEv Event number
774  *
775  * @return true on success
776  */
777  virtual Bool_t PreEvent(Long64_t iEv)
778  {
779  // --- Reset header ----------------------------------------------
780  fShortHead.Reset(fRunNo, iEv);
781  fParticles->Clear();
782  // --- Reset header, etc. ---------------------------------------
783  fHeader->Reset(fRunNo, iEv);
784  fRunLoader->SetEventNumber(iEv);
785  fStack->Reset();
786  fRunLoader->MakeTree("K");
787 
788  TIter next(fCentEstimators);
789  FastCentEstimator* estimator = 0;
790  while ((estimator = static_cast<FastCentEstimator*>(next())))
791  estimator->PreEvent();
792 
793  return true;
794  }
795 
801  virtual Bool_t ProcessHeader()
802  {
803  // --- Copy to short header --------------------------------------
804  fShortHead.fRunNo = fHeader->GetRun();
805  fShortHead.fEventId = fHeader->GetEvent();
806  TArrayF ip;
807  fHeader->GenEventHeader()->PrimaryVertex(ip);
808  fShortHead.fIpX = ip[0];
809  fShortHead.fIpY = ip[1];
810  fShortHead.fIpZ = ip[2];
811 
812  // --- Check header type -----------------------------------------
813  AliGenEventHeader* genHeader = fHeader->GenEventHeader();
814  AliCollisionGeometry* geometry =
815  dynamic_cast<AliCollisionGeometry*>(genHeader);
816  AliGenPythiaEventHeader* pythia =
817  dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
818  AliGenDPMjetEventHeader* dpm =
819  dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
820  AliGenGeVSimEventHeader* gev =
821  dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
822  AliGenHerwigEventHeader* herwig =
823  dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
824  AliGenCocktailEventHeader* cocktail =
825  dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
826  AliGenHijingEventHeader* hijing =
827  dynamic_cast<AliGenHijingEventHeader*>(genHeader);
828  if (cocktail) {
829  TList* headers = cocktail->GetHeaders();
830  if (!headers) Warning("", "No headers in cocktail!");
831  TIter next(headers);
832  AliGenEventHeader* header = 0;
833  AliCollisionGeometry* geom = 0;
834  while ((header = static_cast<AliGenEventHeader*>(next()))) {
835  AliCollisionGeometry* g = dynamic_cast<AliCollisionGeometry*>(header);
836  if (g) geom = g;
837  hijing = dynamic_cast<AliGenHijingEventHeader*>(header);
838  }
839  if (geom) geometry = geom;
840  }
841  if (hijing) fShortHead.fEG = FastShortHeader::kHijing;
842  if (dpm) fShortHead.fEG = FastShortHeader::kDPMJet;
843  if (pythia) fShortHead.fEG = FastShortHeader::kPythia;
844 
845  if (geometry) {
846  fShortHead.fB = geometry->ImpactParameter();
847  fShortHead.fNtgt = geometry->TargetParticipants();
848  fShortHead.fNproj = geometry->ProjectileParticipants();
849  fShortHead.fNbin = geometry->NN();
850  fShortHead.fPhiR = geometry->ReactionPlaneAngle();
851  fShortHead.fNSpecNproj = geometry->ProjSpectatorsn();
852  fShortHead.fNSpecPproj = geometry->ProjSpectatorsp();
853  fShortHead.fNSpecNtgt = geometry->TargSpectatorsn();
854  fShortHead.fNSpecPtgt = geometry->TargSpectatorsp();
855  }
856  // --- Determine diffraction flags -------------------------------
857  Bool_t sd = false;
858  Bool_t dd = false;
859  if (pythia) {
860  Int_t type = pythia->ProcessType();
861  if (type < 100) { // pythia6
862  switch (type) {
863  case 92: case 93: sd = true; break;
864  case 94: dd = true; break;
865  }
866  }
867  else {
868  switch (type) { // Pythia8
869  case 103: case 104: sd = true; break;
870  case 105: dd = true; break;
871  }
872  }
873  fShortHead.fB = pythia->GetImpactParameter();
874  fShortHead.fNtgt = 1;
875  fShortHead.fNproj = 1;
876  fShortHead.fNbin = 1;
877  }
878  if (dpm) {
879  // Try to figure out number of spectators in either side
880  UInt_t nSpecNproj = 0;
881  UInt_t nSpecNtgt = 0;
882  UInt_t nSpecPproj = 0;
883  UInt_t nSpecPtgt = 0;
884  Int_t nPart = fStack->GetNprimary();
885  for (Int_t iPart = 0; iPart < nPart; iPart++) {
886  TParticle* particle = fStack->Particle(iPart);
887  Int_t ks = particle->GetStatusCode();
888  Int_t side = 0;
889  if (ks == 13) side = -1;
890  if (ks == 14) side = +1;
891  if (side == 0) continue;
892  Int_t kf = particle->GetPdgCode();
893  if (kf == kNeutron) {
894  if (side < 0) nSpecNproj++;
895  else nSpecNtgt++;
896  }
897  else if (kf == kProton) {
898  if (side < 0) nSpecPproj++;
899  else nSpecPtgt++;
900  }
901  }
902  fShortHead.fNSpecNtgt = nSpecNtgt;
903  fShortHead.fNSpecNproj = nSpecNproj;
904  fShortHead.fNSpecPtgt = nSpecPtgt;
905  fShortHead.fNSpecPproj = nSpecPproj;
906  // fShortHead.Print();
907 
908  Int_t type = dpm->ProcessType();
909 #ifndef NO_DPMJET_TYPE
910  switch (type) {
911  case 5: case 6: sd = true;
912  case 7: dd = true;
913  }
914 #else
915  static bool first = true;
916 
917  if (first) {
918  Func_t add = gSystem->DynFindSymbol("*", "dtglcp_");
919  if (!add)
920  Warning("", "Didn't find dtglcp_");
921  else
922  _dtglcp = (DtglcpCommon*)add;
923  }
924  // The below - or rather a different implementation with some
925  // errors - was proposed by Cvetan - I don't think it's right
926  // though. See also
927  //
928  // https://cern.ch/twiki/pub/ALICE/PAPaperCentrality/normalization.pdf
929  // https://cern.ch/twiki/bin/view/ALICE/PAMCProductionStudies
930  //
931  Int_t nsd1=0, nsd2=0, ndd=0;
932  Int_t npP = dpm->ProjectileParticipants();
933  Int_t npT = dpm->TargetParticipants();
934  // Get the numbeer of single and double diffractive participants
935  dpm->GetNDiffractive(nsd1,nsd2,ndd);
936  // Check if all partipants are single/double diffractive
937  if ((ndd == 0) && ((npP == nsd1) || (npT == nsd2))) sd = true;
938  else if (ndd == (npP + npT)) dd = true;
939  Int_t ncp = dpm->NN();
940  Int_t nct = dpm->NNw();
941  Int_t nwp = dpm->NwN();
942  Int_t nwt = dpm->NwNw();
943  Int_t nwtacc = _dtglcp->nwtacc;
944  Int_t nwtsam = _dtglcp->nwtsam;
945  if (first) {
946  Printf("@ Npp sd1 Npt sd2 dd tpe Ncp Nct Nwp Nwt acc sam");
947  first = false;
948  }
949  Printf("@ %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d",
950  npP, nsd1, npT, nsd2, ndd, type, ncp, nct, nwp, nwt,
951  nwtacc, nwtsam);
952 #endif
953  }
954  if (gev) fShortHead.fPhiR = gev->GetEventPlane();
955  if (herwig) {
956  Int_t type = herwig->ProcessType();
957  switch (type) {
958  case 5: case 6: sd = true; break;
959  }
960  fShortHead.fNtgt = 1;
961  fShortHead.fNproj = 1;
962  fShortHead.fNbin = 1;
963  }
964  fShortHead.fType = (sd ? 0x1 : 0) | (dd ? 0x2 : 0);
965 
966  // --- Check centrality -----------------------------------------
967  Double_t b = fShortHead.fB;
968  Double_t c = fBEstimator->GetCentrality(b);
969  Info("ProcessHeader", "b=%f isProjA=%d isTgtA=%d cms=%f",
970  b, fIsProjA, fIsTgtA, fGenerator ? fGenerator->GetEnergyCMS() : -1);
971  if (c >= 0) fShortHead.fC = c;
972 
973  // --- Check if within vertex cut -------------------------------
974  Bool_t selected = (fShortHead.fIpZ <= fHIpz->GetXaxis()->GetXmax() &&
975  fShortHead.fIpZ >= fHIpz->GetXaxis()->GetXmin());
976 
977  // --- Only update histograms if within IPz cut ------------------
978  if (selected) {
979  fHPhiR->Fill(fShortHead.fPhiR*TMath::RadToDeg());
980  fHB->Fill(fShortHead.fB);
981  fHIpz->Fill(fShortHead.fIpZ);
982  if (dd) fHType->Fill(3);
983  if (sd) fHType->Fill(2);
984  if (!dd && !sd) fHType->Fill(1);
985  fHCent->Fill(c);
986  fHEta ->AddBinContent(0); // Count events
987  fHY ->AddBinContent(0); // Count events
988  fHTime->AddBinContent(0); // Count events
989  // fShortHead.Print();
990  }
991  TIter next(fCentEstimators);
992  FastCentEstimator* estimator = 0;
993  while ((estimator = static_cast<FastCentEstimator*>(next())))
994  estimator->ProcessHeader(fShortHead);
995 
996  return selected;
997  }
1003  virtual void CheckTrigger(Double_t eta)
1004  {
1005  if (eta < +5.1 && eta > +2.8)
1006  fShortHead.fTrigMask |= FastShortHeader::kV0A;
1007  if (eta < -1.7 && eta > -3.7)
1008  fShortHead.fTrigMask |= FastShortHeader::kV0C;
1009  if (eta < +6.3 && eta > +4.8)
1010  fShortHead.fTrigMask |= FastShortHeader::kADA;
1011  if (eta < -4.9 && eta > -7.0)
1012  fShortHead.fTrigMask |= FastShortHeader::kADC;
1013  if (TMath::Abs(eta) < 1)
1014  fShortHead.fTrigMask |= FastShortHeader::kEta1;
1015  }
1023  virtual Bool_t ProcessParticles(Bool_t selected)
1024  {
1025  Int_t nPart = fStack->GetNprimary();
1026  for (Int_t iPart = 0; iPart < nPart; iPart++) {
1027  TParticle* particle = fStack->Particle(iPart);
1028  TParticlePDG* pdg = particle->GetPDG();
1029  Bool_t primary = fStack->IsPhysicalPrimary(iPart);
1030  Bool_t weakDecay = fStack->IsSecondaryFromWeakDecay(iPart);
1031  Bool_t charged = (pdg && TMath::Abs(pdg->Charge()) > 0);
1032  if (primary) particle->SetBit(BIT(14));
1033  if (weakDecay) particle->SetBit(BIT(15));
1034  if (charged) particle->SetBit(BIT(16));
1035 
1036  new ((*fParticles)[iPart]) TParticle(*particle);
1037 
1038  TIter next(fCentEstimators);
1039  FastCentEstimator* estimator = 0;
1040  while ((estimator = static_cast<FastCentEstimator*>(next())))
1041  estimator->Process(particle);
1042 
1043  if (!selected || !charged || !primary) continue;
1044 
1045  Double_t y = particle->Y();
1046  if (y > fHY->GetXaxis()->GetXmin() &&
1047  y < fHY->GetXaxis()->GetXmax())
1048  // Avoid filling under/overflow bins
1049  fHY->Fill(y);
1050 
1051  Double_t pT = particle->Pt();
1052  if (pT < 1e-10) continue;
1053  Double_t pZ = particle->Pz();
1054  Double_t theta = TMath::ATan2(pT, pZ);
1055  Double_t eta = -TMath::Log(TMath::Tan(theta/2));
1056  CheckTrigger(eta);
1057  if (eta > fHEta->GetXaxis()->GetXmin() &&
1058  eta < fHEta->GetXaxis()->GetXmax())
1059  // Avoid filling under/overflow bins
1060  fHEta->Fill(eta);
1061  }
1062  return true;
1063  }
1068  virtual void PostEvent()
1069  {
1070  fHeader->SetNprimary(fStack->GetNprimary());
1071  fHeader->SetNtrack(fStack->GetNtrack());
1072 
1073  fHTrig->Fill(0);
1074  if (fShortHead.fTrigMask & FastShortHeader::kV0A) fHTrig->Fill(1);
1075  if (fShortHead.fTrigMask & FastShortHeader::kV0C) fHTrig->Fill(2);
1076  if (fShortHead.fTrigMask & FastShortHeader::kADA) fHTrig->Fill(3);
1077  if (fShortHead.fTrigMask & FastShortHeader::kADC) fHTrig->Fill(4);
1078  if (fShortHead.fTrigMask & FastShortHeader::kEta1) fHTrig->Fill(5);
1079 
1080  TIter next(fCentEstimators);
1081  FastCentEstimator* estimator = 0;
1082  while ((estimator = static_cast<FastCentEstimator*>(next())))
1083  estimator->PostEvent();
1084 
1085  fTree->Fill();
1086 
1087  fStack->FinishEvent();
1088  fHeader->SetStack(fStack);
1089 
1090  fRunLoader->TreeE()->Fill();
1091  fRunLoader->WriteKinematics("OVERWRITE");
1092  }
1093  virtual void Generate() { fGenerator->Generate(); }
1094 
1102  virtual Bool_t Process(Long64_t iEv)
1103  {
1104  // --- The stopwatch ---------------------------------------------
1105  TStopwatch timer;
1106  timer.Start();
1107  PreEvent(iEv);
1108  fHTime->Fill(1, timer.RealTime());
1109 
1110  // --- Generate event --------------------------------------------
1111  timer.Start();
1112  Generate();
1113  fHTime->Fill(2, timer.RealTime());
1114 
1115  // --- Process the header ----------------------------------------
1116  timer.Start();
1117  Bool_t selected = ProcessHeader();
1118  fHTime->Fill(3, timer.RealTime());
1119 
1120  // --- Loop over particles ---------------------------------------
1121  timer.Start();
1122  ProcessParticles(selected);
1123  fHTime->Fill(4, timer.RealTime());
1124 
1125  // --- Do final stuff --------------------------------------------
1126  timer.Start();
1127  PostEvent();
1128  fHTime->Fill(5, timer.RealTime());
1129 
1130  return true;
1131  }
1132  virtual void FinishRun()
1133  {
1134  fGenerator->FinishRun();
1135  fRunLoader->WriteHeader("OVERWRITE");
1136  fGenerator->Write();
1137  fRunLoader->Write();
1138 
1139  }
1144  void SlaveTerminate()
1145  {
1146  FinishRun();
1147  if (fFile) {
1148  if (fProofFile) {
1149  if (fVerbose) fProofFile->Print();
1150  fOutput->Add(fProofFile);
1151  fOutput->Add(new TH1F("filename", fFileName.Data(),1,0,1));
1152  }
1153  // Flush out tree
1154  fFile->cd();
1155  fTree->Write(0, TObject::kOverwrite);
1156  fFile->Close();
1157  fFile->Delete();
1158  fFile = 0;
1159  }
1160  if (fAliceFile) {
1161  TFile* galice = GetGAlice();
1162  if (galice) {
1163  if (fVerbose) fAliceFile->Print();
1164  fAliceFile->AdoptFile(galice);
1165  fAliceFile->SetOutputFileName(fAliceFile->GetName());
1166  fOutput->Add(fAliceFile);
1167  galice->Write();
1168  }
1169  }
1170  if (fKineFile) {
1171  TFile* kine = GetKine();
1172  if (kine) {
1173  if (fVerbose) fKineFile->Print();
1174  fKineFile->AdoptFile(kine);
1175  fKineFile->SetOutputFileName(fKineFile->GetName());
1176  fOutput->Add(fKineFile);
1177  kine->Write();
1178  }
1179  }
1180 
1181  if (fVerbose) {
1182  Info("SlaveTerminate", "Content of output list");
1183  gROOT->IncreaseDirLevel();
1184  fOutput->ls();
1185  gROOT->DecreaseDirLevel();
1186 
1187  gSystem->Exec("echo \"Content of working directory\"");
1188  gSystem->Exec("ls -l1 | sed 's/^/ /'");
1189  }
1190  }
1198  void FlushList(TCollection* c, TDirectory* dir)
1199  {
1200  dir->cd();
1201  TIter next(c);
1202  TObject* o = 0;
1203  while ((o = next())) {
1204  if (o->IsA()->InheritsFrom(TCollection::Class())) {
1205  if (fVerbose) Info("FlushList", "Got collection: %s", c->GetName());
1206  TDirectory* cur = dir->mkdir(o->GetName());
1207  FlushList(static_cast<TCollection*>(o), cur);
1208  dir->cd();
1209  continue;
1210  }
1211  o->Write();
1212  }
1213  dir->cd();
1214  }
1215 
1220  void Terminate()
1221  {
1222  if (gProof) gProof->ClearFeedback();
1223 
1224  if (!fList)
1225  fList = static_cast<TList*>(fOutput->FindObject("histograms"));
1226  if (!fList) {
1227  Error("Terminate", "No output list");
1228  return;
1229  }
1230 
1231  if (!fProofFile) {
1232  TObject* fn = fOutput->FindObject("filename");
1233  if (fn) fFileName = fn->GetTitle();
1234  fProofFile =
1235  static_cast<TProofOutputFile*>(fOutput->FindObject(FileName()));
1236  }
1237  if (fProofFile)
1238  fFile = fProofFile->OpenFile("UPDATE");
1239  if (!fFile)
1240  fFile = TFile::Open(FileName(),"UPDATE");
1241 
1242  TList* estimators = static_cast<TList*>(fList->FindObject("estimators"));
1243  TList* histos = static_cast<TList*>(fList->FindObject("histograms"));
1244  if (!histos) {
1245  Warning("Terminate", "No histogram list found in output");
1246  fList->ls();
1247  }
1248  TIter next(fCentEstimators);
1249  FastCentEstimator* estimator = 0;
1250  while ((estimator = static_cast<FastCentEstimator*>(next())))
1251  estimator->Terminate(estimators);
1252 
1253  fHEta = static_cast<TH1*>(histos->FindObject("dNdeta"));
1254  fHY = static_cast<TH1*>(histos->FindObject("dNdy"));
1255  fHIpz = static_cast<TH1*>(histos->FindObject("ipZ"));
1256  fHType = static_cast<TH1*>(histos->FindObject("type"));
1257  fHCent = static_cast<TH1*>(histos->FindObject("cent"));
1258  fHB = static_cast<TH1*>(histos->FindObject("b"));
1259  fHPhiR = static_cast<TH1*>(histos->FindObject("phiR"));
1260  fHTime = static_cast<TH1*>(histos->FindObject("timing"));
1261  fHTrig = static_cast<TH1*>(histos->FindObject("trigger"));
1262 
1263  if (!(fHEta && fHY && fHIpz && fHType && fHB && fHPhiR && fHTime)) {
1264  Warning("Terminate", "Missing histograms (%p,%p,%p,%p,%p,%p,%p)",
1265  fHEta, fHY, fHIpz, fHType, fHB, fHPhiR, fHTime);
1266  return;
1267  }
1268 
1269  Int_t nTotal = fHIpz->GetEntries();
1270  fHEta ->Scale(1./nTotal, "width");
1271  fHY ->Scale(1./nTotal, "width");
1272  fHB ->Scale(1./nTotal, "width");
1273  fHPhiR->Scale(1./nTotal, "width");
1274  fHTime->Scale(1./nTotal, "width");
1275  fHTrig->Scale(1./nTotal, "width");
1276 
1277  if (!fFile){
1278  Warning("Terminate", "No file to write to");
1279  return;
1280  }
1281 
1282  FlushList(fList, fFile); // ->Write();
1283 
1284  fTree = static_cast<TTree*>(fFile->Get("T"));
1285  if (!fTree) Warning("Terminate", "No tree");
1286 
1287  if (fVerbose) fFile->ls();
1288  fFile->Close();
1289 
1290  MoveAliceFiles();
1291  }
1297  TFile* GetGAlice()
1298  {
1299  if (!fRunLoader) return 0;
1300  TString galiceName = fRunLoader->GetFileName();
1301  TFile* file = gROOT->GetFile(galiceName);
1302  if (!file) {
1303  Warning("GetGAlice", "Didn't find galice file \"%s\"", galiceName.Data());
1304  gROOT->GetListOfFiles()->ls();
1305  return 0;
1306  }
1307  return file;
1308  }
1314  TFile* GetKine()
1315  {
1316  if (!fRunLoader) return 0;
1317  TString kineName = "Kinematics.root";
1318  TString galiceName = fRunLoader->GetFileName();
1319  TString dir = gSystem->DirName(galiceName);
1320  if (dir.EqualTo(".")) dir = "";
1321  if (!dir.IsNull() && dir[dir.Length()-1] != '/') dir.Append("/");
1322  kineName.Prepend(dir);
1323 
1324  TFile* file = gROOT->GetFile(kineName);
1325  if (!file) {
1326  Warning("GetKine", "Didn't find kinematics file \"%s\"", kineName.Data());
1327  gROOT->GetListOfFiles()->ls();
1328  return 0;
1329  }
1330  return file;
1331  }
1338  void MoveAliceFiles()
1339  {
1340  if (!fInput) return;
1341 
1342  TObject* save = fInput->FindObject("PROOF_SaveGALICE");
1343  if (!save) return;
1344 
1345  TString sMode = save->GetTitle();
1346  if (!sMode.EqualTo("split", TString::kIgnoreCase)) return;
1347 
1348  TList* lst = new TList;
1349  TSystemDirectory* dir = new TSystemDirectory(".",
1350  gSystem->WorkingDirectory());
1351  TList* files = dir->GetListOfFiles();
1352  TSystemFile* file = 0;
1353  TIter next(files);
1354  while ((file = static_cast<TSystemFile*>(next()))) {
1355  if (file->IsDirectory()) continue;
1356  TString fn(file->GetName());
1357  if (!fn.BeginsWith("galice") && !fn.BeginsWith("Kinematics"))
1358  continue;
1359 
1360  TPRegexp regex("(.*)_([^_]+)\\.root");
1361  TObjArray* matches = regex.MatchS(fn);
1362  if (matches->GetEntriesFast() < 3) {
1363  delete matches;
1364  continue;
1365  }
1366  TString ord = matches->At(2)->GetName();
1367  TString bse = matches->At(1)->GetName();
1368 
1369  if (gSystem->AccessPathName(ord,kFileExists))
1370  gSystem->MakeDirectory(ord);
1371 
1372  if (fVerbose)
1373  Info("MoveAliceFiles", "Moving %s to %s/%s.root",
1374  fn.Data(), ord.Data(), bse.Data());
1375  file->Move(Form("%s/%s.root", ord.Data(), bse.Data()));
1376 
1377  if (!bse.EqualTo("galice")) continue;
1378  TObjString* url = new TObjString(Form("file://%s/%s/%s.root?#TE",
1379  file->GetTitle(),
1380  ord.Data(),
1381  bse.Data()));
1382  if (fVerbose)
1383  Info("MoveAliceFiles", "Adding \"%s\" to file list",
1384  url->GetName());
1385  lst->Add(url);
1386  }
1387  if (lst->GetEntries() <= 0) return;
1388  if (fVerbose) lst->ls();
1389 
1390  TFile* out = TFile::Open("index.root","RECREATE");
1391  lst->Write("TE",TObject::kSingleKey);
1392  out->Write();
1393  out->Close();
1394 
1395  }
1401  Int_t Version() const { return 1; }
1402 
1407  TString fEGName; // Name of event generator
1408  Int_t fRunNo; // Run to simulate
1409  Double_t fBMin; // Least impact parameter
1410  Double_t fBMax; // Largest impact parameter
1411  TObject* fGRP;
1412  TList* fOverrides;
1413  Long64_t fNEvents; // Number of requested events
1414  Bool_t fIsTgtA;
1415  Bool_t fIsProjA;
1416  /* @} */
1421  AliGenerator* fGenerator;
1422  AliRunLoader* fRunLoader;
1423  AliStack* fStack;
1424  AliHeader* fHeader;
1425  /* @} */
1430  TTree* fTree;
1431  TClonesArray* fParticles;
1432 
1436  TList* fList;
1437  TH1* fHEta;
1438  TH1* fHY;
1439  TH1* fHIpz;
1440  TH1* fHType;
1441  TH1* fHCent;
1442  TH1* fHB;
1443  TH1* fHPhiR;
1444  TH1* fHTime;
1445  TH1* fHTrig;
1446  /* @} */
1451  BCentEstimator* fBEstimator; // Always present
1452  TList* fCentEstimators; // Centrality estimators
1453  /* @} */
1458  TProofOutputFile* fProofFile;
1459  TProofOutputFile* fAliceFile;
1460  TProofOutputFile* fKineFile;
1461  TFile* fFile;
1462  mutable TString fFileName;
1463  /* @} */
1464  Bool_t fVerbose; // Verbosity
1465  Int_t fMonitor;
1466 
1467 #ifndef __CINT__
1468  FastShortHeader fShortHead;
1469 #endif
1470 
1484  static Bool_t LocalRun(Long64_t nev,
1485  UInt_t run,
1486  const TString& gen,
1487  Double_t bMin,
1488  Double_t bMax,
1489  Int_t monitor,
1490  Bool_t verbose,
1491  const TString& overrides="")
1492 
1493  {
1494  FastSim* sim = new FastSim(gen,run,bMin,bMax,nev,monitor);
1495  SetOverrides(sim, overrides);
1496  sim->fVerbose = verbose;
1497  sim->Begin(0);
1498  sim->SlaveBegin(0);
1499 
1500  for (Long64_t i=0; i <nev; i++) {
1501  Printf("=== Event # %6lld/%6lld ==========================",
1502  i+1, nev);
1503  sim->Process(i);
1504  }
1505  sim->SlaveTerminate();
1506  sim->Terminate();
1507 
1508  return true;
1509  }
1513  static void ProofLoadLibs()
1514  {
1515  if (!gProof) return;
1516 
1517  // Remember to copy changes to RunFast.C
1518  TList clsLib;
1519  clsLib.Add(new TNamed("TVirtualMC", "libVMC"));
1520  clsLib.Add(new TNamed("TLorentzVector", "libPhysics"));
1521  clsLib.Add(new TNamed("TLinearFitter", "libMinuit"));
1522  clsLib.Add(new TNamed("TTree", "libTree"));
1523  clsLib.Add(new TNamed("TProof", "libProof"));
1524  clsLib.Add(new TNamed("TGFrame", "libGui"));
1525  clsLib.Add(new TNamed("TSAXParser", "libXMLParser"));
1526  clsLib.Add(new TNamed("AliVEvent", "libSTEERBase"));
1527  clsLib.Add(new TNamed("AliESDEvent", "libESD"));
1528  clsLib.Add(new TNamed("AliAODEvent", "libAOD"));
1529  clsLib.Add(new TNamed("AliAnalysisManager", "libANALYSIS"));
1530  clsLib.Add(new TNamed("AliCDBManager", "libCDB"));
1531  clsLib.Add(new TNamed("AliRawVEvent", "libRAWDatabase"));
1532  clsLib.Add(new TNamed("AliHit", "libSTEER"));
1533  clsLib.Add(new TNamed("AliGenMC", "libEVGEN"));
1534  clsLib.Add(new TNamed("AliFastEvent", "libFASTSIM"));
1535 
1536  TIter next(&clsLib);
1537  TObject* obj = 0;
1538  while ((obj = next())) {
1539  gProof->Exec(Form("gROOT->LoadClass(\"%s\",\"%s\");",
1540  obj->GetName(), obj->GetTitle()));
1541  }
1542  }
1560  static Bool_t ProofRun(const TUrl& url,
1561  Long64_t nev,
1562  UInt_t run,
1563  const TString& gen,
1564  Double_t bMin,
1565  Double_t bMax,
1566  Int_t monitor=-1,
1567  Bool_t verbose=false,
1568  const TString& overrides="",
1569  const TString& save="none",
1570  const char* opt="")
1571  {
1572  TProof::Reset(url.GetUrl());
1573  TProof::Open(url.GetUrl());
1574  gProof->ClearCache();
1575 
1576  TString phy = gSystem->ExpandPathName("$(ALICE_PHYSICS)");
1577  TString ali = gSystem->ExpandPathName("$(ALICE_ROOT)");
1578  // TString fwd = gSystem->ExpandPathName("$ANA_SRC");
1579  TString fwd = phy + "/PWGLF/FORWARD/analysis2";
1580 
1581  gProof->AddIncludePath(Form("%s/include", ali.Data()));
1582  gProof->AddIncludePath(Form("%s/include", phy.Data()));
1583  ProofLoadLibs();
1584  gProof->Load(Form("%s/sim/GRP.C",fwd.Data()), true);
1585  gProof->Load(Form("%s/sim/BaseConfig.C",fwd.Data()), true);
1586  gProof->Load(Form("%s/sim/EGConfig.C",fwd.Data()), true);
1587 
1588  // gROOT->ProcessLine("gProof->SetLogLevel(5);");
1589  gProof->Load(Form("%s/sim/FastShortHeader.C", fwd.Data()));
1590  gProof->Load(Form("%s/sim/FastCentEstimators.C+%s",fwd.Data(),opt));
1591  gProof->Load(Form("%s/sim/FastMonitor.C+%s",fwd.Data(),opt));
1592  gProof->Load(Form("%s/sim/FastSim.C+%s", fwd.Data(), opt),true);
1593  gProof->SetParameter("PROOF_SaveGALICE", save);
1594 
1595  FastSim* sim = new FastSim(gen,run,bMin,bMax,nev,monitor);
1596  SetOverrides(sim, overrides);
1597  sim->fVerbose = verbose;
1598  gProof->Process(sim, nev, "");
1599 
1600  return true; // status >= 0;
1601  }
1612  static Bool_t Str2KeyVal(const TString& in,
1613  TString& key,
1614  TString& val,
1615  const char sep='=')
1616  {
1617  Int_t idx = in.Index(sep);
1618  if (idx == kNPOS) return false;
1619 
1620  key = in(0,idx);
1621  val = in(idx+1, in.Length()-idx-1);
1622  return true;
1623  }
1624  static void SetOverrides(FastSim* sim, const TString& override)
1625  {
1626  if (override.IsNull()) return;
1627 
1628  const char* valid[] = { "beamEnergy", // UInt_t [GeV]
1629  "energy", // UInt_t [GeV]
1630  "period", // String
1631  "run", // UInt_t
1632  "beam1.a", // UInt_t
1633  "beam1.z", // UInt_t
1634  "beam2.a", // UInt_t
1635  "beam2.z", // UInt_t
1636  0 };
1637  TObjArray* tokens = override.Tokenize(",");
1638  TObjString* token = 0;
1639  TIter next(tokens);
1640  while ((token = static_cast<TObjString*>(next()))) {
1641  TString& str = token->String();
1642  if (str.IsNull()) continue;
1643 
1644  TString key, val;
1645  if (!Str2KeyVal(str,key,val, ':')) {
1646  Printf("Warning: FastSim::Run: incomplete override '%s'",str.Data());
1647  continue;
1648  }
1649  const char** pvalid = valid;
1650  while (*pvalid) {
1651  if (key.EqualTo(*pvalid, TString::kIgnoreCase)) {
1652  break;
1653  }
1654  pvalid++;
1655  }
1656  if (!*pvalid) {
1657  Printf("Warning: FastSim::Run: Invalid override '%s'", key.Data());
1658  continue;
1659  }
1660  // Special case for a string
1661  if (key.EqualTo("period",TString::kIgnoreCase))
1662  val = Form("\"%s\"", val.Data());
1663  sim->AddOverride(*pvalid, val);
1664  }
1665  // delete tokens;
1666  }
1698  static Bool_t Run(const char* url, const char* opt="")
1699  {
1700  Printf("Will run fast simulation with:\n\n\t%s\n\n",url);
1701  Long64_t nev = 10000;
1702  UInt_t run = 0;
1703  TString eg = "default";
1704  TString override= "";
1705  TString save = "none";
1706  Double_t bMin = 0;
1707  Double_t bMax = 20;
1708  Int_t monitor = -1;
1709  Bool_t verbose = false;
1710  TUrl u(url);
1711  TString out;
1712  TObjArray* opts = TString(u.GetOptions()).Tokenize("&");
1713  TObjString* token = 0;
1714  TIter nextToken(opts);
1715  while ((token = static_cast<TObjString*>(nextToken()))) {
1716  TString& str = token->String();
1717  if (str.IsNull()) continue;
1718 
1719  if (str.EqualTo("verbose")) { verbose = true; str = ""; }
1720 
1721  if (str.IsNull()) continue;
1722 
1723  TString key, val;
1724  if (!Str2KeyVal(str,key,val)) {
1725  if (!out.IsNull()) out.Append("&");
1726  out.Append(str);
1727  continue;
1728  }
1729 
1730  if (key.EqualTo("events")) nev = val.Atoll();
1731  else if (key.EqualTo("run")) run = val.Atoi();
1732  else if (key.EqualTo("eg")) eg = val;
1733  else if (key.EqualTo("override")) override = val;
1734  else if (key.EqualTo("save")) save = val;
1735  else if (key.EqualTo("monitor")) monitor = val.Atoi();
1736  else if (key.EqualTo("b")) {
1737  TString min, max;
1738  if (Str2KeyVal(val, min, max, '-')) {
1739  bMin = min.Atof();
1740  bMax = max.Atof();
1741  }
1742  }
1743  else {
1744  if (!out.IsNull()) out.Append("&");
1745  out.Append(str);
1746  }
1747  }
1748  opts->Delete();
1749  u.SetOptions(out);
1750  if (!u.IsValid()) {
1751  Printf("Error: FastSim::Run: URL %s is invalid", u.GetUrl());
1752  return false;
1753  }
1754 
1755  Bool_t isLocal = TString(u.GetProtocol()).EqualTo("local");
1756 
1757  Printf("Run %s for %lld events anchored at %d\n"
1758  " Impact paramter range: %5.1f-%5.1f fm\n"
1759  " Monitor frequency: %d sec\n"
1760  " Execution url: %s",
1761  eg.Data(), nev, run, bMin, bMax, monitor, u.GetUrl());
1762 
1763 
1764  TStopwatch timer;
1765  timer.Start();
1766 
1767  Bool_t ret = false;
1768  if (isLocal)
1769  ret = LocalRun(nev, run, eg, bMin, bMax, monitor, verbose, override);
1770  else
1771  ret = ProofRun(u, nev, run, eg, bMin, bMax,
1772  monitor, verbose, override, save, opt);
1773  timer.Print();
1774 
1775  return ret;
1776  }
1777 
1778  ClassDef(FastSim,3);
1779 };
1780 
1781 
1782 struct EPosSim : public FastSim
1783 {
1784  EPosSim(UInt_t run=0, Int_t monitor=0)
1785  : FastSim("epos", run, 0, 20, 100000, monitor),
1786  fInTree(0),
1787  fInNTot(0),
1788  fInB(0),
1789  fInPDG(0),
1790  fInStatus(0),
1791  fInPx(0),
1792  fInPy(0),
1793  fInPz(0),
1794  fInE(0),
1795  fInM(0),
1796  fInNcollH(0),
1797  fInNpartP(0),
1798  fInNpartT(0),
1799  fInNcoll(0),
1800  fInNSpcPN(0),
1801  fInNSpcTN(0),
1802  fInNSpcPP(0),
1803  fInNSpcTP(0),
1804  fInPhiR(0)
1805  {
1806  }
1807  Bool_t SetupBranches()
1808  {
1809  fInNTot = fInTree->GetLeaf("nPart");
1810  fInB = fInTree->GetLeaf("ImpactParameter");
1811  fInPDG = fInTree->GetLeaf("pdgid");
1812  fInStatus = fInTree->GetLeaf("status");
1813  fInPx = fInTree->GetLeaf("px");
1814  fInPy = fInTree->GetLeaf("py");
1815  fInPz = fInTree->GetLeaf("pz");
1816  fInE = fInTree->GetLeaf("E");
1817  fInM = fInTree->GetLeaf("m");
1818  // These are probably EPOS-LHC specific
1819  fInNcollH = fInTree->GetLeaf("Ncoll_hard");
1820  fInNpartP = fInTree->GetLeaf("Npart_proj");
1821  fInNpartT = fInTree->GetLeaf("Npart_targ");
1822  fInNcoll = fInTree->GetLeaf("Ncoll");
1823  fInNSpcPN = fInTree->GetLeaf("Nspec_proj_neut");
1824  fInNSpcTN = fInTree->GetLeaf("Nspec_targ_neut");
1825  fInNSpcPP = fInTree->GetLeaf("Nspec_proj_prot");
1826  fInNSpcTP = fInTree->GetLeaf("Nspec_targ_prot");
1827  fInPhiR = fInTree->GetLeaf("phiR");
1828 
1829 
1830 
1831  return (fInNTot && fInB && fInPDG && fInPx && fInPy && fInPz && fInE);
1832  }
1833  void Begin(TTree* tree)
1834  {
1835  FastSim::Begin(tree);
1836  TIter next(fCentEstimators);
1837  FastCentEstimator* estimator = 0;
1838  while ((estimator = static_cast<FastCentEstimator*>(next()))) {
1839  if (!estimator->IsA()->InheritsFrom(V0CentEstimator::Class()))
1840  continue;
1841  V0CentEstimator* v = static_cast<V0CentEstimator*>(estimator);
1842  v->Flip(false); // Flip detector acceptance of V0A/C
1843  // v->Print("nah");
1844  }
1845  }
1846  void Init(TTree* tree)
1847  {
1848  Info("Init", "Initializing with tree %p (%s)",
1849  tree, (tree ? tree->ClassName() : ""));
1850  if (!tree) return;
1851 
1852  TFile* file = tree->GetCurrentFile();
1853  Info("Init", "Current file: (%p) %s", file,
1854  (file ? file->GetName() : ""));
1855 
1856  fInTree = tree;
1857  if (!SetupBranches())
1858  Fatal("Init", "Failed to set-up branches");
1859  // if (!SetupEstimator())
1860  // Fatal("Init", "Failed to set-up estimator");
1861  }
1867  Bool_t Notify()
1868  {
1869  if (!fInTree) {
1870  Warning("Notify", "No tree set yet!");
1871  return false;
1872  }
1873  TFile* file = fInTree->GetCurrentFile();
1874  Info("Notify", "processing file: (%p) %s", file,
1875  (file ? file->GetName() : ""));
1876  if (!file) return true;
1877  if (!SetupBranches()) {
1878  Warning("Notify", "Failed to set-up branches");
1879  return false;
1880  }
1881  return true;
1882  }
1883  const char* GetEGTitle() const { return "EPOS-LHC"; }
1884  void SetupSeed() {}
1885  Bool_t SetupGen()
1886  {
1887  fIsTgtA = gROOT->ProcessLine("grp->beam1.IsA()");
1888  fIsProjA = gROOT->ProcessLine("grp->beam2.IsA()");
1889  return true;
1890 
1891  }
1892  Bool_t SetupRun() { return true; }
1893  Bool_t PreEvent(Long64_t iEv)
1894  {
1895  Int_t read = fInTree->GetTree()->GetEntry(iEv);
1896  if (read <= 0) return false;
1897 
1898  // --- Reset header ----------------------------------------------
1899  fShortHead.Reset(fRunNo, iEv);
1900  fParticles->Clear();
1901  // Reset input
1902 
1903  TIter next(fCentEstimators);
1904  FastCentEstimator* estimator = 0;
1905  while ((estimator = static_cast<FastCentEstimator*>(next())))
1906  estimator->PreEvent();
1907 
1908  return true;
1909  }
1910  Bool_t ProcessHeader()
1911  {
1912  fShortHead.fIpX = 0;
1913  fShortHead.fIpY = 0;
1914  fShortHead.fIpZ = 0;
1915  fShortHead.fB = fInB->GetValue();
1916  fShortHead.fNtgt = (fInNpartT ? fInNpartT->GetValue() : 0);
1917  fShortHead.fNproj = (fInNpartP ? fInNpartP->GetValue() : 0);
1918  fShortHead.fNbin = (fInNcoll ? fInNcoll ->GetValue() : 0);
1919  fShortHead.fPhiR = (fInPhiR ? fInPhiR ->GetValue() : 0);
1920  fShortHead.fNSpecNproj = (fInNSpcPN ? fInNSpcPN->GetValue() : 0);
1921  fShortHead.fNSpecNtgt = (fInNSpcTN ? fInNSpcTN->GetValue() : 0);
1922  fShortHead.fNSpecPproj = (fInNSpcPP ? fInNSpcPP->GetValue() : 0);
1923  fShortHead.fNSpecPtgt = (fInNSpcTP ? fInNSpcTP->GetValue() : 0);
1924  fShortHead.fEG = FastShortHeader::kEPOS;
1925 
1926  Double_t c = fBEstimator->GetCentrality(fShortHead.fB);
1927  if (c >= 0) fShortHead.fC = c;
1928 
1929  // --- Check if within vertex cut -------------------------------
1930  Bool_t selected = (fShortHead.fIpZ <= fHIpz->GetXaxis()->GetXmax() &&
1931  fShortHead.fIpZ >= fHIpz->GetXaxis()->GetXmin());
1932 
1933  // --- Only update histograms if within IPz cut ------------------
1934  if (selected) {
1935  fHPhiR->Fill(fShortHead.fPhiR*TMath::RadToDeg());
1936  fHB->Fill(fShortHead.fB);
1937  fHIpz->Fill(fShortHead.fIpZ);
1938  fHCent->Fill(c);
1939  // fShortHead.Print();
1940  }
1941  TIter next(fCentEstimators);
1942  FastCentEstimator* estimator = 0;
1943  while ((estimator = static_cast<FastCentEstimator*>(next())))
1944  estimator->ProcessHeader(fShortHead);
1945  return selected;
1946  }
1947  virtual Bool_t ProcessParticles(Bool_t selected)
1948  {
1949  Int_t nTot = fInNTot->GetValue();
1950  for (Int_t iPart = 0; iPart < nTot; iPart++) {
1951  Int_t status = fInStatus->GetValue(iPart);
1952  Int_t pdg = fInPDG->GetValue(iPart);
1953  Double_t px = fInPx->GetValue(iPart);
1954  Double_t py = fInPy->GetValue(iPart);
1955  Double_t pz = fInPz->GetValue(iPart);
1956  // Double_t pz = -fInPz->GetValue(iPart); // flip sign on pZ
1957  Double_t e = fInE->GetValue(iPart);
1958  // Double_t m = fInM->GetValue(iPart);
1959 
1960  TParticle* particle =
1961  new ((*fParticles)[iPart]) TParticle(pdg, status,-1,-1,-1,-1,
1962  px, py, pz, e, 0, 0, 0, 0);
1963  TParticlePDG* pdgP = particle->GetPDG();
1964  Bool_t primary = status == 1;
1965  Bool_t weakDecay = false;
1966  Bool_t charged = (pdgP && TMath::Abs(pdgP->Charge()) > 0);
1967  if (primary) particle->SetBit(BIT(14));
1968  if (weakDecay) particle->SetBit(BIT(15));
1969  if (charged) particle->SetBit(BIT(16));
1970 
1971  TIter next(fCentEstimators);
1972  FastCentEstimator* estimator = 0;
1973  while ((estimator = static_cast<FastCentEstimator*>(next())))
1974  estimator->Process(particle);
1975 
1976  if (!selected || !charged || !primary) continue;
1977  fHY ->Fill(particle->Y());
1978  Double_t pT = particle->Pt();
1979  if (pT < 1e-10) continue;
1980  Double_t pZ = particle->Pz();
1981  Double_t theta = TMath::ATan2(pT, pZ);
1982  Double_t eta = -TMath::Log(TMath::Tan(theta/2));
1983  CheckTrigger(eta);
1984  fHEta->Fill(eta);
1985  }
1986  return true;
1987  }
1988  void PostEvent()
1989  {
1990  fTree->Fill();
1991 
1992  TIter next(fCentEstimators);
1993  FastCentEstimator* estimator = 0;
1994  while ((estimator = static_cast<FastCentEstimator*>(next())))
1995  estimator->PostEvent();
1996  }
1997  void Generate() {}
1998  void FinishRun() {}
1999  TTree* fInTree;
2000  TLeaf* fInNTot;
2001  TLeaf* fInB;
2002  TLeaf* fInPDG;
2003  TLeaf* fInStatus;
2004  TLeaf* fInPx;
2005  TLeaf* fInPy;
2006  TLeaf* fInPz;
2007  TLeaf* fInE;
2008  TLeaf* fInM;
2009  TLeaf* fInNcollH;
2010  TLeaf* fInNpartP;
2011  TLeaf* fInNpartT;
2012  TLeaf* fInNcoll;
2013  TLeaf* fInNSpcPN;
2014  TLeaf* fInNSpcTN;
2015  TLeaf* fInNSpcPP;
2016  TLeaf* fInNSpcTP;
2017  TLeaf* fInPhiR;
2018 
2029  static Bool_t LocalRun(Long64_t nev,
2030  UInt_t run,
2031  Int_t monitor,
2032  Bool_t verbose)
2033 
2034  {
2035  EPosSim* sim = new EPosSim(run, monitor);
2036  sim->fVerbose = verbose;
2037  sim->Begin(0);
2038  sim->SlaveBegin(0);
2039 
2040  for (Long64_t i=0; i <nev; i++) {
2041  Printf("=== Event # %6lld/%6lld ==========================",
2042  i+1, nev);
2043  sim->Process(i);
2044  }
2045  sim->SlaveTerminate();
2046  sim->Terminate();
2047 
2048  return true;
2049  }
2058  static Bool_t SetupProof(const TUrl& url,
2059  const char* opt="")
2060  {
2061  TProof::Reset(url.GetUrl());
2062  TProof::Open(url.GetUrl());
2063  gProof->ClearCache();
2064 
2065  TString phy = gSystem->ExpandPathName("$(ALICE_PHYSICS)");
2066  TString ali = gSystem->ExpandPathName("$(ALICE_ROOT)");
2067  // TString fwd = gSystem->ExpandPathName("$ANA_SRC");
2068  TString fwd = phy + "/PWGLF/FORWARD/analysis2";
2069 
2070  gProof->AddIncludePath(Form("%s/include", ali.Data()));
2071  gProof->AddIncludePath(Form("%s/include", phy.Data()));
2072  ProofLoadLibs();
2073  gProof->Load(Form("%s/sim/GRP.C",fwd.Data()), true);
2074  gProof->Load(Form("%s/sim/BaseConfig.C",fwd.Data()), true);
2075  gProof->Load(Form("%s/sim/EGConfig.C",fwd.Data()), true);
2076 
2077  // gROOT->ProcessLine("gProof->SetLogLevel(5);");
2078  gProof->Load(Form("%s/sim/FastMonitor.C+%s",fwd.Data(),opt));
2079  gProof->Load(Form("%s/sim/FastShortHeader.C", fwd.Data()));
2080  gProof->Load(Form("%s/sim/FastCentEstimators.C+%s",fwd.Data(),opt));
2081  gProof->Load(Form("%s/sim/FastSim.C+%s", fwd.Data(), opt),true);
2082 
2083  return true; // status >= 0;
2084  }
2116  static Bool_t Run(const char* url, const char* opt="")
2117  {
2118  Printf("Will run fast simulation with:\n\n\t%s\n\n",url);
2119  UInt_t run = 0;
2120  Long64_t nev = -1;
2121  Int_t monitor = -1;
2122  Bool_t verbose = false;
2123  TUrl u(url);
2124  TString out;
2125  TObjArray* opts = TString(u.GetOptions()).Tokenize("&");
2126  TObjString* token = 0;
2127  TIter nextToken(opts);
2128  while ((token = static_cast<TObjString*>(nextToken()))) {
2129  TString& str = token->String();
2130  if (str.IsNull()) continue;
2131 
2132  if (str.EqualTo("verbose")) { verbose = true; str = ""; }
2133 
2134  if (str.IsNull()) continue;
2135 
2136  TString key, val;
2137  if (!Str2KeyVal(str,key,val)) {
2138  if (!out.IsNull()) out.Append("&");
2139  out.Append(str);
2140  continue;
2141  }
2142 
2143  if (key.EqualTo("run")) run = val.Atoi();
2144  else if (key.EqualTo("monitor")) monitor = val.Atoi();
2145  else if (key.EqualTo("events")) nev = val.Atoi();
2146  else {
2147  if (!out.IsNull()) out.Append("&");
2148  out.Append(str);
2149  }
2150  }
2151  opts->Delete();
2152  u.SetOptions(out);
2153  if (!u.IsValid()) {
2154  Printf("Error: FastSim::Run: URL %s is invalid", u.GetUrl());
2155  return false;
2156  }
2157 
2158  Printf("Run EPos anchored at %d\n"
2159  " Monitor frequency: %d sec\n"
2160  " Execution url: %s",
2161  run, monitor, u.GetUrl());
2162 
2163  TString treeName = u.GetAnchor();
2164  if (treeName.IsNull()) treeName = "Particle";
2165  TFile* file = TFile::Open(u.GetFile(), "READ");
2166  if (!file) {
2167  Printf("Error: FastAnalysis::Run: Failed to open %s",
2168  u.GetFile());
2169  return false;
2170  }
2171 
2172  TChain* chain = new TChain(treeName, treeName);
2173  TObject* o = file->Get(treeName);
2174  if (!o) {
2175  Printf("Error: FastAnalysis::Run: Couldn't get %s from %s",
2176  treeName.Data(), u.GetFile());
2177  file->Close();
2178  return false;
2179  }
2180  Int_t cret = 0;
2181  if (o->IsA()->InheritsFrom(TChain::Class()))
2182  cret = chain->Add(static_cast<TChain*>(o));
2183  else if (o->IsA()->InheritsFrom(TTree::Class()))
2184  cret = chain->AddFile(u.GetFile());
2185  else if (o->IsA()->InheritsFrom(TCollection::Class()))
2186  cret = chain->AddFileInfoList(static_cast<TCollection*>(o));
2187  else if (o->IsA()->InheritsFrom(TFileCollection::Class()))
2188  cret = chain->AddFileInfoList(static_cast<TFileCollection*>(o)
2189  ->GetList());
2190  file->Close();
2191  if (cret <= 0 || chain->GetListOfFiles()->GetEntries() <= 0) {
2192  Printf("Error: FastAnalysis::Run: Failed to create chain");
2193  return false;
2194  }
2195 
2196  TString proto = u.GetProtocol();
2197  Bool_t isProof = (proto.EqualTo("proof") || proto.EqualTo("lite"));
2198  if (isProof) {
2199  if (!SetupProof(u,opt)) return false;
2200  chain->SetProof();
2201  }
2202 
2203  EPosSim* sim = new EPosSim(run, monitor);
2204  sim->fVerbose = verbose;
2205  if (nev < 0) nev = TChain::kBigNumber;
2206 
2207  TStopwatch timer;
2208  timer.Start();
2209 
2210  Long64_t ret = chain->Process(sim, "", nev, 0);
2211  timer.Print();
2212 
2213  return ret > 0;
2214  }
2215  Int_t Version() const { return 2; }
2216  ClassDef(EPosSim,1);
2217 };
2218 
2219 #endif
2220 //
2221 // EOF
2222 //
EPosSim::PreEvent
Bool_t PreEvent(Long64_t iEv)
Definition: FastSim.C:1893
FastSim::ProofLoadLibs
static void ProofLoadLibs()
Definition: FastSim.C:1513
FastShortHeader::Reset
void Reset(UInt_t runNo, UInt_t eventNo)
Definition: FastShortHeader.C:61
FastSim::SlaveBegin
void SlaveBegin(TTree *)
Definition: FastSim.C:762
pdg
Int_t pdg
Definition: charmCutsOptimization.C:34
FastSim::fBMin
Double_t fBMin
Definition: FastSim.C:1409
FastCentEstimators.C
FastSim::SetupGRP
virtual Bool_t SetupGRP()
Definition: FastSim.C:471
FastSim::fIsProjA
Bool_t fIsProjA
True if target beam is nuclei.
Definition: FastSim.C:1415
EPosSim::GetEGTitle
const char * GetEGTitle() const
Definition: FastSim.C:1883
FastSim::fHTrig
TH1 * fHTrig
Timing.
Definition: FastSim.C:1445
EPosSim::Run
static Bool_t Run(const char *url, const char *opt="")
Definition: FastSim.C:2116
Double_t
double Double_t
Definition: External.C:58
FastSim::Init
virtual void Init(TTree *)
Definition: FastSim.C:697
FastSim::Process
virtual Bool_t Process(Long64_t iEv)
Definition: FastSim.C:1102
EPosSim::fInB
TLeaf * fInB
Definition: FastSim.C:2001
url
const char * url
Definition: ExtractNoiseGain.C:9
EPosSim::fInNSpcPN
TLeaf * fInNSpcPN
Definition: FastSim.C:2013
DtglcpCommon::nct
Int_t nct
Definition: FastSim.C:87
FastSim::Begin
virtual void Begin(TTree *)
Definition: FastSim.C:704
FastSim::FlushList
void FlushList(TCollection *c, TDirectory *dir)
Definition: FastSim.C:1198
Long64_t
long long Long64_t
Definition: External.C:43
FastSim::GetGAlice
TFile * GetGAlice()
Definition: FastSim.C:1297
EPosSim::SetupProof
static Bool_t SetupProof(const TUrl &url, const char *opt="")
Definition: FastSim.C:2058
FastSim::MoveAliceFiles
void MoveAliceFiles()
Definition: FastSim.C:1338
EPosSim::fInNcoll
TLeaf * fInNcoll
Definition: FastSim.C:2012
DtglcpCommon::nwbsam
Int_t nwbsam
Definition: FastSim.C:82
DtglcpCommon::ncp
Int_t ncp
Definition: FastSim.C:86
gSystem
TSystem * gSystem
FastCentEstimator::PostEvent
virtual void PostEvent()
Definition: FastCentEstimators.C:84
FastSim::fBMax
Double_t fBMax
Definition: FastSim.C:1410
FastSim::fStack
AliStack * fStack
Loader of trees.
Definition: FastSim.C:1423
EPosSim::ProcessHeader
Bool_t ProcessHeader()
Definition: FastSim.C:1910
EPosSim::LocalRun
static Bool_t LocalRun(Long64_t nev, UInt_t run, Int_t monitor, Bool_t verbose)
Definition: FastSim.C:2029
FastCentEstimator::Terminate
virtual void Terminate(TCollection *out)=0
FastSim::fGenerator
AliGenerator * fGenerator
Definition: FastSim.C:1421
DtglcpCommon::rproj
Double_t rproj
Definition: FastSim.C:77
EPosSim::fInNTot
TLeaf * fInNTot
Definition: FastSim.C:2000
c
TCanvas * c
Definition: TestFitELoss.C:172
FastSim::CheckTrigger
virtual void CheckTrigger(Double_t eta)
Definition: FastSim.C:1003
FastSim::fAliceFile
TProofOutputFile * fAliceFile
Proof output file.
Definition: FastSim.C:1459
FastSim::fOverrides
TList * fOverrides
GRP in one line.
Definition: FastSim.C:1412
EPosSim::Init
void Init(TTree *tree)
Definition: FastSim.C:1846
EPosSim::FinishRun
void FinishRun()
Definition: FastSim.C:1998
FastCentEstimator::ProcessHeader
virtual void ProcessHeader(FastShortHeader &)
Definition: FastCentEstimators.C:72
EPosSim::fInTree
TTree * fInTree
Definition: FastSim.C:1999
FastSim::GetTitle
const char * GetTitle() const
Definition: FastSim.C:218
FastSim::fProofFile
TProofOutputFile * fProofFile
Definition: FastSim.C:1458
FastSim::fTree
TTree * fTree
Definition: FastSim.C:1430
DtglcpCommon::nwtbac
Int_t nwtbac
Definition: FastSim.C:85
kNeutron
Definition: FastCentEstimators.C:24
DtglcpCommon::nwtacc
Int_t nwtacc
Definition: FastSim.C:83
EPosSim::Notify
Bool_t Notify()
Definition: FastSim.C:1867
gRandom
TRandom * gRandom
kProton
Definition: FastCentEstimators.C:25
EPosSim::SetupBranches
Bool_t SetupBranches()
Definition: FastSim.C:1807
FastSim::OverrideGRP
void OverrideGRP()
Definition: FastSim.C:656
TUrl
Definition: External.C:92
FastSim::fVerbose
Bool_t fVerbose
Output file name.
Definition: FastSim.C:1464
FastMonitor::Connect
void Connect(Int_t freq=-1)
Definition: FastMonitor.C:136
EPosSim::fInPhiR
TLeaf * fInPhiR
Definition: FastSim.C:2017
EPosSim::fInPy
TLeaf * fInPy
Definition: FastSim.C:2005
FastCentEstimator::Print
virtual void Print(Option_t *option="") const
Definition: FastCentEstimators.C:91
FastSim::FileName
const char * FileName() const
Definition: FastSim.C:166
FastSim::Str2KeyVal
static Bool_t Str2KeyVal(const TString &in, TString &key, TString &val, const char sep='=')
Definition: FastSim.C:1612
FastSim::Version
Int_t Version() const
Definition: FastSim.C:1401
FastSim::fHB
TH1 * fHB
Event type histogram.
Definition: FastSim.C:1442
FastSim::fIsTgtA
Bool_t fIsTgtA
Definition: FastSim.C:1414
EPosSim::Begin
void Begin(TTree *tree)
Definition: FastSim.C:1833
FastSim::fMonitor
Int_t fMonitor
Definition: FastSim.C:1465
EPosSim::fInNcollH
TLeaf * fInNcollH
Definition: FastSim.C:2009
EPosSim::fInPz
TLeaf * fInPz
Definition: FastSim.C:2006
FastSim::FinishRun
virtual void FinishRun()
Definition: FastSim.C:1132
FastSim::GetKine
TFile * GetKine()
Definition: FastSim.C:1314
Int_t
int Int_t
Definition: External.C:63
FastSim::fBEstimator
BCentEstimator * fBEstimator
Definition: FastSim.C:1451
EPosSim::SetupGen
Bool_t SetupGen()
Definition: FastSim.C:1885
TH1I
Definition: External.C:204
UInt_t
unsigned int UInt_t
Definition: External.C:33
EPosSim::fInNpartP
TLeaf * fInNpartP
Definition: FastSim.C:2010
FastSim::fGRP
TObject * fGRP
Definition: FastSim.C:1411
Float_t
float Float_t
Definition: External.C:68
FastSim::fCentEstimators
TList * fCentEstimators
Definition: FastSim.C:1452
V0CentEstimator::Flip
void Flip(Bool_t onlySign=true)
Definition: FastCentEstimators.C:568
FastSim::fShortHead
FastShortHeader fShortHead
Definition: FastSim.C:1468
_dtglcp
DtglcpCommon * _dtglcp
Definition: FastSim.C:89
FastShortHeader.C
EPosSim::Generate
void Generate()
Definition: FastSim.C:1997
FastSim::fHIpz
TH1 * fHIpz
dN/dy
Definition: FastSim.C:1439
FastSim::fHPhiR
TH1 * fHPhiR
B histogram.
Definition: FastSim.C:1443
BCentEstimator::GetCentrality
Double_t GetCentrality(Double_t b) const
Definition: FastCentEstimators.C:476
EPosSim::SetupRun
Bool_t SetupRun()
Definition: FastSim.C:1892
EPosSim::fInNpartT
TLeaf * fInNpartT
Definition: FastSim.C:2011
EPosSim::SetupSeed
void SetupSeed()
Definition: FastSim.C:1884
TH1D
Definition: External.C:212
DtglcpCommon::nwtaac
Int_t nwtaac
Definition: FastSim.C:84
FastSim::GetName
const char * GetName() const
Definition: FastSim.C:212
FastSim::GetSNN
virtual UShort_t GetSNN() const
Definition: FastSim.C:224
FastSim::FastSim
FastSim(const char *eg="", ULong_t runNo=0, Double_t bMin=0, Double_t bMax=20, Long64_t nEvents=0, Int_t monitor=0)
Definition: FastSim.C:124
FastSim::fHType
TH1 * fHType
IPz histogram.
Definition: FastSim.C:1440
FastSim::SetupSeed
virtual void SetupSeed()
Definition: FastSim.C:460
ULong_t
unsigned long ULong_t
Definition: External.C:38
DtglcpCommon::rtarg
Double_t rtarg
Definition: FastSim.C:78
FastSim::fKineFile
TProofOutputFile * fKineFile
Definition: FastSim.C:1460
FastSim::fHeader
AliHeader * fHeader
Stack of particles.
Definition: FastSim.C:1424
FastSim::fHEta
TH1 * fHEta
List of outputs.
Definition: FastSim.C:1437
FastSim::ReadGRPLine
Bool_t ReadGRPLine()
Definition: FastSim.C:624
FastSim::Generate
virtual void Generate()
Definition: FastSim.C:1093
FastSim::fRunNo
Int_t fRunNo
Definition: FastSim.C:1408
EPosSim::fInPDG
TLeaf * fInPDG
Definition: FastSim.C:2002
FastSim::fList
TList * fList
Definition: FastSim.C:1436
FastSim::Run
static Bool_t Run(const char *url, const char *opt="")
Definition: FastSim.C:1698
FastSim::SetupGen
virtual Bool_t SetupGen()
Definition: FastSim.C:516
FastSim::fHTime
TH1 * fHTime
Reaction plane.
Definition: FastSim.C:1444
EPosSim::Version
Int_t Version() const
Definition: FastSim.C:2215
FastSim::fFile
TFile * fFile
Definition: FastSim.C:1461
EPosSim::fInStatus
TLeaf * fInStatus
Definition: FastSim.C:2003
EPosSim::fInM
TLeaf * fInM
Definition: FastSim.C:2008
fwd
const char * fwd
Definition: ExtractNoiseGain.C:7
FastCentEstimator::SetVerbose
void SetVerbose(Bool_t verb)
Definition: FastCentEstimators.C:54
DtglcpCommon::nwasam
Int_t nwasam
Definition: FastSim.C:81
GetListOfFiles
TList * GetListOfFiles(const char *input=".")
Definition: RunFinalQA.C:114
FastCentEstimator::Process
virtual void Process(const TParticle *p)=0
FastMonitor::Register
void Register(TObject *descr, Bool_t proof=true)
Definition: FastMonitor.C:186
FastShortHeader::fEG
enum FastShortHeader::EG_t fEG
FastSim::Terminate
void Terminate()
Definition: FastSim.C:1220
FastSim::PreEvent
virtual Bool_t PreEvent(Long64_t iEv)
Definition: FastSim.C:777
EPosSim::fInNSpcTN
TLeaf * fInNSpcTN
Definition: FastSim.C:2014
FastMonitor.C
EPosSim::fInNSpcTP
TLeaf * fInNSpcTP
Definition: FastSim.C:2016
nEvents
Float_t nEvents[nProd]
Input train file.
Definition: DrawProductionComparison.C:95
DtglcpCommon::bimpac
Double_t bimpac
Definition: FastSim.C:79
EPosSim::EPosSim
EPosSim(UInt_t run=0, Int_t monitor=0)
Definition: FastSim.C:1784
FastSim::ProofRun
static Bool_t ProofRun(const TUrl &url, Long64_t nev, UInt_t run, const TString &gen, Double_t bMin, Double_t bMax, Int_t monitor=-1, Bool_t verbose=false, const TString &overrides="", const TString &save="none", const char *opt="")
Definition: FastSim.C:1560
EPosSim::fInPx
TLeaf * fInPx
Definition: FastSim.C:2004
FastSimMonitor::FastSimMonitor
FastSimMonitor(TSelector *s=0)
Definition: FastSim.C:94
FastSim::fRunLoader
AliRunLoader * fRunLoader
Event generator.
Definition: FastSim.C:1422
EPosSim::PostEvent
void PostEvent()
Definition: FastSim.C:1988
file
TFile * file
TList with histograms for a given trigger.
Definition: DrawAnaCaloTrackQA.C:77
FastSim::PostEvent
virtual void PostEvent()
Definition: FastSim.C:1068
FastSim::fHY
TH1 * fHY
dN/deta
Definition: FastSim.C:1438
FastSim::fFileName
TString fFileName
Output file.
Definition: FastSim.C:1462
FastSim::SetupOutput
Bool_t SetupOutput()
Definition: FastSim.C:256
FastSim::LocalRun
static Bool_t LocalRun(Long64_t nev, UInt_t run, const TString &gen, Double_t bMin, Double_t bMax, Int_t monitor, Bool_t verbose, const TString &overrides="")
Definition: FastSim.C:1484
UShort_t
unsigned short UShort_t
Definition: External.C:28
FastSim::fParticles
TClonesArray * fParticles
Custom tree.
Definition: FastSim.C:1431
DtglcpCommon::nwtsam
Int_t nwtsam
Definition: FastSim.C:80
FastSim::SetupRun
virtual Bool_t SetupRun()
Definition: FastSim.C:555
Bool_t
bool Bool_t
Definition: External.C:53
FastSim::fEGName
TString fEGName
Definition: FastSim.C:1407
FastCentEstimator::PreEvent
virtual void PreEvent()
Definition: FastCentEstimators.C:71
FastSim::AddOverride
void AddOverride(const TString &field, const TString &value)
Definition: FastSim.C:685
FastSim::GetEGTitle
virtual const char * GetEGTitle() const
Definition: FastSim.C:245
FastSim::ProcessHeader
virtual Bool_t ProcessHeader()
Definition: FastSim.C:801
FastSim::ProcessParticles
virtual Bool_t ProcessParticles(Bool_t selected)
Definition: FastSim.C:1023
EPosSim::ProcessParticles
virtual Bool_t ProcessParticles(Bool_t selected)
Definition: FastSim.C:1947
FastSim::SlaveTerminate
void SlaveTerminate()
Definition: FastSim.C:1144
FastCentEstimator::Setup
virtual void Setup(TCollection *out, TTree *tree, UShort_t sNN, Bool_t tgtA, Bool_t projA)=0
FastSim::fHCent
TH1 * fHCent
Event type histogram.
Definition: FastSim.C:1441
TH1
Definition: External.C:196
EPosSim::fInNSpcPP
TLeaf * fInNSpcPP
Definition: FastSim.C:2015
FastSim::fNEvents
Long64_t fNEvents
GRP setting to override.
Definition: FastSim.C:1413
dir
TDirectoryFile * dir
Definition: DrawAnaCaloTrackQA.C:75
FastSim::SetOverrides
static void SetOverrides(FastSim *sim, const TString &override)
Definition: FastSim.C:1624
EPosSim::fInE
TLeaf * fInE
Definition: FastSim.C:2007