6 # include "AliGenerator.h"
7 # include "AliRunLoader.h"
9 # include "AliHeader.h"
10 # include "AliGenEventHeader.h"
12 # include "AliCollisionGeometry.h"
13 # include "AliGenPythiaEventHeader.h"
14 # include "AliGenDPMjetEventHeader.h"
15 # include "AliGenGeVSimEventHeader.h"
16 # include "AliGenHerwigEventHeader.h"
20 # include <TParticle.h>
23 # include <TClonesArray.h>
26 # include <TParticlePDG.h>
27 # include <TStopwatch.h>
29 # include <TProofOutputFile.h>
38 class AliGenEventHeader;
45 class TProofOutputFile;
51 //====================================================================
53 * Monitor output objects
55 struct FastMonitor : public TObject, public TQObject
63 FastMonitor(TSelector* s=0)
64 : fName("FastMonitor"),
68 if (gROOT->IsBatch()) {
69 Warning("FastMonitor", "Batch processing, no monitoring");
74 fName = gProof->GetSessionTag();
75 gDirectory->Add(this);
76 Bool_t ret = gProof->Connect("Feedback(TList *objs)", "FastMonitor", this,
77 "Feedback(TList *objs)");
79 Warning("FastMonitor", "Failed to connect to Proof");
85 fCanvas = new TCanvas(fName, Form("Monitor %s", fName.Data()), 1000, 800);
86 fCanvas->SetFillColor(0);
87 fCanvas->SetFillStyle(0);
88 fCanvas->SetTopMargin(0.01);
89 fCanvas->SetRightMargin(0.01);
92 RegisterDraw(1, "type", "", 0);
93 RegisterDraw(2, "b", "", 0);
94 RegisterDraw(3, "cent", "", 0);
95 RegisterDraw(4, "dNdeta", "", 0x8);
98 * Register a draw of a an object
100 * @param i Pad number
101 * @param name Name of object
102 * @param option Drawing option
108 * - 0x8 Scale to events and bin width
110 void RegisterDraw(Int_t i,
115 TVirtualPad* p = fCanvas->GetPad(i);
117 Warning("RegisterDraw", "Not enough sub-pads (%d)", i);
122 p->SetTopMargin(0.01);
123 p->SetRightMargin(0.01);
124 p->SetName(Form("p_%s", name));
126 if (flags & 0x1) p->SetLogx();
127 if (flags & 0x2) p->SetLogy();
128 if (flags & 0x4) p->SetLogz();
129 if (flags & 0x8) p->SetBit(BIT(15));
134 virtual ~FastMonitor()
137 gProof->Disconnect("Feedback(TList *objs)",this,
138 "Feedback(TList* objs)");
141 * Set name of this object
145 void SetName(const char* name) { fName = name; }
147 * Get the name of this object
151 const char* GetName() const { return fName.Data(); }
153 * Find pad corresponding to an object
155 * @param name Name of object
157 * @return Pointer to pad or null
159 TVirtualPad* FindPad(const TString& name)
163 TString t = Form("p_%s", name.Data());
164 while ((p = fCanvas->GetPad(i))) {
165 if (t.EqualTo(p->GetName())) return p;
171 * Called when we get notified of
173 * @param objs List of monitored objects
175 void Feedback(TList* objs)
177 // Info("FeedBack", "List is %p", objs);
178 // if (objs) objs->ls();
179 if (!fCanvas) return;
181 TList* l = static_cast<TList*>(objs->FindObject("histograms"));
183 Warning("Feedback", "No histograms");
187 TObject* oIpz = l->FindObject("ipZ");
188 if (oIpz && oIpz->IsA()->InheritsFrom(TH1::Class()))
189 nEvents = static_cast<TH1*>(oIpz)->GetEntries();
191 Warning("Feedback", "Histogram ipZ not found");
195 while ((o = next())) {
196 TVirtualPad* p = FindPad(o->GetName());
198 // Info("FeedBack", "no pad for %s", o->GetName());
202 if (o->IsA()->InheritsFrom(TH1::Class())) {
203 TH1* h = static_cast<TH1*>(o);
204 TH1* c = h->DrawCopy(p->GetTitle());
206 c->SetBit(TObject::kCanDelete);
207 if (p->TestBit(BIT(15))) {
208 Info("Feedback", "Scaling %s by 1./%d and width",
209 c->GetName(), nEvents);
210 c->Scale(1./nEvents, "width");
214 TObject* c = o->DrawClone(p->GetTitle());
215 c->SetBit(TObject::kCanDelete);
224 * Function to handle connect signals
232 * Function to handle timer events
234 Bool_t HandleTimer(TTimer*)
236 Info("HandleTimer", "Selector=%p", fSelector);
237 if (!fSelector) return false;
238 Feedback(fSelector->GetOutputList());
245 /** Possibly link to selector */
246 TSelector* fSelector;
247 ClassDef(FastMonitor,1);
251 //====================================================================
253 * Run a event generator simulation
255 struct FastSim : public TSelector
260 * @param eg Event generator
261 * @param runNo Run number to simulate
262 * @param bMin Lease impact parameter
263 * @param bMax Largest impact parameter
265 FastSim(const char* eg="",
295 const char* FileName() const
299 if (!fFileName.IsNull()) fn = fFileName;
301 const char* egName = (fGenerator ?
302 fGenerator->GetName() :
304 fn = Form("%s_%09d", egName, fRunNo);
306 if (fNEvents >= 1000000)
307 fn.Append(Form("_%lldM", fNEvents/1000000));
308 else if (fNEvents >= 1000)
309 fn.Append(Form("_%lldk", fNEvents/1000));
311 fn.Append(Form("_%lld", fNEvents));
319 if (fFileName.IsNull())
320 fFileName = Form("%s_%09d.root", fEGName.Data(), fRunNo);
321 return fFileName.Data();*/
323 const char* GetName() const { return "FastSim"; }
324 const char* GetTitle() const { return "ALICE Event Generator simulation"; }
329 * @return true on success
333 Info("SetupOutput", "First the file");
334 Bool_t isProof = false;
335 if (fInput && fInput->FindObject("PROOF_Ordinal"))
338 Info("SetupOutput", "Making Proof File");
339 fProofFile = new TProofOutputFile(FileName(), "M");
340 // TProofOutputFile::kMerge,
341 // TProofOutputFile::kRemote);
342 fFile = fProofFile->OpenFile("RECREATE");
345 fFile = TFile::Open(FileName(), "RECREATE");
347 Info("SetupOutput", "Making our tree");
348 fTree = new TTree("T", "T");
349 fParticles = new TClonesArray("TParticle");
350 fTree->Branch("header", &fShortHead,
351 "run/i:event:npart:nbin:type:ipx/D:ipy:ipz:b:c:phir");
352 fTree->Branch("particles", &fParticles);
354 fTree->SetDirectory(fFile);
355 fTree->SetAlias("primary", "(particles.fBits&(1<<14))");
356 fTree->SetAlias("weak", "(particles.fBits&(1<<15))");
357 fTree->SetAlias("charged", "(particles.fBits&(1<<16))");
358 fTree->SetAlias("pt", "(sqrt(pow(particles.fPx,2)+"
359 /* */"pow(particles.fPy,2)))");
360 fTree->SetAlias("eta", "(pt<1e-10?1024:"
361 "-log(tan(atan2(particles.Pt(),particles.fPz)/2)))");
362 fTree->SetAlias("good", "(primary&&charged&&abs(eta)<1000)");
363 fTree->SetAlias("sd", "(header.fType & 0x1)");
364 fTree->SetAlias("dd", "(header.fType & 0x2)");
365 fTree->SetAlias("pion", "(abs(particles.fPdgCode)==211)");
366 fTree->SetAlias("kaon", "(abs(particles.fPdgCode)==321)");
367 fTree->SetAlias("proton", "(abs(particles.fPdgCode)==2212)");
368 fTree->SetAlias("electron","(abs(particles.fPdgCode)==11)");
369 fTree->SetAlias("other", "(!pion&&!kaon&&!proton&&!electron)");
370 fTree->SetAlias("beta", "(particles.P()/particle.Energy())");
371 fTree->SetAlias("gamma", "(1./sqrt(1-beta*beta))");
373 Info("SetupOutput", "Making histograms");
374 Double_t maxEta = 10;
375 Double_t dEta = 10./200;
376 fHEta = new TH1D("dNdeta", "Charged particle pseudo-rapidity density",
377 Int_t(2*maxEta/dEta+.5), -maxEta, +maxEta);
379 fHEta->SetXTitle("#it{#eta}");
380 fHEta->SetYTitle("1/N d#it{N}_{ch}/d#it{#eta}");
381 fHEta->SetMarkerColor(kRed+2);
382 fHEta->SetMarkerStyle(20);
383 fHEta->SetDirectory(0);
385 fHIpz = new TH1D("ipZ", "Z-coordinate of interaction point",
387 fHIpz->SetMarkerColor(kGreen+2);
388 fHIpz->SetFillColor(kGreen+2);
389 fHIpz->SetFillStyle(3001);
390 fHIpz->SetXTitle("IP_{#it{z}} [cm]");
391 fHIpz->SetDirectory(0);
393 fHType = new TH1D("type", "Diffractive", 3, .5, 3.5);
394 fHType->SetMarkerColor(kOrange+2);
395 fHType->SetFillColor(kOrange+2);
396 fHType->SetFillStyle(3001);
397 fHType->SetDirectory(0);
398 fHType->GetXaxis()->SetBinLabel(1, "Non");
399 fHType->GetXaxis()->SetBinLabel(2, "Single");
400 fHType->GetXaxis()->SetBinLabel(3, "Double");
402 fHCent = new TH1D("cent", "Centrality", 20, 0, 100);
403 fHCent->SetMarkerColor(kPink+2);
404 fHCent->SetFillColor(kPink+2);
405 fHCent->SetFillStyle(3001);
406 fHCent->SetDirectory(0);
407 fHCent->SetXTitle("Centrality [%]");
409 fHB = new TH1D("b", "Impact parameter", 20, 0, 20);
410 fHB->SetMarkerColor(kYellow+2);
411 fHB->SetFillColor(kYellow+2);
412 fHB->SetFillStyle(3001);
413 fHB->SetXTitle("#it{b} [fm]");
414 fHB->SetDirectory(0);
416 fHPhiR = new TH1D("phiR", "Event plane angle", 360, 0, 360);
417 fHPhiR->SetMarkerColor(kMagenta+2);
418 fHPhiR->SetFillColor(kMagenta+2);
419 fHPhiR->SetFillStyle(3001);
420 fHPhiR->SetXTitle("#it{#Phi}_{R} [degrees]");
421 fHPhiR->SetDirectory(0);
423 fHTime = new TH1D("timing", "Timing of processing", 5,0.5,5.5);
424 fHTime->SetMarkerColor(kBlue+2);
425 fHTime->SetFillColor(kBlue+2);
426 fHTime->SetFillStyle(3001);
427 fHTime->GetXaxis()->SetBinLabel(1,"Reset");
428 fHTime->GetXaxis()->SetBinLabel(2,"Generate");
429 fHTime->GetXaxis()->SetBinLabel(3,"Header");
430 fHTime->GetXaxis()->SetBinLabel(4,"Particles");
431 fHTime->GetXaxis()->SetBinLabel(5,"Filling");
432 fHTime->SetDirectory(0);
435 fList->SetName("histograms");
436 fList->SetOwner(true);
445 Info("SetupOutput", "Adding list ot outputs");
452 Printf(" === Setup ==============================");
453 Printf(" Run #: %6d", fRunNo);
454 Printf(" EG: %30s", fEGName.Data());
455 Printf(" B range: %5.1ffm - %5.1ffm", fBMin, fBMax);
456 Printf(" ========================================");
457 Printf("Macro path: %s", gROOT->GetMacroPath());
459 // --- Check if we shoud get the GRP line ------------------------
460 if (!fGRP && fInput) {
461 fGRP = fInput->FindObject("GRP");
462 std::ofstream* pout = new std::ofstream("grp.dat");
464 Info("SetupGen", "Writing GRP line '%s' to \"grp.dat\"",
466 std::ostream& out = *pout;
467 out << fGRP->GetTitle() << std::endl;
472 // --- Load our settings -----------------------------------------
473 Info("SetupGen", "Loading scripts");
474 gROOT->Macro(Form("GRP.C(%d)", fRunNo));
475 gROOT->Macro("BaseConfig.C");
476 gROOT->Macro("EGConfig.C");
478 gROOT->ProcessLine(Form("VirtualEGCfg::LoadGen(\"%s\")",fEGName.Data()));
480 // --- Make our generator ----------------------------------------
481 Info("SetupGen", "Creating generator");
482 TString egMk = Form("egCfg->MakeGenerator(\"%s\",%f,%f)",
483 fEGName.Data(), fBMin, fBMax);
484 Long64_t egPtr = gROOT->ProcessLine(egMk);
486 Error("Setup", "Failed to make generator");
489 fGenerator = reinterpret_cast<AliGenerator*>(egPtr);
492 if (fFileName.IsNull()) FileName();
493 Info("SetupRun", "File name is '%s'", fFileName.Data());
498 * Setup the generator etc. of the job
500 * @param nev Maximum number of events per file
502 * @return true on success
506 // --- gAlice (bare ROOT) ----------------------------------------
508 new AliRun("gAlice", "The ALICE Off-line framework");
510 Long64_t nev = (fNEvents <= 0 ? 0xFFFFFFFF : fNEvents);
511 // --- Run-loader, stack, etc -----------------------------------
512 Info("SetupRun", "Set-up run Loader");
513 fRunLoader = AliRunLoader::Open("galice.root", "FASTRUN", "RECREATE");
514 fRunLoader->SetCompressionLevel(2);
515 fRunLoader->SetNumberOfEventsPerFile(nev);
516 fRunLoader->LoadKinematics("RECREATE");
517 fRunLoader->MakeTree("E");
518 gAlice->SetRunLoader(fRunLoader);
519 fRunLoader->MakeStack();
520 fStack = fRunLoader->Stack();
521 fHeader = fRunLoader->GetHeader();
523 // --- Initialize generator --------------------------------------
524 Info("SetupRun", "Initializing generator");
526 fGenerator->SetStack(fStack);
531 * Read the previously created grp.dat file
536 std::ifstream* pin = new std::ifstream("grp.dat");
538 Warning("ReadGRPLine", "Failed to open \"grp.dat\"");
541 std::istream& in = *pin;
546 if (line.IsNull()) continue;
547 if (line.BeginsWith("#")) continue;
554 Warning("ReadGRPLine", "Got no line from \"grp.dat\"");
558 fGRP = new TNamed("GRP",env.Data());
576 Info("Begin", "gProof=%p Nomonitor=%p",
577 gProof, (gProof ? gProof->GetParameter("NOMONITOR") : 0));
579 if (gProof && !gProof->GetParameter("NOMONITOR")) {
581 gProof->AddFeedback("histograms");
582 Info("Begin", "Adding monitoring");
584 gROOT->Macro(Form("GRP.C(%d)", fRunNo));
587 gProof->AddInput(fGRP);
592 * Set-up this sub-job
595 void SlaveBegin(TTree*)
601 /* Reset internal caches etc.
603 * @param iEv Event number
605 * @return true on success
607 Bool_t PreEvent(Long64_t iEv)
609 // --- Reset header ----------------------------------------------
610 fShortHead.fRunNo = fRunNo;
611 fShortHead.fEventId = iEv;
612 fShortHead.fIpX = 1024;
613 fShortHead.fIpY = 1024;
614 fShortHead.fIpZ = 1024;
615 fShortHead.fNpart = -1;
616 fShortHead.fNbin = -1;
617 fShortHead.fPhiR = -1;
621 // --- Reset header, etc. ---------------------------------------
622 fHeader->Reset(fRunNo, iEv);
623 fRunLoader->SetEventNumber(iEv);
625 fRunLoader->MakeTree("K");
630 * Process the event header
632 * @return true if the event should be diagnosed
634 Bool_t ProcessHeader()
636 // --- Copy to short header --------------------------------------
637 fShortHead.fRunNo = fHeader->GetRun();
638 fShortHead.fEventId = fHeader->GetEvent();
640 fHeader->GenEventHeader()->PrimaryVertex(ip);
641 fShortHead.fIpX = ip[0];
642 fShortHead.fIpY = ip[1];
643 fShortHead.fIpZ = ip[2];
645 // --- Check header type -----------------------------------------
646 AliGenEventHeader* genHeader = fHeader->GenEventHeader();
647 AliCollisionGeometry* geometry =
648 dynamic_cast<AliCollisionGeometry*>(genHeader);
649 AliGenPythiaEventHeader* pythia =
650 dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
651 AliGenDPMjetEventHeader* dpm =
652 dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
653 AliGenGeVSimEventHeader* gev =
654 dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
655 AliGenHerwigEventHeader* herwig =
656 dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
658 fShortHead.fB = geometry->ImpactParameter();
659 fShortHead.fNpart = (geometry->ProjectileParticipants() +
660 geometry->TargetParticipants());
661 fShortHead.fNbin = geometry->NN();
662 fShortHead.fPhiR = geometry->ReactionPlaneAngle();
664 // --- Determine diffraction flags -------------------------------
668 Int_t type = pythia->ProcessType();
669 if (type < 100) { // pythia6
671 case 92: case 93: sd = true; break;
672 case 94: dd = true; break;
676 switch (type) { // Pythia8
677 case 103: case 104: sd = true; break;
678 case 105: dd = true; break;
681 fShortHead.fB = pythia->GetImpactParameter();
682 fShortHead.fNpart = 2;
683 fShortHead.fNbin = 1;
686 Int_t type = dpm->ProcessType();
688 case 5: case 6: sd = true;
693 if (gev) fShortHead.fPhiR = gev->GetEventPlane();
695 Int_t type = herwig->ProcessType();
697 case 5: case 6: sd = true; break;
699 fShortHead.fNpart = 2;
700 fShortHead.fNbin = 1;
702 fShortHead.fType = (sd ? 0x1 : 0) | (dd ? 0x2 : 0);
704 // --- Check centrality -----------------------------------------
706 // Updated 4th of November 2014 from
707 // cern.ch/twiki/bin/view/ALICE/CentStudies#Tables_with_centrality_bins_AN1
711 Double_t b = fShortHead.fB;
713 if (0.00 >= b && b < 1.57) { c=0.5; np=403.8; nc=1861; }
714 else if (1.57 >= b && b < 2.22) { c=1.5; np=393.6; nc=1766; }
715 else if (2.22 >= b && b < 2.71) { c=2.5; np=382.9; nc=1678; }
716 else if (2.71 >= b && b < 3.13) { c=3.5; np=372; nc=1597; }
717 else if (3.13 >= b && b < 3.50) { c=4.5; np=361.1; nc=1520; }
718 else if (3.50 >= b && b < 4.94) { c=7.5; np=329.4; nc=1316; }
719 else if (4.94 >= b && b < 6.05) { c=12.5; np=281.2; nc=1032; }
720 else if (6.05 >= b && b < 6.98) { c=17.5; np=239; nc=809.8; }
721 else if (6.98 >= b && b < 7.81) { c=22.5; np=202.1; nc=629.6; }
722 else if (7.81 >= b && b < 8.55) { c=27.5; np=169.5; nc=483.7; }
723 else if (8.55 >= b && b < 9.23) { c=32.5; np=141; nc=366.7; }
724 else if (9.23 >= b && b < 9.88) { c=37.5; np=116; nc=273.4; }
725 else if (9.88 >= b && b < 10.47) { c=42.5; np=94.11; nc=199.4; }
726 else if (10.47 >= b && b < 11.04) { c=47.5; np=75.3; nc=143.1; }
727 else if (11.04 >= b && b < 11.58) { c=52.5; np=59.24; nc=100.1; }
728 else if (11.58 >= b && b < 12.09) { c=57.5; np=45.58; nc=68.46; }
729 else if (12.09 >= b && b < 12.58) { c=62.5; np=34.33; nc=45.79; }
730 else if (12.58 >= b && b < 13.05) { c=67.5; np=25.21; nc=29.92; }
731 else if (13.05 >= b && b < 13.52) { c=72.5; np=17.96; nc=19.08; }
732 else if (13.52 >= b && b < 13.97) { c=77.5; np=12.58; nc=12.07; }
733 else if (13.97 >= b && b < 14.43) { c=82.5; np=8.812; nc=7.682; }
734 else if (14.43 >= b && b < 14.96) { c=87.5; np=6.158; nc=4.904; }
735 else if (14.96 >= b && b < 15.67) { c=92.5; np=4.376; nc=3.181; }
736 else if (15.67 >= b && b < 20.00) { c=97.5; np=3.064; nc=1.994; }
738 // Be careful to round off
739 if (fShortHead.fNpart <= 0) fShortHead.fNpart = Int_t(np+.5);
740 if (fShortHead.fNbin <= 0) fShortHead.fNbin = Int_t(nc+.5)/2;
743 // --- Check if within vertex cut -------------------------------
744 Bool_t selected = (fShortHead.fIpZ <= fHIpz->GetXaxis()->GetXmax() &&
745 fShortHead.fIpZ >= fHIpz->GetXaxis()->GetXmin());
747 // --- Only update histograms if within IPz cut ------------------
749 fHPhiR->Fill(fShortHead.fPhiR*TMath::RadToDeg());
750 fHB->Fill(fShortHead.fB);
751 fHIpz->Fill(fShortHead.fIpZ);
752 fHType->Fill(dd ? 3 : sd ? 2 : 1);
758 * Process all particles
760 * @param selected True if particle information should be diagnosed
762 * @return true on success
764 Bool_t ProcessParticles(Bool_t selected)
766 Int_t nPart = fStack->GetNprimary();
767 for (Int_t iPart = 0; iPart < nPart; iPart++) {
768 TParticle* particle = fStack->Particle(iPart);
769 TParticlePDG* pdg = particle->GetPDG();
770 Bool_t primary = fStack->IsPhysicalPrimary(iPart);
771 Bool_t weakDecay = fStack->IsSecondaryFromWeakDecay(iPart);
772 Bool_t charged = (pdg && TMath::Abs(pdg->Charge()) > 0);
773 if (primary) particle->SetBit(BIT(14));
774 if (weakDecay) particle->SetBit(BIT(15));
775 if (charged) particle->SetBit(BIT(16));
777 new ((*fParticles)[iPart]) TParticle(*particle);
779 if (!selected || !charged || !primary) continue;
780 Double_t pT = particle->Pt();
781 if (pT < 1e-10) continue; /// Along beam axis
782 Double_t pZ = particle->Pz();
783 Double_t theta = TMath::ATan2(pT, pZ);
784 Double_t eta = -TMath::Log(TMath::Tan(theta/2));
790 * Do final event processing (fill output)
795 fHeader->SetNprimary(fStack->GetNprimary());
796 fHeader->SetNtrack(fStack->GetNtrack());
800 fStack->FinishEvent();
801 fHeader->SetStack(fStack);
803 fRunLoader->TreeE()->Fill();
804 fRunLoader->WriteKinematics("OVERWRITE");
809 * @param iEv Event number
811 * @return true on success, false otherwize
813 Bool_t Process(Long64_t iEv)
815 // --- The stopwatch ---------------------------------------------
819 fHTime->Fill(1, timer.RealTime());
821 // --- Generate event --------------------------------------------
823 fGenerator->Generate();
824 fHTime->Fill(2, timer.RealTime());
826 // --- Process the header ----------------------------------------
828 Bool_t selected = ProcessHeader();
829 fHTime->Fill(3, timer.RealTime());
831 // --- Loop over particles ---------------------------------------
833 ProcessParticles(selected);
834 fHTime->Fill(4, timer.RealTime());
836 // --- Do final stuff --------------------------------------------
839 fHTime->Fill(5, timer.RealTime());
844 * Finalize this sub-job
847 void SlaveTerminate()
849 fGenerator->FinishRun();
850 fRunLoader->WriteHeader("OVERWRITE");
856 fOutput->Add(fProofFile);
857 fOutput->Add(new TH1F("filename", fFileName.Data(),1,0,1));
861 fTree->Write(0, TObject::kOverwrite);
868 * Final processing of the data
873 if (gProof) gProof->ClearFeedback();
876 fList = static_cast<TList*>(fOutput->FindObject("histograms"));
878 Error("Terminate", "No output list");
883 TObject* fn = fOutput->FindObject("filename");
884 if (fn) fFileName = fn->GetTitle();
886 static_cast<TProofOutputFile*>(fOutput->FindObject(FileName()));
889 fFile = fProofFile->OpenFile("UPDATE");
891 fFile = TFile::Open(FileName(),"UPDATE");
894 fHEta = static_cast<TH1*>(fList->FindObject("dNdeta"));
895 fHIpz = static_cast<TH1*>(fList->FindObject("ipZ"));
896 fHType = static_cast<TH1*>(fList->FindObject("type"));
897 fHCent = static_cast<TH1*>(fList->FindObject("cent"));
898 fHB = static_cast<TH1*>(fList->FindObject("b"));
899 fHPhiR = static_cast<TH1*>(fList->FindObject("phiR"));
900 fHTime = static_cast<TH1*>(fList->FindObject("timing"));
902 if (!(fHEta && fHIpz && fHType && fHB && fHPhiR && fHTime)) {
903 Warning("Terminate", "Missing histograms (%p,%p,%p,%p,%p,%p)",
904 fHEta, fHIpz, fHType, fHB, fHPhiR, fHTime);
908 Int_t nTotal = fHIpz->GetEntries();
909 fHEta ->Scale(1./nTotal, "width");
910 fHB ->Scale(1./nTotal, "width");
911 fHPhiR->Scale(1./nTotal, "width");
912 fHTime->Scale(1./nTotal, "width");
915 Warning("Terminate", "No file to write to");
927 fTree = static_cast<TTree*>(fFile->Get("T"));
928 if (!fTree) Warning("Terminate", "No tree");
933 * Interface version used
937 Int_t Version() const { return 1; }
943 TString fEGName; // Name of event generator
944 Int_t fRunNo; // Run to simulate
945 Double_t fBMin; // Least impact parameter
946 Double_t fBMax; // Largest impact parameter
947 TObject* fGRP; //! GRP in one line
948 Long64_t fNEvents; // Number of requested events
952 * @name ALICE EG interface
954 AliGenerator* fGenerator; //! Event generator
955 AliRunLoader* fRunLoader; //! Loader of trees
956 AliStack* fStack; //! Stack of particles
957 AliHeader* fHeader; //! Header handler
961 * @name Custom output
963 TTree* fTree; //! Custom tree
964 TClonesArray* fParticles; //! List of particles
969 TList* fList; //! List of outputs
970 TH1* fHEta; //! dN/deta
971 TH1* fHIpz; //! IPz histogram
972 TH1* fHType; //! Event type histogram
973 TH1* fHCent; //! Event type histogram
974 TH1* fHB; //! B histogram
975 TH1* fHPhiR; //! Reaction plane
976 TH1* fHTime; //! Timing
982 TProofOutputFile* fProofFile; //! Proof output file
983 TFile* fFile; //! Output file
984 mutable TString fFileName; //! Output file name
1004 * Run this selector as a normal process
1006 * @param nev Number of events
1007 * @param run Run number to anchor in
1008 * @param gen Generator
1009 * @param bMin Least impact parameter [fm]
1010 * @param bMax Largest impact parameter [fm]
1011 * @param monitor Monitor frequency [s]
1013 * @return true on succes
1015 static Bool_t Run(Long64_t nev,
1022 TStopwatch stopwatch;
1025 FastSim* sim = new FastSim(gen,run,bMin,bMax,nev);
1031 // timer = new TTimer(new FastMonitor(sim), monitor*1000,true);
1032 timer = new TTimer(1000);
1033 timer->Connect("Timeout()","FastMonitor",
1034 new FastMonitor(sim), "Handle()");
1035 ::Info("Run", "Turning on monitoring");
1036 timer->Start(-1,false);
1039 for (Long64_t i=0; i <nev; i++) {
1040 Printf("=== Event # %6lld/%6lld ==========================",
1043 if (timer && (i > 0) && (i % 500 == 0)) {
1044 if (timer->CheckTimer(gSystem->Now()))
1045 Printf("Fired timer");
1048 if (timer) timer->TurnOff();
1049 sim->SlaveTerminate();
1055 static void ProofLoadLibs()
1057 if (!gProof) return;
1059 // Remember to copy changes to RunFast.C
1061 clsLib.Add(new TNamed("TVirtualMC", "libVMC"));
1062 clsLib.Add(new TNamed("TLorentzVector", "libPhysics"));
1063 clsLib.Add(new TNamed("TLinearFitter", "libMinuit"));
1064 clsLib.Add(new TNamed("TTree", "libTree"));
1065 clsLib.Add(new TNamed("TProof", "libProof"));
1066 clsLib.Add(new TNamed("TGFrame", "libGui"));
1067 clsLib.Add(new TNamed("TSAXParser", "libXMLParser"));
1068 clsLib.Add(new TNamed("AliVEvent", "libSTEERBase"));
1069 clsLib.Add(new TNamed("AliESDEvent", "libESD"));
1070 clsLib.Add(new TNamed("AliAODEvent", "libAOD"));
1071 clsLib.Add(new TNamed("AliAnalysisManager", "libANALYSIS"));
1072 clsLib.Add(new TNamed("AliCDBManager", "libCDB"));
1073 clsLib.Add(new TNamed("AliRawVEvent", "libRAWDatabase"));
1074 clsLib.Add(new TNamed("AliHit", "libSTEER"));
1075 clsLib.Add(new TNamed("AliGenMC", "libEVGEN"));
1076 clsLib.Add(new TNamed("AliFastEvent", "libFASTSIM"));
1078 TIter next(&clsLib);
1080 while ((obj = next())) {
1081 gProof->Exec(Form("gROOT->LoadClass(\"%s\",\"%s\");",
1082 obj->GetName(), obj->GetTitle()));
1086 * Run this selector in PROOF(Lite)
1088 * @param url Proof URL
1089 * @param nev Number of events
1090 * @param run Run number to anchor in
1091 * @param gen Generator
1092 * @param bMin Least impact parameter [fm]
1093 * @param bMax Largest impact parameter [fm]
1094 * @param monitor Monitor frequency [s]
1095 * @param opt Compilation options
1097 * @return true on succes
1099 static Bool_t Proof(const char* url,
1108 Printf("# events: %lld", nev);
1109 Printf("Run #: %u", run);
1110 Printf("Generator: %s", gen);
1111 Printf("b range: %5.1f-%5.1f", bMin, bMax);
1112 Printf("monitor: %ds", monitor);
1119 gProof->ClearCache();
1121 TString ali = gSystem->ExpandPathName("$(ALICE_ROOT)");
1122 // TString fwd = gSystem->ExpandPathName("$ANA_SRC");
1123 TString fwd = ali + "/PWGLF/FORWARD/analysis2";
1125 gProof->AddIncludePath(Form("%s/include", ali.Data()));
1127 gProof->Load(Form("%s/sim/GRP.C",fwd.Data()), true);
1128 gProof->Load(Form("%s/sim/BaseConfig.C",fwd.Data()), true);
1129 gProof->Load(Form("%s/sim/EGConfig.C",fwd.Data()), true);
1131 // gROOT->ProcessLine("gProof->SetLogLevel(5);");
1132 gProof->Load(Form("%s/sim/FastSim.C+%s", fwd.Data(), opt),true);
1134 if (monitor <= 0) gProof->SetParameter("NOMONITOR", true/*ignored*/);
1135 else gProof->SetParameter("PROOF_FeedbackPeriod",
1136 monitor*1000/*ms*/);
1138 FastSim* sim = new FastSim(gen,run,bMin,bMax,nev);
1139 gProof->Process(sim, nev, "");
1142 return true; // status >= 0;
1144 ClassDef(FastSim,1);