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