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