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