merging post-CVS changes into SVN:
[u/mrichter/AliRoot.git] / PWG0 / highMultiplicity / AliHighMultiplicitySelector.cxx
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 #include <TMath.h>
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 //
34 // Selector that produces plots needed for the high multiplicity analysis with the
35 // pixel detector
36 // Needs cluster file
37 //
38
39 ClassImp(AliHighMultiplicitySelector)
40
41 AliHighMultiplicitySelector::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),
55   fCentralRegion(kFALSE)
56 {
57   //
58   // Constructor. Initialization of pointers
59   //
60 }
61
62 AliHighMultiplicitySelector::~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
72 void AliHighMultiplicitySelector::SlaveBegin(TTree *tree)
73 {
74   // create histograms
75
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
87 void 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
110 Bool_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
126   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, AliPWG0Helper::kMB1);
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
160     if (fCentralRegion)
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);
201   if (digits);
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   {
225     if (fCentralRegion)
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
392 Bool_t AliHighMultiplicitySelector::Notify()
393 {
394   // get next ITS runloader
395
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
411 void 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
433 void 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
456 void AliHighMultiplicitySelector::WriteHistograms(const char* filename)
457 {
458   // write the histograms to a file
459
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
471 void AliHighMultiplicitySelector::ReadHistograms(const char* filename)
472 {
473   // read the data from a file and fill histograms
474
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
533 TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut)
534 {
535   //
536   // returns the trigger efficiency as function of multiplicity with a given cut
537   //
538
539   //cut and multiply with x-section
540   TH1* allEvents = multVsLayer->ProjectionX("fMvsL_x_total", 1, 1001);
541   //allEvents->Sumw2();
542
543   //cut and multiply with x-section
544   TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, 1001);
545   //proj->Sumw2();
546
547   //new TCanvas; allEvents->DrawCopy(); gPad->SetLogy();
548   //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
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
554   return proj;
555 }
556
557 TH1* 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();
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
585   //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
586
587   return proj;
588 }
589
590 void AliHighMultiplicitySelector::MakeGraphs2(const char* title, TH1* xSection, TH2* fMvsL)
591 {
592   // creates graphs
593
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();
609
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
652 void AliHighMultiplicitySelector::MakeGraphs(const char* title, TH1* xSection, TH2* fMvsL, Int_t limit)
653 {
654   // relative x-section, once we have a collision
655
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
721 void AliHighMultiplicitySelector::JPRPlots()
722 {
723   // more plots
724
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
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
765     xSection->SetStats(kFALSE);
766     xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
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);
769     //xSection->GetYaxis()->SetTitle("relative cross section");
770     xSection->GetYaxis()->SetTitleOffset(1.2);
771
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);
807     //canvas2->SetTopMargin(0.05);
808     //canvas2->SetRightMargin(0.05);
809     canvas2->SetLogy();
810     xSection->DrawCopy("HIST");
811
812     TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.3);
813     legend2->SetFillColor(0);
814     legend2->AddEntry(xSection, "no trigger");
815
816     TH1* proj2 = (TH1*) proj->Clone("random");
817     proj2->Reset();
818     // MB lengthRun s 100 Hz
819     Int_t nTrigger = (Int_t) (100 * lengthRun * proj->Integral(1, 1000));
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
827     /*
828     proj2->DrawCopy("SAME");
829     legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio)));
830     */
831
832     proj2 = (TH1*) proj->Clone("random2");
833     proj2->Reset();
834     // 10^31 lum --> 1200 kHz; lengthRun s
835     nTrigger = (Int_t) (rate * proj->Integral(1, 1000) * lengthRun);
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);
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));
849
850     legend2->Draw();
851     line->Draw();
852
853     canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
854     canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
855   }
856 }
857
858 void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
859 {
860   //
861   // produces a spectrum created with N triggers
862   // number of triggers and thresholds for the moment fixed
863   //
864
865   /*
866
867   gSystem->Load("libANALYSIS");
868   gSystem->Load("libPWG0base");
869   .L AliHighMultiplicitySelector.cxx+g
870   x = new AliHighMultiplicitySelector();
871   x->ReadHistograms("highmult_hijing100k.root");
872   x->Ntrigger();
873
874   gSystem->Load("libPWG0base");
875   .L AliHighMultiplicitySelector.cxx+g
876   x = new AliHighMultiplicitySelector();
877   x->ReadHistograms("highmult_hijing100k.root");
878   x->Ntrigger(kFALSE);
879
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
891   // 10^28 lum --> 1.2 kHz
892   // 10^31 lum --> 1200 kHz
893   //Float_t rate = 1200e3;
894   Float_t rate = 250e3;
895
896   // time in s
897   Float_t lengthRun = 1e6;
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
911     //Int_t nCuts = 6;
912     //Int_t cuts[] = { 0, 164, 178, 190, 204, 216 };
913     //Int_t nCuts = 4;
914     //Int_t cuts[] = { 0, 164, 190, 216 };
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;
921     Int_t cuts[] = { 0, 126, 162 };
922     Float_t ratePerTrigger[] = { 60, 20.0, 20.0 };
923
924     // desired trigger rate in Hz
925     //Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 };
926
927     xSection->SetStats(kFALSE);
928     xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
929     xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5));
930     xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 450 : 350);
931     //xSection->GetYaxis()->SetTitle("relative cross section");
932     xSection->GetYaxis()->SetTitleOffset(1.2);
933
934     // relative x-section, once we have a collision
935     xSection->Scale(1.0 / xSection->Integral());
936
937     TCanvas* canvas2 = new TCanvas(Form("HMPlots2_%d_Random", i), Form("HMPlots2_%d_Random", i), 800, 600);
938     canvas2->SetTopMargin(0.05);
939     canvas2->SetRightMargin(0.05);
940     canvas2->SetLogy();
941
942     if (relative)
943       xSection->DrawCopy("HIST");
944
945     TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.4);
946     legend2->SetFillColor(0);
947     legend2->AddEntry(xSection, "cross-section");
948
949     for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut)
950     {
951       Int_t cut = cuts[currentCut];
952
953       TH1* triggerEff = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
954
955       TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
956
957       Float_t triggerLimit = 0;
958       for (Int_t bin = 1; bin <= triggerEff->GetNbinsX(); bin++)
959         if (triggerEff->GetBinContent(bin) < 0.5)
960           triggerLimit = triggerEff->GetXaxis()->GetBinCenter(bin);
961
962       Printf("Efficiency limit (50%%) is at multiplicity %f", triggerLimit);
963
964       Double_t total = 0;
965       if (proj->Integral(1, 1000) > 0)
966         total = proj->Integral(1, 1000);
967
968       printf("Cut at %d: rel. x-section = %e\n", cut, total);
969
970       TH1* proj2 = (TH1*) proj->Clone("random2");
971       proj2->Reset();
972
973       // calculate downscale factor
974       Float_t normalRate = rate * proj->Integral(1, 1000);
975       Float_t downScale = normalRate / ratePerTrigger[currentCut];
976       if (downScale < 1)
977         downScale = 1;
978       Long64_t nTrigger = (Long64_t) (normalRate / downScale * lengthRun);
979       nTrigger = TMath::Nint(((Float_t) nTrigger) / 1000) * 1000;
980
981       Printf("Normal rate is %f, downscale: %f, Simulating %lld triggers", normalRate, downScale, nTrigger);
982       proj2->FillRandom(proj, nTrigger);
983
984       Int_t removed = 0;
985       for (Int_t bin=1; bin<proj2->GetNbinsX(); ++bin)
986         if (proj2->GetBinContent(bin) < 5)
987         {
988           removed += (Int_t) proj2->GetBinContent(bin);
989           proj2->SetBinContent(bin, 0);
990         }
991
992       Printf("Removed %d events", removed);
993
994       if (relative)
995         proj2->Scale(1.0 / nTrigger * proj->Integral(1, 1000));
996
997       proj2->SetLineColor(colors[currentCut]);
998       proj2->SetMarkerStyle(markers[currentCut]);
999       proj2->SetMarkerColor(colors[currentCut]);
1000
1001       if (relative || currentCut > 0) {
1002         proj2->DrawCopy("SAME P");
1003       } else
1004         proj2->DrawCopy(" P");
1005
1006       legend2->AddEntry(proj2, Form("%lld evts, FO > %d chips", nTrigger, cut));
1007
1008       TLine* line = new TLine(triggerLimit, proj2->GetMinimum(), triggerLimit, proj2->GetMaximum());
1009       line->SetLineColor(colors[currentCut]);
1010       line->Draw();
1011     }
1012
1013     legend2->Draw();
1014
1015     canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
1016     canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
1017   }
1018 }
1019
1020 void AliHighMultiplicitySelector::Contamination()
1021 {
1022   //
1023   // produces a spectrum created with N triggers
1024   // number of triggers and thresholds for the moment fixed
1025   //
1026
1027   /*
1028
1029   gSystem->Load("libANALYSIS");
1030   gSystem->Load("libPWG0base");
1031   .L AliHighMultiplicitySelector.cxx+g
1032   x = new AliHighMultiplicitySelector();
1033   x->ReadHistograms("highmult_hijing100k.root");
1034   x->Contamination();
1035
1036   */
1037
1038   // get x-sections
1039   TFile* file = TFile::Open("crosssectionEx.root");
1040   if (!file)
1041     return;
1042
1043   TH1* xSections[2];
1044   xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1045   xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1046
1047   // rate = L * sigma
1048   // sigma ~ 80 mb (Pythia 14 TeV)
1049   // 10^28 lum --> 8e2 Hz
1050   // 10^31 lum --> 8e5 Hz
1051   Double_t rates[] = { 8e2, 8e3, 8e4, 8e5 };
1052
1053   Int_t nCuts = 4;
1054   Int_t cuts[] = { 104, 134, 154, 170 };
1055
1056   // put to 2 for second layer
1057   for (Int_t i=0; i<1; ++i)
1058   {
1059     if (!xSections[i])
1060       continue;
1061
1062     // relative x-section, once we have a collision
1063     xSections[i]->Scale(1.0 / xSections[i]->Integral());
1064
1065     Int_t max = xSections[i]->GetNbinsX();
1066     max = 500;
1067
1068     Float_t* xSection = new Float_t[max];
1069     for (Int_t mult = 0; mult < max; mult++)
1070       xSection[mult] = xSections[i]->GetBinContent(mult+1);
1071
1072     TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
1073
1074     TGraph* graph = new TGraph;
1075
1076     for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut)
1077     {
1078       Int_t cut = cuts[currentCut];
1079       Double_t rate = rates[currentCut];
1080       //Double_t rate = rates[3];
1081
1082       // coll. in 100 ns window
1083       //Double_t windowSize = 100e-9;
1084       Double_t windowSize = 25e-9;
1085       Double_t collPerWindow = windowSize * rate;
1086       Printf("coll/window = %f", collPerWindow);
1087       Double_t windowsPerSecond = 1.0 / windowSize;
1088
1089       TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
1090       Float_t* triggerEff = new Float_t[max];
1091       for (Int_t mult = 0; mult < max; mult++)
1092         triggerEff[mult] = triggerEffHist->GetBinContent(mult+1);
1093
1094       Double_t triggerRate = 0;
1095       for (Int_t mult = 0; mult < max; mult++)
1096         triggerRate += xSection[mult] * triggerEff[mult];
1097
1098       triggerRate *= TMath::Poisson(1, collPerWindow) * windowsPerSecond;
1099
1100       Printf("Rate for 1 collision is %f Hz", triggerRate);
1101
1102       Double_t triggerRate2 = 0;
1103       for (Int_t mult = 0; mult < max; mult++)
1104         for (Int_t mult2 = mult; mult2 < max; mult2++)
1105           if (mult+mult2 < max)
1106             triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2];
1107
1108       triggerRate2 *= TMath::Poisson(2, collPerWindow) * windowsPerSecond;
1109
1110       Printf("Rate for 2 collisions is %f Hz --> %.1f%%", triggerRate2, triggerRate2 / triggerRate * 100);
1111
1112       Double_t triggerRate3 = 0;
1113
1114       for (Int_t mult = 0; mult < max; mult++)
1115         for (Int_t mult2 = mult; mult2 < max-mult; mult2++)
1116           for (Int_t mult3 = 0; mult3 < max-mult-mult2; mult3++)
1117             //if (mult+mult2+mult3 < max)
1118               triggerRate3 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * xSection[mult3] * triggerEff[mult+mult2+mult3];
1119
1120       triggerRate3 *= TMath::Poisson(3, collPerWindow) * windowsPerSecond;
1121       //triggerRate3 *= collPerWindow * collPerWindow * rate;
1122
1123       Printf("Rate for 3 collisions is %f Hz --> %.1f%%", triggerRate3, triggerRate3 / triggerRate * 100);
1124
1125       Float_t totalContamination = (triggerRate2 + triggerRate3) / (triggerRate + triggerRate2 + triggerRate3);
1126
1127       Printf("Total contamination is %.1f%%", totalContamination * 100);
1128
1129       graph->SetPoint(graph->GetN(), cut, totalContamination);
1130
1131       continue;
1132
1133       Double_t triggerRate4 = 0;
1134       for (Int_t mult = 0; mult < max; mult++)
1135         for (Int_t mult2 = mult; mult2 < max-mult; mult2++)
1136           for (Int_t mult3 = 0; mult3 < max-mult-mult2; mult3++)
1137             for (Int_t mult4 = 0; mult4 < max-mult-mult2-mult3; mult4++)
1138               //if (mult+mult2+mult3+mult4 < max)
1139                 triggerRate4 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * xSection[mult3] * xSection[mult4] * triggerEff[mult+mult2+mult3+mult4];
1140
1141       //triggerRate4 *= collPerWindow * collPerWindow * collPerWindow * rate;
1142       triggerRate4 *= TMath::Poisson(4, collPerWindow) * windowsPerSecond;
1143
1144       Printf("Rate for 4 collisions is %f Hz --> %.1f%%", triggerRate4, triggerRate4 / triggerRate * 100);
1145
1146       // general code for n collisions follows, however much slower...
1147       /*
1148       const Int_t maxdepth = 4;
1149       for (Int_t depth = 1; depth <= maxdepth; depth++) {
1150         Double_t triggerRate = 0;
1151
1152         Int_t m[maxdepth];
1153         for (Int_t d=0; d<maxdepth; d++)
1154           m[d] = 0;
1155
1156         while (m[0] < max) {
1157           Double_t value = 1;
1158           Int_t sum = 0;
1159           for (Int_t d=0; d<depth; d++) {
1160             value *= xSection[m[d]];
1161             sum += m[d];
1162           }
1163
1164           if (sum < max) {
1165             value *= triggerEff[sum];
1166             triggerRate += value;
1167           }
1168
1169           Int_t increase = depth-1;
1170           ++m[increase];
1171           while (m[increase] == max && increase > 0) {
1172             m[increase] = 0;
1173             --increase;
1174             ++m[increase];
1175           }
1176         }
1177
1178         triggerRate *= rate * TMath::Power(collPerWindow, depth - 1);
1179
1180         Printf("Rate for %d collisions is %f Hz", depth, triggerRate);
1181       }*/
1182     }
1183
1184     new TCanvas; graph->Draw("AP*");
1185   }
1186 }
1187
1188 void AliHighMultiplicitySelector::Contamination2()
1189 {
1190   //
1191   // produces a spectrum created with N triggers
1192   // number of triggers and thresholds for the moment fixed
1193   //
1194
1195   /*
1196
1197   gSystem->Load("libANALYSIS");
1198   gSystem->Load("libPWG0base");
1199   .L AliHighMultiplicitySelector.cxx+g
1200   x = new AliHighMultiplicitySelector();
1201   x->ReadHistograms("highmult_hijing100k.root");
1202   x->Contamination2();
1203
1204   */
1205
1206   // get x-sections
1207   TFile* file = TFile::Open("crosssectionEx.root");
1208   if (!file)
1209     return;
1210
1211   TH1* xSections[2];
1212   xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1213   xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1214
1215   Int_t nCuts = 4;
1216   Int_t cuts[] = { 104, 134, 154, 170 };
1217
1218   new TCanvas;
1219
1220   Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
1221   Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
1222
1223   // put to 2 for second layer
1224   for (Int_t i=0; i<1; ++i)
1225   {
1226     if (!xSections[i])
1227       continue;
1228
1229     // relative x-section, once we have a collision
1230     xSections[i]->Scale(1.0 / xSections[i]->Integral());
1231
1232     Int_t max = xSections[i]->GetNbinsX();
1233     max = 500;
1234
1235     Float_t* xSection = new Float_t[max];
1236     for (Int_t mult = 0; mult < max; mult++)
1237       xSection[mult] = xSections[i]->GetBinContent(mult+1);
1238
1239     TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
1240
1241     for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut)
1242     {
1243       TGraph* graph = new TGraph;
1244       graph->SetMarkerColor(colors[currentCut]);
1245       graph->SetMarkerStyle(markers[currentCut]);
1246
1247       Int_t cut = cuts[currentCut];
1248
1249       TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
1250       Float_t* triggerEff = new Float_t[max];
1251       for (Int_t mult = 0; mult < max; mult++)
1252         triggerEff[mult] = triggerEffHist->GetBinContent(mult+1);
1253
1254       Double_t triggerRate = 0;
1255       for (Int_t mult = 0; mult < max; mult++)
1256         triggerRate += xSection[mult] * triggerEff[mult];
1257
1258       Printf("Raw value for 1 collision is %e", triggerRate);
1259
1260       Double_t triggerRate2 = 0;
1261       for (Int_t mult = 0; mult < max; mult++)
1262         for (Int_t mult2 = mult; mult2 < max; mult2++)
1263           if (mult+mult2 < max)
1264             triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2];
1265
1266       Printf("Raw value for 2 collisions is %e", triggerRate2);
1267
1268       for (Double_t doubleRate = 0; doubleRate <= 0.2; doubleRate += 0.005)
1269       {
1270         Float_t totalContamination = (triggerRate2 * doubleRate) / (triggerRate + triggerRate2 * doubleRate);
1271
1272         //Printf("Total contamination is %.1f%%", totalContamination * 100);
1273
1274         graph->SetPoint(graph->GetN(), doubleRate, totalContamination);
1275       }
1276
1277       graph->Draw((currentCut == 0) ? "A*" : "* SAME");
1278       graph->GetXaxis()->SetRangeUser(0, 1);
1279     }
1280   }
1281 }
1282
1283 void AliHighMultiplicitySelector::DrawHistograms()
1284 {
1285   // draws the histograms
1286
1287   /*
1288
1289   gSystem->Load("libPWG0base");
1290   .L AliHighMultiplicitySelector.cxx+
1291   x = new AliHighMultiplicitySelector();
1292   x->ReadHistograms("highmult_pythia.root");
1293   x->DrawHistograms();
1294
1295
1296   gSystem->Load("libPWG0base");
1297   .L AliHighMultiplicitySelector.cxx+
1298   x = new AliHighMultiplicitySelector();
1299   x->ReadHistograms("highmult_hijing.root");
1300   x->DrawHistograms();
1301
1302   gSystem->Load("libPWG0base");
1303   .L AliHighMultiplicitySelector.cxx+
1304   x = new AliHighMultiplicitySelector();
1305   x->ReadHistograms("highmult_central.root");
1306   x->DrawHistograms();
1307
1308   gSystem->Load("libPWG0base");
1309   .L AliHighMultiplicitySelector.cxx+
1310   x = new AliHighMultiplicitySelector();
1311   x->ReadHistograms("highmult_hijing100k.root");
1312   x->DrawHistograms();
1313
1314   */
1315
1316   /*TCanvas* canvas = new TCanvas("chips", "chips", 600, 400);
1317
1318   fChipsLayer2->SetLineColor(2);
1319   fChipsLayer2->SetStats(kFALSE);
1320   fChipsLayer1->SetStats(kFALSE);
1321   fChipsLayer2->SetTitle("");
1322   fChipsLayer2->DrawCopy();
1323   fChipsLayer1->DrawCopy("SAME");
1324   canvas->SaveAs("chips.gif");
1325
1326   canvas = new TCanvas("L1vsL2", "L1vsL2", 600, 400);
1327   fL1vsL2->SetStats(kFALSE);
1328   fL1vsL2->DrawCopy("COLZ");
1329   gPad->SetLogz();
1330   canvas->SaveAs("L1vsL2.gif");*/
1331
1332   TCanvas *canvas = new TCanvas("L1", "L1", 800, 600);
1333   canvas->SetTopMargin(0.05);
1334   canvas->SetRightMargin(0.12);
1335   fMvsL1->SetStats(kFALSE);
1336   fMvsL1->DrawCopy("COLZ");
1337   gPad->SetLogz();
1338
1339   canvas->SaveAs("L1NoCurve.gif");
1340   canvas->SaveAs("L1NoCurve.eps");
1341
1342   // draw corresponding theoretical curve
1343   TF1* func = new TF1("func", "[0]*(1-(1-1/[0])**x)", 1, 1000);
1344   func->SetParameter(0, 400-5*2);
1345   func->DrawCopy("SAME");
1346
1347   canvas->SaveAs("L1.gif");
1348
1349   canvas = new TCanvas("L2", "L2", 600, 400);
1350   //fMvsL2->GetYaxis()->SetRangeUser(0, 150);
1351   fMvsL2->SetStats(kFALSE);
1352   fMvsL2->DrawCopy("COLZ");
1353   gPad->SetLogz();
1354   func->SetParameter(0, 800-5*4);
1355   func->DrawCopy("SAME");
1356   canvas->SaveAs("L2.gif");
1357
1358   // get x-sections
1359   TFile* file = TFile::Open("crosssectionEx.root");
1360   if (file)
1361   {
1362     TH1* xSection2 = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1363     TH1* xSection15 = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1364
1365     MakeGraphs2("Layer1", xSection2, fMvsL1);
1366     return;
1367
1368     // 5 times the mean
1369     //MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 5)); //75 * 2 * 2);
1370     //MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 5)); //(Int_t) (75 * 1.5 * 2));
1371
1372     MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 8));
1373     MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 8));
1374     MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 9));
1375     MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 9));
1376     MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 10));
1377     MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 10));
1378
1379     file->Close();
1380   }
1381
1382   // make spread hists
1383   TGraph* spread = new TGraph;
1384   spread->SetTitle("Spread L1;true multiplicity;RMS");
1385
1386   for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
1387   {
1388     TH1* proj = fMvsL1->ProjectionY("proj", i, i);
1389     spread->SetPoint(spread->GetN(), i, proj->GetRMS());
1390   }
1391
1392   canvas = new TCanvas("SpreadL1", "SpreadL1", 600, 400);
1393   spread->Draw("A*");
1394   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1395
1396   TF1* log = new TF1("log", "[0]*log([1]*x)", 1, 150);
1397   log->SetParLimits(0, 0, 10);
1398   log->SetParLimits(1, 1e-5, 10);
1399
1400   spread->Fit(log, "", "", 1, 150);
1401   log->DrawCopy("SAME");
1402
1403   TGraph* spread2 = new TGraph;
1404   spread2->SetTitle("Spread L2;true multiplicity;RMS");
1405
1406   for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
1407   {
1408     TH1* proj = fMvsL2->ProjectionY("proj", i, i);
1409     spread2->SetPoint(spread2->GetN(), i, proj->GetRMS());
1410   }
1411
1412   canvas = new TCanvas("SpreadL2", "SpreadL2", 600, 400);
1413   spread2->Draw("A*");
1414   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1415
1416   spread2->Fit(log, "", "", 1, 150);
1417   log->DrawCopy("SAME");
1418
1419   canvas = new TCanvas("Clusters_L1", "Clusters_L1", 600, 400);
1420   fClvsL1->SetStats(kFALSE);
1421   fClvsL1->DrawCopy("COLZ");
1422   gPad->SetLogz();
1423
1424   func->SetParameter(0, 400-5*2);
1425   func->DrawCopy("SAME");
1426
1427   canvas->SaveAs("Clusters_L1.gif");
1428
1429   canvas = new TCanvas("Clusters_L2", "Clusters_L2", 600, 400);
1430   fClvsL2->SetStats(kFALSE);
1431   fClvsL2->DrawCopy("COLZ");
1432   gPad->SetLogz();
1433   func->SetParameter(0, 800-5*4);
1434   func->DrawCopy("SAME");
1435   canvas->SaveAs("Clusters_L2.gif");
1436
1437   canvas = new TCanvas("ChipsFired", "ChipsFired", 600, 400);
1438   //fChipsFired->GetYaxis()->SetRangeUser(0, 150);
1439   fChipsFired->SetStats(kFALSE);
1440   fChipsFired->DrawCopy("COLZ");
1441   canvas->SaveAs("ChipsFired.gif");
1442
1443   /*TH1F* tresholdHistL1 = new TH1F("tresholdHistL1", ";chip treshold;<n>", BINNING_LAYER1);
1444   TH1F* tresholdHistL2 = new TH1F("tresholdHistL2", ";chip treshold;<n>", BINNING_LAYER2);
1445
1446   for (Int_t treshold = 0; treshold < 800; treshold++)
1447   {
1448     if (fPrimaryL1->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
1449     {
1450       TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL1->GetHistogram());
1451       if (mult)
1452         tresholdHistL1->Fill(treshold, mult->GetMean());
1453     }
1454     if (fPrimaryL2->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
1455     {
1456       TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL2->GetHistogram());
1457       if (mult)
1458         tresholdHistL2->Fill(treshold, mult->GetMean());
1459     }
1460   }
1461
1462   canvas = new TCanvas("TresholdL1", "TresholdL1", 600, 400);
1463   tresholdHistL1->Draw();
1464   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1465
1466   canvas = new TCanvas("TresholdL2", "TresholdL2", 600, 400);
1467   tresholdHistL2->Draw();
1468   canvas->SaveAs(Form("%s.gif", canvas->GetName()));*/
1469
1470   fPrimaryL1->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff1", "", "prof goff");
1471   fPrimaryL2->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff2", "", "prof goff");
1472
1473   canvas = new TCanvas("Efficiency", "Efficiency", 600, 400);
1474   fPrimaryL1->GetHistogram()->SetStats(kFALSE);
1475   fPrimaryL1->GetHistogram()->Draw();
1476   fPrimaryL2->GetHistogram()->SetLineColor(2);
1477   fPrimaryL2->GetHistogram()->SetStats(kFALSE);
1478   fPrimaryL2->GetHistogram()->Draw("SAME");
1479   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1480
1481   canvas = new TCanvas("ClustersZL1", "ClustersZL1", 600, 400);
1482   fClusterZL1->Rebin(2);
1483   fClusterZL1->Draw();
1484   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1485
1486   canvas = new TCanvas("ClustersZL2", "ClustersZL2", 600, 400);
1487   fClusterZL2->Draw();
1488   fClusterZL2->Rebin(2);
1489   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1490 }
1491
1492 TGraph* AliHighMultiplicitySelector::IntFractRate()
1493 {
1494   // A plot which shows the fractional rate above any threshold
1495   // as function of threshold (i.e. the integral of dSigma/dN as function of
1496   // N, normalised to 1 for N=0)
1497
1498   /*
1499   gSystem->Load("libANALYSIS");
1500   gSystem->Load("libPWG0base");
1501   .L AliHighMultiplicitySelector.cxx+
1502   x = new AliHighMultiplicitySelector();
1503   x->ReadHistograms("highmult_hijing100k.root");
1504   x->IntFractRate();
1505   */
1506
1507   // get x-sections
1508   TFile* file = TFile::Open("crosssectionEx.root");
1509   if (!file)
1510     return 0;
1511
1512   TH1* xSection;
1513   xSection = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1514
1515   TGraph* result = new TGraph;
1516
1517   for (Int_t threshold = 0; threshold < 300; threshold += 2)
1518   {
1519     TH1* proj = GetXSectionCut(xSection, fMvsL1, threshold);
1520
1521     //new TCanvas; proj->DrawCopy();
1522
1523     Double_t integral = proj->Integral();
1524
1525     Printf("Cut at %d, integral is %e", threshold, integral);
1526
1527     result->SetPoint(result->GetN(), threshold, integral);
1528   }
1529
1530   TCanvas* canvas = new TCanvas("IntFractRate", "IntFractRate", 600, 400);
1531   gPad->SetLogy();
1532
1533   result->Draw("A*");
1534   result->GetXaxis()->SetTitle("threshold");
1535   result->GetYaxis()->SetTitle("integrated fractional rate above threshold");
1536
1537   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1538
1539   return result;
1540 }
1541
1542 void AliHighMultiplicitySelector::MBComparison()
1543 {
1544   //
1545   // finds the threshold from which onwards the number of found events above N times the mean
1546   // is higher using a high mult. trigger than just triggering with MB
1547   //
1548
1549   /*
1550
1551   gSystem->Load("libANALYSIS");
1552   gSystem->Load("libPWG0base");
1553   .L AliHighMultiplicitySelector.cxx+g
1554   x = new AliHighMultiplicitySelector();
1555   x->ReadHistograms("highmult_hijing100k.root");
1556   x->MBComparison();
1557
1558   */
1559
1560   // get x-sections
1561   TFile* file = TFile::Open("crosssectionEx.root");
1562   if (!file)
1563     return;
1564
1565   TH1* xSections[2];
1566   xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1567   xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1568
1569   // rate = L * sigma
1570   // sigma ~ 80 mb (Pythia 14 TeV)
1571   // 10^28 lum --> 8e2 Hz
1572   // 10^31 lum --> 8e5 Hz
1573   Int_t nRates = 4;
1574   Double_t rates[] = { 8e2, 8e3, 8e4, 8e5 };
1575
1576   // threshold in number of fired chips for corresponding rate
1577   //Int_t cuts[] = { 104, 134, 154, 170 }; // values for 20 Hz
1578   Int_t cuts[] = { 82, 124, 147, 164 };    // values for 50 Hz
1579
1580   // bandwidth, fractions (for MB, high mult.)
1581   Float_t bandwidth = 1e3;
1582   Float_t fractionMB = 0.5;
1583   Float_t fractionHM = 0.05;
1584
1585   // different limits to define "interesting events"
1586   Int_t nLimits = 9;
1587   Int_t limits[] = { 0, 1, 2, 4, 6, 7, 8, 9, 10 };
1588
1589   // put to 2 for second layer
1590   for (Int_t i=0; i<1; ++i)
1591   {
1592     if (!xSections[i])
1593       continue;
1594
1595     TH1* xSection = xSections[i];
1596     TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
1597
1598     // relative x-section, once we have a collision
1599     xSection->Scale(1.0 / xSection->Integral());
1600
1601     xSection->SetStats(kFALSE);
1602     xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
1603     xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5));
1604     xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 450 : 350);
1605     xSection->GetYaxis()->SetTitleOffset(1.2);
1606
1607     TCanvas* canvas = new TCanvas("MBComparison", "MBComparison", 1000, 800);
1608     canvas->Divide(3, 3);
1609
1610     for (Int_t currentLimit = 0; currentLimit<nLimits; currentLimit++)
1611     {
1612       // limit is N times the mean
1613       Int_t limit = (Int_t) (xSection->GetMean() * limits[currentLimit]);
1614       if (limit < 1)
1615         limit = 1;
1616
1617       TGraph* graphMB = new TGraph;
1618       graphMB->SetTitle(Form("Events with %d times above <n> (i.e. n >= %d)", limits[currentLimit], limit));
1619       graphMB->SetMarkerStyle(20);
1620
1621       TGraph* graphBoth = new TGraph;
1622       graphBoth->SetMarkerStyle(21);
1623
1624       Float_t min = bandwidth;
1625       Float_t max = 0;
1626
1627       for (Int_t current = 0; current<nRates; ++current)
1628       {
1629         Float_t rate = rates[current];
1630         Int_t cut = cuts[current];
1631
1632         TH1* triggerEff = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
1633         TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
1634
1635         Float_t downScaleMB1 = rate / bandwidth;
1636         if (downScaleMB1 < 1)
1637           downScaleMB1 = 1;
1638
1639         Float_t downScaleMB2 = rate / (bandwidth * fractionMB);
1640         if (downScaleMB2 < 1)
1641           downScaleMB2 = 1;
1642
1643         Float_t downScaleHM = rate * proj->Integral(1, xSection->GetNbinsX()) / (bandwidth * fractionHM);
1644         if (downScaleHM < 1)
1645           downScaleHM = 1;
1646
1647         Float_t rateMB1 = rate / downScaleMB1 * xSection->Integral(limit, xSection->GetNbinsX());
1648         Float_t rateMB2 = rate / downScaleMB2 * xSection->Integral(limit, xSection->GetNbinsX());
1649         Float_t rateHM = rate / downScaleHM * proj->Integral(limit, xSection->GetNbinsX());
1650         Float_t combinedRate = rateMB2 + rateHM;
1651
1652         graphMB->SetPoint(graphMB->GetN(), rate, rateMB1);
1653         graphBoth->SetPoint(graphBoth->GetN(), rate, combinedRate);
1654
1655         min = TMath::Min(min, TMath::Min(rateMB1, combinedRate));
1656         max = TMath::Max(min, TMath::Max(rateMB1, combinedRate));
1657
1658         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);
1659         Printf("   %.2e Hz in MB-only mode", rateMB1);
1660         Printf("   %.2e Hz = %.2e Hz + %.2e Hz in MB + high mult. mode", combinedRate, rateMB2, rateHM);
1661
1662         Printf("   The downscale factors are: %.2f %.2f %.2f", downScaleMB1, downScaleMB2, downScaleHM);
1663
1664         Int_t triggerLimit = 0;
1665         for (Int_t bin = 1; bin <= triggerEff->GetNbinsX(); bin++)
1666           if (triggerEff->GetBinContent(bin) < 0.5)
1667             triggerLimit = (Int_t) triggerEff->GetXaxis()->GetBinCenter(bin);
1668
1669         Printf("   Efficiency limit (50%%) is at multiplicity %d", triggerLimit);
1670         Float_t fractionGood = proj->Integral(triggerLimit, proj->GetNbinsX()) / proj->Integral();
1671         Printf("   %.2f %% of the events are above the trigger limit", fractionGood * 100);
1672
1673         if (triggerLimit > limit)
1674           Printf("   WARNING: interesting events also counted inside the trigger limit");
1675
1676         Printf("");
1677       }
1678
1679       canvas->cd(currentLimit+1)->SetLogx();
1680       canvas->cd(currentLimit+1)->SetLogy();
1681
1682       graphMB->Draw("AP");
1683       graphBoth->Draw("P SAME");
1684
1685       graphMB->GetYaxis()->SetRangeUser(0.5 * min, 2 * max);
1686       graphMB->GetXaxis()->SetTitle("Raw rate in Hz");
1687       graphMB->GetYaxis()->SetTitle("Event rate in Hz");
1688     }
1689   }
1690 }