]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG0/highMultiplicity/AliHighMultiplicitySelector.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
adding offline triggers
[u/mrichter/AliRoot.git] / PWG0 / highMultiplicity / AliHighMultiplicitySelector.cxx
CommitLineData
9608df41 1/* $Id$ */
2
3#include "AliHighMultiplicitySelector.h"
4
5#include <TVector3.h>
6#include <TFile.h>
7#include <TH1F.h>
8#include <TH2F.h>
9#include <TNtuple.h>
10#include <TProfile.h>
11#include <TParticle.h>
12#include <TCanvas.h>
13#include <TTree.h>
14#include <TGeoManager.h>
15#include <TF1.h>
16#include <TGraph.h>
17#include <TLegend.h>
18#include <TLine.h>
dd367a14 19#include <TMath.h>
9608df41 20
21#include <AliLog.h>
22#include <AliESD.h>
23#include <AliRunLoader.h>
24#include <AliStack.h>
25
26#include <AliITSgeom.h>
27#include <AliITSLoader.h>
28#include <AliITSdigitSPD.h>
29#include <AliITSRecPoint.h>
30
31#include "AliPWG0Helper.h"
32
33//
a9aed06c 34// Selector that produces plots needed for the high multiplicity analysis with the
35// pixel detector
36// Needs cluster file
9608df41 37//
38
39ClassImp(AliHighMultiplicitySelector)
40
41AliHighMultiplicitySelector::AliHighMultiplicitySelector() :
42 AliSelectorRL(),
43 fChipsLayer1(0),
44 fChipsLayer2(0),
45 fL1vsL2(0),
46 fMvsL1(0),
47 fMvsL2(0),
48 fChipsFired(0),
49 fPrimaryL1(0),
50 fPrimaryL2(0),
51 fClusterZL1(0),
52 fClusterZL2(0),
53 fClvsL1(0),
54 fClvsL2(0),
a9aed06c 55 fCentralRegion(kFALSE)
9608df41 56{
57 //
58 // Constructor. Initialization of pointers
59 //
60}
61
62AliHighMultiplicitySelector::~AliHighMultiplicitySelector()
63{
64 //
65 // Destructor
66 //
67
68 // histograms are in the output list and deleted when the output
69 // list is deleted by the TSelector dtor
70}
71
72void AliHighMultiplicitySelector::SlaveBegin(TTree *tree)
73{
a9aed06c 74 // create histograms
75
9608df41 76 AliSelectorRL::SlaveBegin(tree);
77
78 fChipsFired = new TH2F("fChipsFired", ";Module;Chip;Count", 240, -0.5, 239.5, 5, -0.5, 4.5);
79
80 fPrimaryL1 = new TNtuple("fPrimaryL1", "", "multiplicity:firedChips:chipsByPrimaries:clusters");
81 fPrimaryL2 = new TNtuple("fPrimaryL2", "", "multiplicity:firedChips:chipsByPrimaries:clusters");
82
83 fClusterZL1 = new TH1F("fClusterZL1", ";z", 400, -20, 20);
84 fClusterZL2 = new TH1F("fClusterZL2", ";z", 400, -20, 20);
85}
86
87void AliHighMultiplicitySelector::Init(TTree* tree)
88{
89 // read the user objects
90
91 AliSelectorRL::Init(tree);
92
93 // enable only the needed branches
94 if (tree)
95 {
96 tree->SetBranchStatus("*", 0);
97 tree->SetBranchStatus("fTriggerMask", 1);
98
99 /*if (fTree->GetCurrentFile())
100 {
101 TString fileName(fTree->GetCurrentFile()->GetName());
102 fileName.ReplaceAll("AliESDs", "geometry");
103
104 // load geometry
105 TGeoManager::Import(fileName);
106 }*/
107 }
108}
109
110Bool_t AliHighMultiplicitySelector::Process(Long64_t entry)
111{
112 //
113 // processing
114 //
115
116 if (AliSelectorRL::Process(entry) == kFALSE)
117 return kFALSE;
118
119 // Check prerequisites
120 if (!fESD)
121 {
122 AliDebug(AliLog::kError, "ESD branch not available");
123 return kFALSE;
124 }
125
69b09e3b 126 Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB1);
9608df41 127
128 if (!eventTriggered)
129 {
130 AliDebug(AliLog::kDebug, "Event not triggered. Skipping.");
131 return kTRUE;
132 }
133
134 AliStack* stack = GetStack();
135 if (!stack)
136 {
137 AliDebug(AliLog::kError, "Stack not available");
138 return kFALSE;
139 }
140
141 Int_t nPrim = stack->GetNprimary();
142 Int_t multiplicity21 = 0;
143 Int_t multiplicity16 = 0;
144
145 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
146 {
147 AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
148
149 TParticle* particle = stack->Particle(iMc);
150
151 if (!particle)
152 {
153 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
154 continue;
155 }
156
157 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
158 continue;
159
a9aed06c 160 if (fCentralRegion)
9608df41 161 {
162 if (TMath::Abs(particle->Eta()) < 1.05)
163 multiplicity21++;
164 if (TMath::Abs(particle->Eta()) < 0.8)
165 multiplicity16++;
166 }
167 else
168 {
169 if (TMath::Abs(particle->Eta()) < 2.1)
170 multiplicity21++;
171 if (TMath::Abs(particle->Eta()) < 1.6)
172 multiplicity16++;
173 }
174 }// end of mc particle
175
176 AliRunLoader* runLoader = GetRunLoader();
177 if (!runLoader)
178 {
179 AliDebug(AliLog::kError, "runloader not available");
180 return kFALSE;
181 }
182
183 // TDirectory::TContext restores the current directory is restored when the scope ends.
184 // This helps around ROOT bug #26025 and is good behaviour anyway
185 TDirectory::TContext context(0);
186 AliITSLoader* loader = (AliITSLoader*) runLoader->GetLoader( "ITSLoader" );
187 loader->LoadDigits("READ");
188 TTree* treeD = loader->TreeD();
189 if (!treeD)
190 {
191 AliDebug(AliLog::kError, "Could not retrieve TreeD of ITS");
192 return kFALSE;
193 }
194
195 treeD->SetBranchStatus("*", 0);
196 treeD->SetBranchStatus("ITSDigitsSPD.fTracks*", 1);
197 treeD->SetBranchStatus("ITSDigitsSPD.fCoord1", 1);
198
199 TClonesArray* digits = 0;
200 treeD->SetBranchAddress("ITSDigitsSPD", &digits);
a6e0ebfe 201 if (digits)
9608df41 202 digits->Clear();
203
204 // each value for both layers
205 Int_t totalNumberOfFO[2];
206 Int_t chipsHitByPrimaries[2];
207 //Int_t chipsHitBySecondaries[2];
208
209 for (Int_t i=0; i<2; ++i)
210 {
211 totalNumberOfFO[i] = 0;
212 chipsHitByPrimaries[i] = 0;
213 //chipsHitBySecondaries[i] = 0;
214 }
215
216 Int_t startSPD = 0; //geom->GetStartSPD();
217 Int_t lastSPD = 239; //geom->GetLastSPD();
218
219 //printf("%d %d\n", startSPD, lastSPD);
220// for (Int_t l=0; l<240; ++l) { AliITSgeomTGeo::GetModuleId(l, i, j, k); printf("%d --> %d\n", l, i); }
221
222 // loop over modules (ladders)
223 for (Int_t moduleIndex=startSPD; moduleIndex<lastSPD+1; moduleIndex++)
224 {
a9aed06c 225 if (fCentralRegion)
9608df41 226 {
227 if ((moduleIndex % 4) == 0 || (moduleIndex % 4) == 3)
228 continue;
229 }
230
231 Int_t currentLayer = 0;
232 if (moduleIndex >= 80)
233 currentLayer = 1;
234
235 treeD->GetEvent(moduleIndex);
236
237 // get number of digits in this module
238 Int_t ndigitsInModule = digits->GetEntriesFast();
239
240 // get number of digits in each chip of the module
241 Int_t ndigitsInChip[5];
242 Bool_t hitByPrimary[5];
243 for( Int_t iChip=0; iChip<5; iChip++)
244 {
245 ndigitsInChip[iChip]=0;
246 hitByPrimary[iChip] = kFALSE;
247 }
248
249 // loop over digits in this module
250 for (Int_t iDig=0; iDig<ndigitsInModule; iDig++)
251 {
252 AliITSdigitSPD* dp = (AliITSdigitSPD*) digits->At(iDig);
253 Int_t column = dp->GetCoord1();
254 Int_t isChip = column / 32;
255
256 //printf("Digit %d has column %d which translates to chip %d\n", iDig, column, isChip);
257
258 fChipsFired->Fill(moduleIndex, isChip);
259
260 ndigitsInChip[isChip]++;
261
262 Bool_t debug = kFALSE;
263
264 // find out which particle caused this chip to fire
265 // if we find at least one primary we consider this chip being fired by a primary
266 for (Int_t track=0; track<10; ++track)
267 {
268 Int_t label = dp->GetTrack(track);
269
270 if (label < 0)
271 continue;
272
273 if (debug)
274 printf("track %d contains label %d\n", track, label);
275
276 TParticle* particle = stack->Particle(label);
277
278 if (!particle)
279 {
280 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (digit loop).", label));
281 continue;
282 }
283
284 if (debug)
285 {
286 particle->Print();
287 printf("statuscode = %d, p = %f, m = %d\n", particle->GetStatusCode(), particle->P(), particle->GetFirstMother());
288 }
289
290 // TODO delta electrons should be traced back to their mother. this is e.g. solved in AliITSClusterFinderV2::CheckLabels2
291 while (particle->P() < 0.02 && particle->GetStatusCode() == 0 && particle->GetFirstMother() >= 0)
292 {
293 label = particle->GetFirstMother();
294 particle = stack->Particle(label);
295
296 if (!particle)
297 break;
298
299 if (debug)
300 {
301 printf("-->\n");
302 printf("statuscode = %d, p = %f, m = %d\n", particle->GetStatusCode(), particle->P(), particle->GetFirstMother());
303 particle->Print();
304 }
305 }
306
307 if (!particle)
308 {
309 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (digit loop, finding delta electrons).", label));
310 continue;
311 }
312
313 if (label > nPrim)
314 continue;
315
316 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
317 continue;
318
319 if (debug)
320 printf("This was a primary (or delta-electron of a primary)!\n");
321
322 hitByPrimary[isChip] = kTRUE;
323 }
324 }
325
326 // get number of FOs in the module
327 for (Int_t ifChip=0; ifChip<5; ifChip++)
328 if( ndigitsInChip[ifChip] >= 1 )
329 {
330 totalNumberOfFO[currentLayer]++;
331 if (hitByPrimary[ifChip])
332 {
333 chipsHitByPrimaries[currentLayer]++;
334 }
335 //else
336 // chipsHitBySecondaries[currentLayer]++;
337 }
338 }
339
340 //printf("Fired chips: %d %d\n", totalNumberOfFOLayer1, totalNumberOfFOLayer2);
341
342 // now find clusters
343 Int_t clustersLayer[2];
344 clustersLayer[0] = 0;
345 clustersLayer[1] = 0;
346
347 loader->LoadRecPoints("READ");
348 TTree* treeR = loader->TreeR();
349 if (!treeR)
350 {
351 AliDebug(AliLog::kError, "Could not retrieve TreeR of ITS");
352 return kFALSE;
353 }
354
355 // TODO branches!
356 //treeR->SetBranchStatus("*", 0);
357
358 TClonesArray* itsClusters = 0;
359 treeR->SetBranchAddress("ITSRecPoints", &itsClusters);
360
361 Int_t nTreeEntries = treeR->GetEntries();
362 for (Int_t iEntry = 0; iEntry < nTreeEntries; ++iEntry)
363 {
364 if (!treeR->GetEvent(iEntry))
365 continue;
366
367 Int_t nClusters = itsClusters->GetEntriesFast();
368
369 while(nClusters--)
370 {
371 AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters);
372
373 if (cluster->GetLayer() == 0)
374 {
375 clustersLayer[0]++;
376 fClusterZL1->Fill(cluster->GetZ());
377 }
378 else if (cluster->GetLayer() == 1)
379 {
380 clustersLayer[1]++;
381 fClusterZL2->Fill(cluster->GetZ());
382 }
383 }
384 }
385
386 fPrimaryL1->Fill(multiplicity21, totalNumberOfFO[0], chipsHitByPrimaries[0], clustersLayer[0]);
387 fPrimaryL2->Fill(multiplicity16, totalNumberOfFO[1], chipsHitByPrimaries[1], clustersLayer[1]);
388
389 return kTRUE;
390}
391
392Bool_t AliHighMultiplicitySelector::Notify()
393{
a9aed06c 394 // get next ITS runloader
395
9608df41 396 AliRunLoader* runLoader = GetRunLoader();
397
398 if (runLoader)
399 {
400 AliITSLoader* loader = (AliITSLoader* )runLoader->GetLoader( "ITSLoader" );
401 if (loader)
402 {
403 loader->UnloadDigits();
404 loader->UnloadRecPoints();
405 }
406 }
407
408 return AliSelectorRL::Notify();
409}
410
411void AliHighMultiplicitySelector::SlaveTerminate()
412{
413 // The SlaveTerminate() function is called after all entries or objects
414 // have been processed. When running with PROOF SlaveTerminate() is called
415 // on each slave server.
416
417 AliSelectorRL::SlaveTerminate();
418
419 // Add the histograms to the output on each slave server
420 if (!fOutput)
421 {
422 AliDebug(AliLog::kError, "ERROR: Output list not initialized.");
423 return;
424 }
425
426 fOutput->Add(fChipsFired);
427 fOutput->Add(fPrimaryL1);
428 fOutput->Add(fPrimaryL2);
429 fOutput->Add(fClusterZL1);
430 fOutput->Add(fClusterZL2);
431}
432
433void AliHighMultiplicitySelector::Terminate()
434{
435 // The Terminate() function is the last function to be called during
436 // a query. It always runs on the client, it can be used to present
437 // the results graphically or save the results to file.
438
439 AliSelectorRL::Terminate();
440
441 fChipsFired = dynamic_cast<TH2F*> (fOutput->FindObject("fChipsFired"));
442 fPrimaryL1 = dynamic_cast<TNtuple*> (fOutput->FindObject("fPrimaryL1"));
443 fPrimaryL2 = dynamic_cast<TNtuple*> (fOutput->FindObject("fPrimaryL2"));
444 fClusterZL1 = dynamic_cast<TH1F*> (fOutput->FindObject("fClusterZL1"));
445 fClusterZL2 = dynamic_cast<TH1F*> (fOutput->FindObject("fClusterZL2"));
446
447 if (!fClusterZL1 || !fClusterZL2 || !fChipsFired || !fPrimaryL1 || !fPrimaryL2)
448 {
449 AliError("Histograms not available");
450 return;
451 }
452
453 WriteHistograms();
454}
455
456void AliHighMultiplicitySelector::WriteHistograms(const char* filename)
457{
a9aed06c 458 // write the histograms to a file
459
9608df41 460 TFile* file = TFile::Open(filename, "RECREATE");
461
462 fChipsFired->Write();
463 fPrimaryL1->Write();
464 fPrimaryL2->Write();
465 fClusterZL1->Write();
466 fClusterZL2->Write();
467
468 file->Close();
469}
470
471void AliHighMultiplicitySelector::ReadHistograms(const char* filename)
472{
a9aed06c 473 // read the data from a file and fill histograms
474
9608df41 475 TFile* file = TFile::Open(filename);
476
477 if (!file)
478 return;
479
480 fPrimaryL1 = dynamic_cast<TNtuple*> (file->Get("fPrimaryL1"));
481 fPrimaryL2 = dynamic_cast<TNtuple*> (file->Get("fPrimaryL2"));
482 fChipsFired = dynamic_cast<TH2F*> (file->Get("fChipsFired"));
483 fClusterZL1 = dynamic_cast<TH1F*> (file->Get("fClusterZL1"));
484 fClusterZL2 = dynamic_cast<TH1F*> (file->Get("fClusterZL2"));
485
486 #define MULT 1001, -0.5, 1000.5
487 #define BINNING_LAYER1 401, -0.5, 400.5
488 #define BINNING_LAYER2 801, -0.5, 800.5
489
490 fChipsLayer1 = new TH1F("fChipsLayer1", "Layer 1;Fired Chips;Count", BINNING_LAYER1);
491 fChipsLayer2 = new TH1F("fChipsLayer2", "Layer 2;Fired Chips;Count", BINNING_LAYER2);
492
493 fL1vsL2 = new TH2F("fL1vsL2", ";Fired Chips Layer 1;Fired Chips Layer 2", BINNING_LAYER1, BINNING_LAYER2);
494 fMvsL1 = new TH2F("fMvsL1", ";true multiplicity;fired chips layer1", MULT, BINNING_LAYER1);
495 fMvsL2 = new TH2F("fMvsL2", ";true multiplicity;fired chips layer2", MULT, BINNING_LAYER2);
496
497 fClvsL1 = new TH2F("fClvsL1", ";clusters layer1;fired chips layer1", MULT, BINNING_LAYER1);
498 fClvsL2 = new TH2F("fClvsL2", ";clusters layer2;fired chips layer2", MULT, BINNING_LAYER2);
499
500 for (Int_t i = 0; i < fPrimaryL1->GetEntries(); i++)
501 {
502 fPrimaryL1->GetEvent(i);
503 fPrimaryL2->GetEvent(i);
504
505 Int_t multiplicity21 = (Int_t) fPrimaryL1->GetArgs()[0];
506 Int_t multiplicity16 = (Int_t) fPrimaryL2->GetArgs()[0];
507
508 Int_t totalNumberOfFO[2];
509 totalNumberOfFO[0] = (Int_t) fPrimaryL1->GetArgs()[1];
510 totalNumberOfFO[1] = (Int_t) fPrimaryL2->GetArgs()[1];
511
512 Int_t chipsHitByPrimaries[2];
513 chipsHitByPrimaries[0] = (Int_t) fPrimaryL1->GetArgs()[2];
514 chipsHitByPrimaries[1] = (Int_t) fPrimaryL2->GetArgs()[2];
515
516 Int_t clustersLayer[2];
517 clustersLayer[0] = (Int_t) fPrimaryL1->GetArgs()[3];
518 clustersLayer[1] = (Int_t) fPrimaryL2->GetArgs()[3];
519
520 fChipsLayer1->Fill(totalNumberOfFO[0]);
521 fChipsLayer2->Fill(totalNumberOfFO[1]);
522
523 fL1vsL2->Fill(totalNumberOfFO[0], totalNumberOfFO[1]);
524
525 fMvsL1->Fill(multiplicity21, totalNumberOfFO[0]);
526 fMvsL2->Fill(multiplicity16, totalNumberOfFO[1]);
527
528 fClvsL1->Fill(clustersLayer[0], totalNumberOfFO[0]);
529 fClvsL2->Fill(clustersLayer[1], totalNumberOfFO[1]);
530 }
531}
532
f6093269 533TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut)
9608df41 534{
f6093269 535 //
536 // returns the trigger efficiency as function of multiplicity with a given cut
537 //
9608df41 538
539 //cut and multiply with x-section
f6093269 540 TH1* allEvents = multVsLayer->ProjectionX("fMvsL_x_total", 1, 1001);
9608df41 541 //allEvents->Sumw2();
542
543 //cut and multiply with x-section
f6093269 544 TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, 1001);
9608df41 545 //proj->Sumw2();
546
f6093269 547 //new TCanvas; allEvents->DrawCopy(); gPad->SetLogy();
548 //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
9608df41 549
550 // make probability distribution out of it
551 // TODO binomial errors do not work??? weird...
552 proj->Divide(proj, allEvents, 1, 1, "B");
553
f6093269 554 return proj;
555}
556
557TH1* AliHighMultiplicitySelector::GetXSectionCut(TH1* xSection, TH2* multVsLayer, Int_t cut)
558{
559 // returns the rel. cross section of the true spectrum that is measured when a cut at <cut> is performed
560
561 TH1* proj = GetTriggerEfficiency(multVsLayer, cut);
562
563 //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
9608df41 564
565 for (Int_t i=1; i<=proj->GetNbinsX(); i++)
566 {
567 if (i <= xSection->GetNbinsX())
568 {
569 Double_t value = proj->GetBinContent(i) * xSection->GetBinContent(i);
570 Double_t error = 0;
571
572 if (value != 0)
573 error = value * (proj->GetBinError(i) / proj->GetBinContent(i) + xSection->GetBinError(i) / xSection->GetBinContent(i));
574
575 proj->SetBinContent(i, value);
576 proj->SetBinError(i, error);
577 }
578 else
579 {
580 proj->SetBinContent(i, 0);
581 proj->SetBinError(i, 0);
582 }
583 }
584
f6093269 585 //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
9608df41 586
587 return proj;
588}
589
f6093269 590void AliHighMultiplicitySelector::MakeGraphs2(const char* title, TH1* xSection, TH2* fMvsL)
9608df41 591{
a9aed06c 592 // creates graphs
593
f6093269 594 TGraph* effGraph = new TGraph;
595 effGraph->SetTitle(Form("%s;Cut on fired chips;mult. where eff. >50%%", title));
596 TGraph* ratioGraph = new TGraph;
597 ratioGraph->SetTitle(Form("%s;Cut on fired chips;x-section_(>=eff. limit) / x-section_(total)", title));
598 TGraph* totalGraph = new TGraph;
599 totalGraph->SetTitle(Form("%s;Cut on fired chips;rel x-section_(>=eff. limit)", title));
600
601 for (Int_t cut = 0; cut <= 300; cut+=50)
602 {
603 TH1* proj = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("clone");
604
605 //proj->Rebin(3);
606 //proj->Scale(1.0 / 3);
607
608 new TCanvas; proj->DrawCopy();
9608df41 609
f6093269 610 Int_t limitBin = proj->GetNbinsX()+1;
611 while (limitBin > 1 && proj->GetBinContent(limitBin-1) > 0.5)
612 limitBin--;
613
614 Float_t limit = proj->GetXaxis()->GetBinCenter(limitBin);
615
616 effGraph->SetPoint(effGraph->GetN(), cut, limit);
617
618 proj = GetXSectionCut(xSection, fMvsL, cut);
619
620 Double_t ratio = 0;
621 Double_t total = 0;
622 if (proj->Integral(1, 1001) > 0)
623 {
624 ratio = proj->Integral(proj->FindBin(limit), 1001) / proj->Integral(1, 1001);
625 total = proj->Integral(proj->FindBin(limit), 1001);
626 }
627
628 ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio);
629 totalGraph->SetPoint(totalGraph->GetN(), cut, total);
630
631 Printf("Cut at %d --> trigger eff. is > 0.5 for mult. >= %.2f. That is the case for %f of the triggered, %e of all events", cut, limit, ratio, total);
632 }
633
634 TCanvas* canvas = new TCanvas(Form("%s_Efficiency", title), Form("%s_Efficiency", title), 1200, 800);
635 canvas->Divide(2, 2);
636
637 canvas->cd(1);
638 effGraph->Draw("A*");
639
640 for (Int_t i=8; i<=10; ++i)
641 {
642 TLine* line = new TLine(0, xSection->GetMean() * i, 300, xSection->GetMean() * i);
643 line->Draw();
644 }
645
646 canvas->cd(2); ratioGraph->Draw("A*");
647 canvas->cd(3); gPad->SetLogy(); totalGraph->Draw("A*");
648
649 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
650}
651
652void AliHighMultiplicitySelector::MakeGraphs(const char* title, TH1* xSection, TH2* fMvsL, Int_t limit)
653{
9608df41 654 // relative x-section, once we have a collision
a9aed06c 655
9608df41 656 xSection->Scale(1.0 / xSection->Integral());
657
658 TGraph* ratioGraph = new TGraph;
659 ratioGraph->SetTitle(Form("%s;Cut on fired chips;x-section_(>=%d) / x-section_(total)", title, limit));
660 TGraph* totalGraph = new TGraph;
661 totalGraph->SetTitle(Form("%s;Cut on fired chips;rel x-section_(>=%d)", title, limit));
662
663 Double_t max = 0;
664 Int_t bestCut = -1;
665 Double_t bestRatio = -1;
666 Double_t bestTotal = -1;
667 Int_t fullCut = -1;
668 Double_t fullRatio = -1;
669 Double_t fullTotal = -1;
670
671 fMvsL->Sumw2();
672
673 for (Int_t cut = 50; cut <= 300; cut+=2)
674 {
675 TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
676
677 Double_t ratio = 0;
678 Double_t total = 0;
679 if (proj->Integral(1, 1000) > 0)
680 {
681 ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000);
682 total = proj->Integral(limit, 1000);
683 }
684
685 max = TMath::Max(max, total);
686
687 //printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio);
688
689 if (total < max * 0.9 && bestCut == -1)
690 {
691 bestCut = cut;
692 bestRatio = ratio;
693 bestTotal = total;
694 }
695
696 if (ratio == 1 && fullCut == -1)
697 {
698 fullCut = cut;
699 fullRatio = ratio;
700 fullTotal = total;
701 }
702
703 ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio);
704 totalGraph->SetPoint(totalGraph->GetN(), cut, total);
705 }
706
707 if (bestCut != -1)
708 printf("Best cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", bestCut, limit, bestTotal, limit, bestRatio);
709 if (fullCut != -1)
710 printf("100%% cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", fullCut, limit, fullTotal, limit, fullRatio);
711
712 TCanvas* canvas = new TCanvas(Form("%s_RatioXSection_%d", title, limit), Form("%s_RatioXSection_%d", title, limit), 600, 400);
713 ratioGraph->Draw("A*");
714 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
715
716 canvas = new TCanvas(Form("%s_TotalXSection_%d", title, limit), Form("%s_TotalXSection_%d", title, limit), 600, 400);
717 totalGraph->Draw("A*");
718 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
719}
720
721void AliHighMultiplicitySelector::JPRPlots()
722{
a9aed06c 723 // more plots
724
9608df41 725 /*
726
727 gSystem->Load("libPWG0base");
728 .L AliHighMultiplicitySelector.cxx+
729 x = new AliHighMultiplicitySelector();
730 x->ReadHistograms("highmult_hijing100k.root");
731 x->JPRPlots();
732
733 */
734
735 // get x-sections
736 TFile* file = TFile::Open("crosssectionEx.root");
737 if (!file)
738 return;
739
740 TH1* xSections[2];
741 xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
742 xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
743
744 for (Int_t i=0; i<2; ++i)
745 {
746 if (!xSections[i])
747 continue;
748
749 TH1* xSection = xSections[i];
750 TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
751 //Int_t cut = (i == 0) ? 164 : 150; // 8 times the mean
752 //Int_t cut = (i == 0) ? 178 : 166; // 9 times the mean
753 Int_t cut = (i == 0) ? 190 : 182; // 10 times the mean
754
f6093269 755 // limit is N times the mean
756 Int_t limit = (Int_t) (xSection->GetMean() * 10);
757
758 // 10^28 lum --> 1.2 kHz
759 // 10^31 lum --> 1200 kHz
760 Float_t rate = 1200e3;
761
762 // time in s
763 Float_t lengthRun = 1e6;
764
9608df41 765 xSection->SetStats(kFALSE);
f6093269 766 xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
9608df41 767 xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5));
768 xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 400 : 350);
f6093269 769 //xSection->GetYaxis()->SetTitle("relative cross section");
9608df41 770 xSection->GetYaxis()->SetTitleOffset(1.2);
771
9608df41 772 // relative x-section, once we have a collision
773 xSection->Scale(1.0 / xSection->Integral());
774
775 TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
776
777 Double_t ratio = 0;
778 Double_t total = 0;
779 if (proj->Integral(1, 1000) > 0)
780 {
781 ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000);
782 total = proj->Integral(limit, 1000);
783 }
784
785 printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio);
786
787 TCanvas* canvas = new TCanvas(Form("HMPlots_%d", i), Form("HMPlots_%d", i), 800, 600);
788 canvas->SetLogy();
789 xSection->DrawCopy();
790 proj->SetLineColor(2);
791 proj->SetStats(kFALSE);
792 proj->DrawCopy("SAME");
793
794 TLegend* legend = new TLegend(0.15, 0.15, 0.45, 0.3);
795 legend->SetFillColor(0);
796 legend->AddEntry(xSection, "no trigger");
797 legend->AddEntry(proj, Form("FO trigger > %d chips", cut));
798 legend->Draw();
799
800 TLine* line = new TLine(limit, xSection->GetMinimum() * 0.5, limit, xSection->GetMaximum() * 2);
801 line->SetLineWidth(2);
802 line->Draw();
803
804 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
805
806 TCanvas* canvas2 = new TCanvas(Form("HMPlots_%d_Random", i), Form("HMPlots_%d_Random", i), 800, 600);
f6093269 807 //canvas2->SetTopMargin(0.05);
808 //canvas2->SetRightMargin(0.05);
9608df41 809 canvas2->SetLogy();
f6093269 810 xSection->DrawCopy("HIST");
9608df41 811
f6093269 812 TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.3);
9608df41 813 legend2->SetFillColor(0);
814 legend2->AddEntry(xSection, "no trigger");
815
816 TH1* proj2 = (TH1*) proj->Clone("random");
817 proj2->Reset();
f6093269 818 // MB lengthRun s 100 Hz
819 Int_t nTrigger = (Int_t) (100 * lengthRun * proj->Integral(1, 1000));
9608df41 820 proj2->FillRandom(proj, nTrigger);
821
822 TH1* proj3 = (TH1*) proj2->Clone("random_clone");
823 proj3->Divide(proj);
824 proj3->Fit("pol0", "0", "");
825 proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0));
826
f6093269 827 /*
9608df41 828 proj2->DrawCopy("SAME");
829 legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio)));
f6093269 830 */
9608df41 831
832 proj2 = (TH1*) proj->Clone("random2");
833 proj2->Reset();
f6093269 834 // 10^31 lum --> 1200 kHz; lengthRun s
835 nTrigger = (Int_t) (rate * proj->Integral(1, 1000) * lengthRun);
9608df41 836 proj2->FillRandom(proj, nTrigger);
837
838 proj3 = (TH1*) proj2->Clone("random_clone2");
839 proj3->Divide(proj);
840 proj3->Fit("pol0", "0", "");
841 proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0));
842
843 proj2->SetLineColor(4);
f6093269 844 proj2->SetMarkerStyle(7);
845 proj2->SetMarkerColor(4);
846 proj2->DrawCopy("SAME P");
847 //legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio)));
848 legend2->AddEntry(proj2, Form("FO trigger > %d chips", cut));
9608df41 849
850 legend2->Draw();
851 line->Draw();
852
853 canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
f6093269 854 canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
9608df41 855 }
856}
857
d1594f8a 858void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
f91bc12d 859{
860 //
861 // produces a spectrum created with N triggers
862 // number of triggers and thresholds for the moment fixed
863 //
864
865 /*
866
dd367a14 867 gSystem->Load("libANALYSIS");
f91bc12d 868 gSystem->Load("libPWG0base");
869 .L AliHighMultiplicitySelector.cxx+g
870 x = new AliHighMultiplicitySelector();
871 x->ReadHistograms("highmult_hijing100k.root");
872 x->Ntrigger();
873
d1594f8a 874 gSystem->Load("libPWG0base");
875 .L AliHighMultiplicitySelector.cxx+g
876 x = new AliHighMultiplicitySelector();
877 x->ReadHistograms("highmult_hijing100k.root");
878 x->Ntrigger(kFALSE);
879
f91bc12d 880 */
881
882 // get x-sections
883 TFile* file = TFile::Open("crosssectionEx.root");
884 if (!file)
885 return;
886
887 TH1* xSections[2];
888 xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
889 xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
890
51f6de65 891 // 10^28 lum --> 1.4 kHz
892 // 10^31 lum --> 1400 kHz
893 //Float_t rate = 1400e3;
894 Float_t rate = 1.4e3;
f91bc12d 895
896 // time in s
51f6de65 897 Float_t lengthRun = 1e5;
f91bc12d 898
899 Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
900 Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
901
902 // put to 2 for second layer
903 for (Int_t i=0; i<1; ++i)
904 {
905 if (!xSections[i])
906 continue;
907
908 TH1* xSection = xSections[i];
909 TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
910
a9aed06c 911 //Int_t nCuts = 6;
912 //Int_t cuts[] = { 0, 164, 178, 190, 204, 216 };
f91bc12d 913 //Int_t nCuts = 4;
914 //Int_t cuts[] = { 0, 164, 190, 216 };
d1594f8a 915
916 //Int_t nCuts = 4;
917 //Int_t cuts[] = { 0, 114, 145, 165 };
918 //Float_t ratePerTrigger[] = { 60, 13.3, 13.3, 13.3 };
919
920 Int_t nCuts = 3;
51f6de65 921 Int_t cuts[] = { 0, 114, 148 };
922
923 //Int_t nCuts = 3;
924 //Int_t cuts[] = { 0, 126, 162 };
925 //Float_t ratePerTrigger[] = { 60, 20.0, 20.0 };
a9aed06c 926
f91bc12d 927 // desired trigger rate in Hz
51f6de65 928 Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 };
f91bc12d 929
930 xSection->SetStats(kFALSE);
931 xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
932 xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5));
933 xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 450 : 350);
934 //xSection->GetYaxis()->SetTitle("relative cross section");
935 xSection->GetYaxis()->SetTitleOffset(1.2);
936
937 // relative x-section, once we have a collision
938 xSection->Scale(1.0 / xSection->Integral());
939
940 TCanvas* canvas2 = new TCanvas(Form("HMPlots2_%d_Random", i), Form("HMPlots2_%d_Random", i), 800, 600);
941 canvas2->SetTopMargin(0.05);
942 canvas2->SetRightMargin(0.05);
943 canvas2->SetLogy();
d1594f8a 944
945 if (relative)
946 xSection->DrawCopy("HIST");
f91bc12d 947
948 TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.4);
949 legend2->SetFillColor(0);
950 legend2->AddEntry(xSection, "cross-section");
951
952 for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut)
953 {
954 Int_t cut = cuts[currentCut];
955
d1594f8a 956 TH1* triggerEff = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
957
f91bc12d 958 TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
959
d1594f8a 960 Float_t triggerLimit = 0;
961 for (Int_t bin = 1; bin <= triggerEff->GetNbinsX(); bin++)
962 if (triggerEff->GetBinContent(bin) < 0.5)
963 triggerLimit = triggerEff->GetXaxis()->GetBinCenter(bin);
964
965 Printf("Efficiency limit (50%%) is at multiplicity %f", triggerLimit);
966
f91bc12d 967 Double_t total = 0;
968 if (proj->Integral(1, 1000) > 0)
969 total = proj->Integral(1, 1000);
970
971 printf("Cut at %d: rel. x-section = %e\n", cut, total);
972
973 TH1* proj2 = (TH1*) proj->Clone("random2");
974 proj2->Reset();
975
976 // calculate downscale factor
977 Float_t normalRate = rate * proj->Integral(1, 1000);
978 Float_t downScale = normalRate / ratePerTrigger[currentCut];
979 if (downScale < 1)
980 downScale = 1;
981 Long64_t nTrigger = (Long64_t) (normalRate / downScale * lengthRun);
dd367a14 982 nTrigger = TMath::Nint(((Float_t) nTrigger) / 1000) * 1000;
f91bc12d 983
984 Printf("Normal rate is %f, downscale: %f, Simulating %lld triggers", normalRate, downScale, nTrigger);
51f6de65 985 if (nTrigger == 0)
986 continue;
987
f91bc12d 988 proj2->FillRandom(proj, nTrigger);
989
990 Int_t removed = 0;
991 for (Int_t bin=1; bin<proj2->GetNbinsX(); ++bin)
992 if (proj2->GetBinContent(bin) < 5)
993 {
994 removed += (Int_t) proj2->GetBinContent(bin);
995 proj2->SetBinContent(bin, 0);
996 }
997
998 Printf("Removed %d events", removed);
999
d1594f8a 1000 if (relative)
1001 proj2->Scale(1.0 / nTrigger * proj->Integral(1, 1000));
f91bc12d 1002
1003 proj2->SetLineColor(colors[currentCut]);
1004 proj2->SetMarkerStyle(markers[currentCut]);
1005 proj2->SetMarkerColor(colors[currentCut]);
d1594f8a 1006
1007 if (relative || currentCut > 0) {
1008 proj2->DrawCopy("SAME P");
1009 } else
1010 proj2->DrawCopy(" P");
1011
745d6088 1012 TString eventStr;
1013 if (nTrigger > 1e6)
1014 {
1015 eventStr.Form("%lld M", nTrigger / 1000 / 1000);
1016 }
1017 else if (nTrigger > 1e3)
1018 {
1019 eventStr.Form("%lld K", nTrigger / 1000);
1020 }
1021 else
1022 eventStr.Form("%lld", nTrigger);
1023
1024 TString triggerStr;
1025 if (cut == 0)
1026 {
1027 triggerStr = "minimum bias";
1028 }
1029 else
1030 triggerStr.Form("FO > %d chips", cut);
d1594f8a 1031
745d6088 1032 legend2->AddEntry(proj2, Form("%s evts, %s", eventStr.Data(), triggerStr.Data()));
1033
1034 if (triggerLimit > 1)
1035 {
1036 TLine* line = new TLine(triggerLimit, proj2->GetMinimum(), triggerLimit, proj2->GetMaximum());
1037 line->SetLineColor(colors[currentCut]);
1038 line->Draw();
1039 }
f91bc12d 1040 }
1041
1042 legend2->Draw();
1043
1044 canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
1045 canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
1046 }
1047}
1048
d1594f8a 1049void AliHighMultiplicitySelector::Contamination()
1050{
1051 //
d1594f8a 1052 //
1053
1054 /*
1055
1056 gSystem->Load("libANALYSIS");
1057 gSystem->Load("libPWG0base");
1058 .L AliHighMultiplicitySelector.cxx+g
1059 x = new AliHighMultiplicitySelector();
1060 x->ReadHistograms("highmult_hijing100k.root");
1061 x->Contamination();
1062
1063 */
1064
1065 // get x-sections
1066 TFile* file = TFile::Open("crosssectionEx.root");
1067 if (!file)
1068 return;
1069
1070 TH1* xSections[2];
1071 xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1072 xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1073
1074 // rate = L * sigma
1075 // sigma ~ 80 mb (Pythia 14 TeV)
1076 // 10^28 lum --> 8e2 Hz
1077 // 10^31 lum --> 8e5 Hz
1078 Double_t rates[] = { 8e2, 8e3, 8e4, 8e5 };
1079
1080 Int_t nCuts = 4;
1081 Int_t cuts[] = { 104, 134, 154, 170 };
1082
1083 // put to 2 for second layer
1084 for (Int_t i=0; i<1; ++i)
1085 {
1086 if (!xSections[i])
1087 continue;
1088
1089 // relative x-section, once we have a collision
1090 xSections[i]->Scale(1.0 / xSections[i]->Integral());
1091
1092 Int_t max = xSections[i]->GetNbinsX();
1093 max = 500;
1094
1095 Float_t* xSection = new Float_t[max];
1096 for (Int_t mult = 0; mult < max; mult++)
1097 xSection[mult] = xSections[i]->GetBinContent(mult+1);
1098
1099 TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
1100
1101 TGraph* graph = new TGraph;
1102
1103 for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut)
1104 {
1105 Int_t cut = cuts[currentCut];
1106 Double_t rate = rates[currentCut];
1107 //Double_t rate = rates[3];
1108
1109 // coll. in 100 ns window
69b09e3b 1110 Double_t windowSize = 100e-9;
1111 //Double_t windowSize = 25e-9;
d1594f8a 1112 Double_t collPerWindow = windowSize * rate;
1113 Printf("coll/window = %f", collPerWindow);
1114 Double_t windowsPerSecond = 1.0 / windowSize;
1115
1116 TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
1117 Float_t* triggerEff = new Float_t[max];
1118 for (Int_t mult = 0; mult < max; mult++)
1119 triggerEff[mult] = triggerEffHist->GetBinContent(mult+1);
1120
1121 Double_t triggerRate = 0;
1122 for (Int_t mult = 0; mult < max; mult++)
1123 triggerRate += xSection[mult] * triggerEff[mult];
1124
1125 triggerRate *= TMath::Poisson(1, collPerWindow) * windowsPerSecond;
1126
1127 Printf("Rate for 1 collision is %f Hz", triggerRate);
1128
1129 Double_t triggerRate2 = 0;
1130 for (Int_t mult = 0; mult < max; mult++)
1131 for (Int_t mult2 = mult; mult2 < max; mult2++)
1132 if (mult+mult2 < max)
1133 triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2];
1134
1135 triggerRate2 *= TMath::Poisson(2, collPerWindow) * windowsPerSecond;
d1594f8a 1136
1137 Printf("Rate for 2 collisions is %f Hz --> %.1f%%", triggerRate2, triggerRate2 / triggerRate * 100);
1138
1139 Double_t triggerRate3 = 0;
144ff489 1140
d1594f8a 1141 for (Int_t mult = 0; mult < max; mult++)
1142 for (Int_t mult2 = mult; mult2 < max-mult; mult2++)
1143 for (Int_t mult3 = 0; mult3 < max-mult-mult2; mult3++)
1144 //if (mult+mult2+mult3 < max)
1145 triggerRate3 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * xSection[mult3] * triggerEff[mult+mult2+mult3];
1146
1147 triggerRate3 *= TMath::Poisson(3, collPerWindow) * windowsPerSecond;
1148 //triggerRate3 *= collPerWindow * collPerWindow * rate;
1149
1150 Printf("Rate for 3 collisions is %f Hz --> %.1f%%", triggerRate3, triggerRate3 / triggerRate * 100);
1151
1152 Float_t totalContamination = (triggerRate2 + triggerRate3) / (triggerRate + triggerRate2 + triggerRate3);
1153
1154 Printf("Total contamination is %.1f%%", totalContamination * 100);
1155
1156 graph->SetPoint(graph->GetN(), cut, totalContamination);
1157
1158 continue;
1159
1160 Double_t triggerRate4 = 0;
1161 for (Int_t mult = 0; mult < max; mult++)
1162 for (Int_t mult2 = mult; mult2 < max-mult; mult2++)
1163 for (Int_t mult3 = 0; mult3 < max-mult-mult2; mult3++)
1164 for (Int_t mult4 = 0; mult4 < max-mult-mult2-mult3; mult4++)
1165 //if (mult+mult2+mult3+mult4 < max)
1166 triggerRate4 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * xSection[mult3] * xSection[mult4] * triggerEff[mult+mult2+mult3+mult4];
1167
1168 //triggerRate4 *= collPerWindow * collPerWindow * collPerWindow * rate;
1169 triggerRate4 *= TMath::Poisson(4, collPerWindow) * windowsPerSecond;
1170
1171 Printf("Rate for 4 collisions is %f Hz --> %.1f%%", triggerRate4, triggerRate4 / triggerRate * 100);
1172
1173 // general code for n collisions follows, however much slower...
1174 /*
1175 const Int_t maxdepth = 4;
1176 for (Int_t depth = 1; depth <= maxdepth; depth++) {
1177 Double_t triggerRate = 0;
1178
1179 Int_t m[maxdepth];
1180 for (Int_t d=0; d<maxdepth; d++)
1181 m[d] = 0;
1182
1183 while (m[0] < max) {
1184 Double_t value = 1;
1185 Int_t sum = 0;
1186 for (Int_t d=0; d<depth; d++) {
1187 value *= xSection[m[d]];
1188 sum += m[d];
1189 }
1190
1191 if (sum < max) {
1192 value *= triggerEff[sum];
1193 triggerRate += value;
1194 }
1195
1196 Int_t increase = depth-1;
1197 ++m[increase];
1198 while (m[increase] == max && increase > 0) {
1199 m[increase] = 0;
1200 --increase;
1201 ++m[increase];
1202 }
1203 }
1204
1205 triggerRate *= rate * TMath::Power(collPerWindow, depth - 1);
1206
1207 Printf("Rate for %d collisions is %f Hz", depth, triggerRate);
1208 }*/
1209 }
1210
1211 new TCanvas; graph->Draw("AP*");
1212 }
1213}
1214
144ff489 1215void AliHighMultiplicitySelector::Contamination2()
1216{
1217 //
1218 // produces a spectrum created with N triggers
1219 // number of triggers and thresholds for the moment fixed
1220 //
1221
1222 /*
1223
1224 gSystem->Load("libANALYSIS");
1225 gSystem->Load("libPWG0base");
1226 .L AliHighMultiplicitySelector.cxx+g
1227 x = new AliHighMultiplicitySelector();
1228 x->ReadHistograms("highmult_hijing100k.root");
1229 x->Contamination2();
1230
1231 */
1232
1233 // get x-sections
1234 TFile* file = TFile::Open("crosssectionEx.root");
1235 if (!file)
1236 return;
1237
1238 TH1* xSections[2];
1239 xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1240 xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1241
1242 Int_t nCuts = 4;
1243 Int_t cuts[] = { 104, 134, 154, 170 };
1244
1245 new TCanvas;
1246
1247 Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
1248 Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
1249
1250 // put to 2 for second layer
1251 for (Int_t i=0; i<1; ++i)
1252 {
1253 if (!xSections[i])
1254 continue;
1255
1256 // relative x-section, once we have a collision
1257 xSections[i]->Scale(1.0 / xSections[i]->Integral());
1258
1259 Int_t max = xSections[i]->GetNbinsX();
1260 max = 500;
1261
1262 Float_t* xSection = new Float_t[max];
1263 for (Int_t mult = 0; mult < max; mult++)
1264 xSection[mult] = xSections[i]->GetBinContent(mult+1);
1265
1266 TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
1267
1268 for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut)
1269 {
1270 TGraph* graph = new TGraph;
1271 graph->SetMarkerColor(colors[currentCut]);
1272 graph->SetMarkerStyle(markers[currentCut]);
1273
1274 Int_t cut = cuts[currentCut];
1275
1276 TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
1277 Float_t* triggerEff = new Float_t[max];
1278 for (Int_t mult = 0; mult < max; mult++)
1279 triggerEff[mult] = triggerEffHist->GetBinContent(mult+1);
1280
1281 Double_t triggerRate = 0;
1282 for (Int_t mult = 0; mult < max; mult++)
1283 triggerRate += xSection[mult] * triggerEff[mult];
1284
1285 Printf("Raw value for 1 collision is %e", triggerRate);
1286
1287 Double_t triggerRate2 = 0;
1288 for (Int_t mult = 0; mult < max; mult++)
1289 for (Int_t mult2 = mult; mult2 < max; mult2++)
1290 if (mult+mult2 < max)
1291 triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2];
1292
1293 Printf("Raw value for 2 collisions is %e", triggerRate2);
1294
69b09e3b 1295 for (Double_t doubleRate = 0; doubleRate <= 0.3; doubleRate += 0.005)
144ff489 1296 {
1297 Float_t totalContamination = (triggerRate2 * doubleRate) / (triggerRate + triggerRate2 * doubleRate);
1298
1299 //Printf("Total contamination is %.1f%%", totalContamination * 100);
1300
1301 graph->SetPoint(graph->GetN(), doubleRate, totalContamination);
1302 }
1303
1304 graph->Draw((currentCut == 0) ? "A*" : "* SAME");
1305 graph->GetXaxis()->SetRangeUser(0, 1);
1306 }
1307 }
1308}
1309
69b09e3b 1310void AliHighMultiplicitySelector::Contamination3()
1311{
1312 //
1313 //
1314
1315 /*
1316
1317 gSystem->Load("libANALYSIS");
1318 gSystem->Load("libPWG0base");
1319 .L AliHighMultiplicitySelector.cxx+g
1320 x = new AliHighMultiplicitySelector();
1321 x->ReadHistograms("highmult_hijing100k.root");
1322 x->Contamination3();
1323
1324 */
1325
1326 // get x-sections
1327 TFile* file = TFile::Open("crosssectionEx.root");
1328 if (!file)
1329 return;
1330
1331 TH1* xSections[2];
1332 xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1333 xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1334
1335 // prob for a collision in a bunch crossing
1336 Int_t nRates = 3;
1337 Float_t rates[] = {0.07, 0.1, 0.2};
1338
1339 new TCanvas;
1340
1341 Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
1342 Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
1343
1344 // put to 2 for second layer
1345 for (Int_t i=0; i<1; ++i)
1346 {
1347 if (!xSections[i])
1348 continue;
1349
1350 // relative x-section, once we have a collision
1351 xSections[i]->Scale(1.0 / xSections[i]->Integral());
1352
1353 Int_t max = xSections[i]->GetNbinsX();
1354 max = 500;
1355
1356 Float_t* xSection = new Float_t[max];
1357 for (Int_t mult = 0; mult < max; mult++)
1358 xSection[mult] = xSections[i]->GetBinContent(mult+1);
1359
1360 TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
1361
1362 for (Int_t currentRate = 0; currentRate<nRates; ++currentRate)
1363 {
1364 TGraph* graph = new TGraph;
1365 graph->SetMarkerColor(colors[currentRate]);
1366 graph->SetMarkerStyle(markers[currentRate]);
1367
1368 Float_t rate = rates[currentRate];
1369
1370 for (Int_t cut = 100; cut <= 201; cut += 10)
1371 {
1372 Printf("cut at %d", cut);
1373
1374 TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
1375 Float_t* triggerEff = new Float_t[max];
1376 for (Int_t mult = 0; mult < max; mult++)
1377 triggerEff[mult] = triggerEffHist->GetBinContent(mult+1);
1378
1379 Double_t triggerRate = 0;
1380 for (Int_t mult = 0; mult < max; mult++)
1381 triggerRate += xSection[mult] * triggerEff[mult];
1382
1383 Printf("Raw value for 1 collision is %e", triggerRate);
1384
1385 Double_t triggerRate2 = 0;
1386 for (Int_t mult = 0; mult < max; mult++)
1387 for (Int_t mult2 = mult; mult2 < max; mult2++)
1388 if (mult+mult2 < max)
1389 triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2];
1390
1391 Printf("Raw value for 2 collisions is %e", triggerRate2);
1392
1393 Double_t triggerRate3 = 0;
1394 for (Int_t mult = 0; mult < max; mult++)
1395 for (Int_t mult2 = 0; mult2 < max; mult2++)
1396 for (Int_t mult3 = 0; mult3 < max; mult3++)
1397 if (mult+mult2+mult3 < max)
1398 triggerRate3 += xSection[mult] * xSection[mult2] * xSection[mult3] * triggerEff[mult+mult2+mult3];
1399
1400 Printf("Raw value for 3 collisions is %e", triggerRate3);
1401
1402 Double_t singleRate = TMath::Poisson(1, rate);
1403 Double_t doubleRate = TMath::Poisson(2, rate);
1404 Double_t tripleRate = TMath::Poisson(3, rate);
1405
1406 Printf("single = %f, double = %f, triple = %f", singleRate, doubleRate, tripleRate);
1407
1408 Float_t totalContamination = (triggerRate2 * doubleRate + triggerRate3 * tripleRate) / (triggerRate * singleRate + triggerRate2 * doubleRate + triggerRate3 * tripleRate);
1409
1410 //Printf("Total contamination is %.1f%%", totalContamination * 100);
1411
1412 graph->SetPoint(graph->GetN(), cut, totalContamination);
1413 }
1414
1415 graph->Draw((currentRate == 0) ? "A*" : "* SAME");
1416 graph->GetXaxis()->SetTitle("layer 1 threshold");
1417 graph->GetYaxis()->SetTitle("contamination");
1418 graph->GetYaxis()->SetRangeUser(0, 1);
1419 }
1420 }
1421}
1422
9608df41 1423void AliHighMultiplicitySelector::DrawHistograms()
1424{
a9aed06c 1425 // draws the histograms
1426
9608df41 1427 /*
1428
1429 gSystem->Load("libPWG0base");
1430 .L AliHighMultiplicitySelector.cxx+
1431 x = new AliHighMultiplicitySelector();
1432 x->ReadHistograms("highmult_pythia.root");
1433 x->DrawHistograms();
1434
1435
1436 gSystem->Load("libPWG0base");
1437 .L AliHighMultiplicitySelector.cxx+
1438 x = new AliHighMultiplicitySelector();
1439 x->ReadHistograms("highmult_hijing.root");
1440 x->DrawHistograms();
1441
1442 gSystem->Load("libPWG0base");
1443 .L AliHighMultiplicitySelector.cxx+
1444 x = new AliHighMultiplicitySelector();
1445 x->ReadHistograms("highmult_central.root");
1446 x->DrawHistograms();
1447
745d6088 1448 gSystem->Load("libANALYSIS");
9608df41 1449 gSystem->Load("libPWG0base");
1450 .L AliHighMultiplicitySelector.cxx+
1451 x = new AliHighMultiplicitySelector();
1452 x->ReadHistograms("highmult_hijing100k.root");
1453 x->DrawHistograms();
1454
1455 */
1456
1457 /*TCanvas* canvas = new TCanvas("chips", "chips", 600, 400);
1458
1459 fChipsLayer2->SetLineColor(2);
1460 fChipsLayer2->SetStats(kFALSE);
1461 fChipsLayer1->SetStats(kFALSE);
1462 fChipsLayer2->SetTitle("");
1463 fChipsLayer2->DrawCopy();
1464 fChipsLayer1->DrawCopy("SAME");
1465 canvas->SaveAs("chips.gif");
1466
1467 canvas = new TCanvas("L1vsL2", "L1vsL2", 600, 400);
1468 fL1vsL2->SetStats(kFALSE);
1469 fL1vsL2->DrawCopy("COLZ");
1470 gPad->SetLogz();
1471 canvas->SaveAs("L1vsL2.gif");*/
1472
f6093269 1473 TCanvas *canvas = new TCanvas("L1", "L1", 800, 600);
1474 canvas->SetTopMargin(0.05);
1475 canvas->SetRightMargin(0.12);
9608df41 1476 fMvsL1->SetStats(kFALSE);
1477 fMvsL1->DrawCopy("COLZ");
1478 gPad->SetLogz();
1479
1480 canvas->SaveAs("L1NoCurve.gif");
f6093269 1481 canvas->SaveAs("L1NoCurve.eps");
9608df41 1482
745d6088 1483 TLine* line = new TLine(fMvsL1->GetXaxis()->GetXmin(), 150, fMvsL1->GetXaxis()->GetXmax(), 150);
1484 line->SetLineWidth(2);
1485 line->SetLineColor(kRed);
1486 line->Draw();
1487
1488 canvas->SaveAs("L1NoCurveCut.gif");
1489 canvas->SaveAs("L1NoCurveCut.eps");
1490
1491 return;
1492
9608df41 1493 // draw corresponding theoretical curve
1494 TF1* func = new TF1("func", "[0]*(1-(1-1/[0])**x)", 1, 1000);
1495 func->SetParameter(0, 400-5*2);
1496 func->DrawCopy("SAME");
1497
1498 canvas->SaveAs("L1.gif");
1499
1500 canvas = new TCanvas("L2", "L2", 600, 400);
1501 //fMvsL2->GetYaxis()->SetRangeUser(0, 150);
1502 fMvsL2->SetStats(kFALSE);
1503 fMvsL2->DrawCopy("COLZ");
1504 gPad->SetLogz();
1505 func->SetParameter(0, 800-5*4);
1506 func->DrawCopy("SAME");
1507 canvas->SaveAs("L2.gif");
1508
f6093269 1509 // get x-sections
1510 TFile* file = TFile::Open("crosssectionEx.root");
1511 if (file)
1512 {
1513 TH1* xSection2 = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1514 TH1* xSection15 = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1515
1516 MakeGraphs2("Layer1", xSection2, fMvsL1);
1517 return;
1518
1519 // 5 times the mean
1520 //MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 5)); //75 * 2 * 2);
1521 //MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 5)); //(Int_t) (75 * 1.5 * 2));
1522
1523 MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 8));
1524 MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 8));
1525 MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 9));
1526 MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 9));
1527 MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 10));
1528 MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 10));
1529
1530 file->Close();
1531 }
1532
9608df41 1533 // make spread hists
1534 TGraph* spread = new TGraph;
1535 spread->SetTitle("Spread L1;true multiplicity;RMS");
1536
1537 for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
1538 {
1539 TH1* proj = fMvsL1->ProjectionY("proj", i, i);
1540 spread->SetPoint(spread->GetN(), i, proj->GetRMS());
1541 }
1542
1543 canvas = new TCanvas("SpreadL1", "SpreadL1", 600, 400);
1544 spread->Draw("A*");
1545 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1546
1547 TF1* log = new TF1("log", "[0]*log([1]*x)", 1, 150);
1548 log->SetParLimits(0, 0, 10);
1549 log->SetParLimits(1, 1e-5, 10);
1550
1551 spread->Fit(log, "", "", 1, 150);
1552 log->DrawCopy("SAME");
1553
1554 TGraph* spread2 = new TGraph;
1555 spread2->SetTitle("Spread L2;true multiplicity;RMS");
1556
1557 for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
1558 {
1559 TH1* proj = fMvsL2->ProjectionY("proj", i, i);
1560 spread2->SetPoint(spread2->GetN(), i, proj->GetRMS());
1561 }
1562
1563 canvas = new TCanvas("SpreadL2", "SpreadL2", 600, 400);
1564 spread2->Draw("A*");
1565 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1566
1567 spread2->Fit(log, "", "", 1, 150);
1568 log->DrawCopy("SAME");
1569
9608df41 1570 canvas = new TCanvas("Clusters_L1", "Clusters_L1", 600, 400);
1571 fClvsL1->SetStats(kFALSE);
1572 fClvsL1->DrawCopy("COLZ");
1573 gPad->SetLogz();
1574
1575 func->SetParameter(0, 400-5*2);
1576 func->DrawCopy("SAME");
1577
1578 canvas->SaveAs("Clusters_L1.gif");
1579
1580 canvas = new TCanvas("Clusters_L2", "Clusters_L2", 600, 400);
1581 fClvsL2->SetStats(kFALSE);
1582 fClvsL2->DrawCopy("COLZ");
1583 gPad->SetLogz();
1584 func->SetParameter(0, 800-5*4);
1585 func->DrawCopy("SAME");
1586 canvas->SaveAs("Clusters_L2.gif");
1587
1588 canvas = new TCanvas("ChipsFired", "ChipsFired", 600, 400);
1589 //fChipsFired->GetYaxis()->SetRangeUser(0, 150);
1590 fChipsFired->SetStats(kFALSE);
1591 fChipsFired->DrawCopy("COLZ");
1592 canvas->SaveAs("ChipsFired.gif");
1593
1594 /*TH1F* tresholdHistL1 = new TH1F("tresholdHistL1", ";chip treshold;<n>", BINNING_LAYER1);
1595 TH1F* tresholdHistL2 = new TH1F("tresholdHistL2", ";chip treshold;<n>", BINNING_LAYER2);
1596
1597 for (Int_t treshold = 0; treshold < 800; treshold++)
1598 {
1599 if (fPrimaryL1->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
1600 {
1601 TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL1->GetHistogram());
1602 if (mult)
1603 tresholdHistL1->Fill(treshold, mult->GetMean());
1604 }
1605 if (fPrimaryL2->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
1606 {
1607 TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL2->GetHistogram());
1608 if (mult)
1609 tresholdHistL2->Fill(treshold, mult->GetMean());
1610 }
1611 }
1612
1613 canvas = new TCanvas("TresholdL1", "TresholdL1", 600, 400);
1614 tresholdHistL1->Draw();
1615 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1616
1617 canvas = new TCanvas("TresholdL2", "TresholdL2", 600, 400);
1618 tresholdHistL2->Draw();
1619 canvas->SaveAs(Form("%s.gif", canvas->GetName()));*/
1620
1621 fPrimaryL1->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff1", "", "prof goff");
1622 fPrimaryL2->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff2", "", "prof goff");
1623
1624 canvas = new TCanvas("Efficiency", "Efficiency", 600, 400);
1625 fPrimaryL1->GetHistogram()->SetStats(kFALSE);
1626 fPrimaryL1->GetHistogram()->Draw();
1627 fPrimaryL2->GetHistogram()->SetLineColor(2);
1628 fPrimaryL2->GetHistogram()->SetStats(kFALSE);
1629 fPrimaryL2->GetHistogram()->Draw("SAME");
1630 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1631
1632 canvas = new TCanvas("ClustersZL1", "ClustersZL1", 600, 400);
1633 fClusterZL1->Rebin(2);
1634 fClusterZL1->Draw();
1635 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1636
1637 canvas = new TCanvas("ClustersZL2", "ClustersZL2", 600, 400);
1638 fClusterZL2->Draw();
1639 fClusterZL2->Rebin(2);
1640 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1641}
a9aed06c 1642
1643TGraph* AliHighMultiplicitySelector::IntFractRate()
1644{
1645 // A plot which shows the fractional rate above any threshold
1646 // as function of threshold (i.e. the integral of dSigma/dN as function of
1647 // N, normalised to 1 for N=0)
1648
1649 /*
d1594f8a 1650 gSystem->Load("libANALYSIS");
a9aed06c 1651 gSystem->Load("libPWG0base");
1652 .L AliHighMultiplicitySelector.cxx+
1653 x = new AliHighMultiplicitySelector();
1654 x->ReadHistograms("highmult_hijing100k.root");
1655 x->IntFractRate();
1656 */
1657
1658 // get x-sections
1659 TFile* file = TFile::Open("crosssectionEx.root");
1660 if (!file)
1661 return 0;
1662
1663 TH1* xSection;
1664 xSection = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1665
1666 TGraph* result = new TGraph;
1667
1668 for (Int_t threshold = 0; threshold < 300; threshold += 2)
1669 {
1670 TH1* proj = GetXSectionCut(xSection, fMvsL1, threshold);
1671
1672 //new TCanvas; proj->DrawCopy();
1673
1674 Double_t integral = proj->Integral();
1675
d1594f8a 1676 Printf("Cut at %d, integral is %e", threshold, integral);
a9aed06c 1677
1678 result->SetPoint(result->GetN(), threshold, integral);
1679 }
1680
1681 TCanvas* canvas = new TCanvas("IntFractRate", "IntFractRate", 600, 400);
1682 gPad->SetLogy();
1683
1684 result->Draw("A*");
1685 result->GetXaxis()->SetTitle("threshold");
1686 result->GetYaxis()->SetTitle("integrated fractional rate above threshold");
1687
1688 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1689
1690 return result;
1691}
d1594f8a 1692
1693void AliHighMultiplicitySelector::MBComparison()
1694{
1695 //
1696 // finds the threshold from which onwards the number of found events above N times the mean
1697 // is higher using a high mult. trigger than just triggering with MB
1698 //
1699
1700 /*
1701
1702 gSystem->Load("libANALYSIS");
1703 gSystem->Load("libPWG0base");
1704 .L AliHighMultiplicitySelector.cxx+g
1705 x = new AliHighMultiplicitySelector();
1706 x->ReadHistograms("highmult_hijing100k.root");
1707 x->MBComparison();
1708
1709 */
1710
1711 // get x-sections
1712 TFile* file = TFile::Open("crosssectionEx.root");
1713 if (!file)
1714 return;
1715
1716 TH1* xSections[2];
1717 xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1718 xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1719
1720 // rate = L * sigma
1721 // sigma ~ 80 mb (Pythia 14 TeV)
1722 // 10^28 lum --> 8e2 Hz
1723 // 10^31 lum --> 8e5 Hz
1724 Int_t nRates = 4;
1725 Double_t rates[] = { 8e2, 8e3, 8e4, 8e5 };
1726
1727 // threshold in number of fired chips for corresponding rate
1728 //Int_t cuts[] = { 104, 134, 154, 170 }; // values for 20 Hz
1729 Int_t cuts[] = { 82, 124, 147, 164 }; // values for 50 Hz
1730
1731 // bandwidth, fractions (for MB, high mult.)
1732 Float_t bandwidth = 1e3;
1733 Float_t fractionMB = 0.5;
1734 Float_t fractionHM = 0.05;
1735
1736 // different limits to define "interesting events"
1737 Int_t nLimits = 9;
1738 Int_t limits[] = { 0, 1, 2, 4, 6, 7, 8, 9, 10 };
1739
1740 // put to 2 for second layer
1741 for (Int_t i=0; i<1; ++i)
1742 {
1743 if (!xSections[i])
1744 continue;
1745
1746 TH1* xSection = xSections[i];
1747 TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
1748
1749 // relative x-section, once we have a collision
1750 xSection->Scale(1.0 / xSection->Integral());
1751
1752 xSection->SetStats(kFALSE);
1753 xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
1754 xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5));
1755 xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 450 : 350);
1756 xSection->GetYaxis()->SetTitleOffset(1.2);
1757
1758 TCanvas* canvas = new TCanvas("MBComparison", "MBComparison", 1000, 800);
1759 canvas->Divide(3, 3);
1760
1761 for (Int_t currentLimit = 0; currentLimit<nLimits; currentLimit++)
1762 {
1763 // limit is N times the mean
1764 Int_t limit = (Int_t) (xSection->GetMean() * limits[currentLimit]);
1765 if (limit < 1)
1766 limit = 1;
1767
1768 TGraph* graphMB = new TGraph;
1769 graphMB->SetTitle(Form("Events with %d times above <n> (i.e. n >= %d)", limits[currentLimit], limit));
1770 graphMB->SetMarkerStyle(20);
1771
1772 TGraph* graphBoth = new TGraph;
1773 graphBoth->SetMarkerStyle(21);
1774
1775 Float_t min = bandwidth;
1776 Float_t max = 0;
1777
1778 for (Int_t current = 0; current<nRates; ++current)
1779 {
1780 Float_t rate = rates[current];
1781 Int_t cut = cuts[current];
1782
1783 TH1* triggerEff = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
1784 TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
1785
1786 Float_t downScaleMB1 = rate / bandwidth;
1787 if (downScaleMB1 < 1)
1788 downScaleMB1 = 1;
1789
1790 Float_t downScaleMB2 = rate / (bandwidth * fractionMB);
1791 if (downScaleMB2 < 1)
1792 downScaleMB2 = 1;
1793
1794 Float_t downScaleHM = rate * proj->Integral(1, xSection->GetNbinsX()) / (bandwidth * fractionHM);
1795 if (downScaleHM < 1)
1796 downScaleHM = 1;
1797
1798 Float_t rateMB1 = rate / downScaleMB1 * xSection->Integral(limit, xSection->GetNbinsX());
1799 Float_t rateMB2 = rate / downScaleMB2 * xSection->Integral(limit, xSection->GetNbinsX());
1800 Float_t rateHM = rate / downScaleHM * proj->Integral(limit, xSection->GetNbinsX());
1801 Float_t combinedRate = rateMB2 + rateHM;
1802
1803 graphMB->SetPoint(graphMB->GetN(), rate, rateMB1);
1804 graphBoth->SetPoint(graphBoth->GetN(), rate, combinedRate);
1805
1806 min = TMath::Min(min, TMath::Min(rateMB1, combinedRate));
1807 max = TMath::Max(min, TMath::Max(rateMB1, combinedRate));
1808
1809 Printf("The rates for events with %d times above <n> (i.e. n >= %d) at a rate of %.2e Hz is:", limits[currentLimit], limit, rate);
1810 Printf(" %.2e Hz in MB-only mode", rateMB1);
1811 Printf(" %.2e Hz = %.2e Hz + %.2e Hz in MB + high mult. mode", combinedRate, rateMB2, rateHM);
1812
1813 Printf(" The downscale factors are: %.2f %.2f %.2f", downScaleMB1, downScaleMB2, downScaleHM);
1814
1815 Int_t triggerLimit = 0;
1816 for (Int_t bin = 1; bin <= triggerEff->GetNbinsX(); bin++)
1817 if (triggerEff->GetBinContent(bin) < 0.5)
1818 triggerLimit = (Int_t) triggerEff->GetXaxis()->GetBinCenter(bin);
1819
1820 Printf(" Efficiency limit (50%%) is at multiplicity %d", triggerLimit);
144ff489 1821 Float_t fractionGood = proj->Integral(triggerLimit, proj->GetNbinsX()) / proj->Integral();
1822 Printf(" %.2f %% of the events are above the trigger limit", fractionGood * 100);
d1594f8a 1823
1824 if (triggerLimit > limit)
1825 Printf(" WARNING: interesting events also counted inside the trigger limit");
1826
1827 Printf("");
1828 }
1829
1830 canvas->cd(currentLimit+1)->SetLogx();
1831 canvas->cd(currentLimit+1)->SetLogy();
1832
1833 graphMB->Draw("AP");
1834 graphBoth->Draw("P SAME");
1835
1836 graphMB->GetYaxis()->SetRangeUser(0.5 * min, 2 * max);
1837 graphMB->GetXaxis()->SetTitle("Raw rate in Hz");
1838 graphMB->GetYaxis()->SetTitle("Event rate in Hz");
1839 }
1840 }
1841}