1ecd9c1ac4a8609c533b04d4fa2b68ab1c902d34
[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
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 //
33 // Selector that produces plots needed for the high multiplicity analysis with the
34 // pixel detector
35 // Needs cluster file
36 //
37
38 ClassImp(AliHighMultiplicitySelector)
39
40 AliHighMultiplicitySelector::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),
54   fCentralRegion(kFALSE)
55 {
56   //
57   // Constructor. Initialization of pointers
58   //
59 }
60
61 AliHighMultiplicitySelector::~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
71 void AliHighMultiplicitySelector::SlaveBegin(TTree *tree)
72 {
73   // create histograms
74
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
86 void 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
109 Bool_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
159     if (fCentralRegion)
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   {
224     if (fCentralRegion)
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
391 Bool_t AliHighMultiplicitySelector::Notify()
392 {
393   // get next ITS runloader
394
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
410 void 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
432 void 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
455 void AliHighMultiplicitySelector::WriteHistograms(const char* filename)
456 {
457   // write the histograms to a file
458
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
470 void AliHighMultiplicitySelector::ReadHistograms(const char* filename)
471 {
472   // read the data from a file and fill histograms
473
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
532 TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut)
533 {
534   //
535   // returns the trigger efficiency as function of multiplicity with a given cut
536   //
537
538   //cut and multiply with x-section
539   TH1* allEvents = multVsLayer->ProjectionX("fMvsL_x_total", 1, 1001);
540   //allEvents->Sumw2();
541
542   //cut and multiply with x-section
543   TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, 1001);
544   //proj->Sumw2();
545
546   //new TCanvas; allEvents->DrawCopy(); gPad->SetLogy();
547   //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
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
553   return proj;
554 }
555
556 TH1* 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();
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
584   //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
585
586   return proj;
587 }
588
589 void AliHighMultiplicitySelector::MakeGraphs2(const char* title, TH1* xSection, TH2* fMvsL)
590 {
591   // creates graphs
592
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();
608
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
651 void AliHighMultiplicitySelector::MakeGraphs(const char* title, TH1* xSection, TH2* fMvsL, Int_t limit)
652 {
653   // relative x-section, once we have a collision
654
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
720 void AliHighMultiplicitySelector::JPRPlots()
721 {
722   // more plots
723
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
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
764     xSection->SetStats(kFALSE);
765     xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
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);
768     //xSection->GetYaxis()->SetTitle("relative cross section");
769     xSection->GetYaxis()->SetTitleOffset(1.2);
770
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);
806     //canvas2->SetTopMargin(0.05);
807     //canvas2->SetRightMargin(0.05);
808     canvas2->SetLogy();
809     xSection->DrawCopy("HIST");
810
811     TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.3);
812     legend2->SetFillColor(0);
813     legend2->AddEntry(xSection, "no trigger");
814
815     TH1* proj2 = (TH1*) proj->Clone("random");
816     proj2->Reset();
817     // MB lengthRun s 100 Hz
818     Int_t nTrigger = (Int_t) (100 * lengthRun * proj->Integral(1, 1000));
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
826     /*
827     proj2->DrawCopy("SAME");
828     legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio)));
829     */
830
831     proj2 = (TH1*) proj->Clone("random2");
832     proj2->Reset();
833     // 10^31 lum --> 1200 kHz; lengthRun s
834     nTrigger = (Int_t) (rate * proj->Integral(1, 1000) * lengthRun);
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);
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));
848
849     legend2->Draw();
850     line->Draw();
851
852     canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
853     canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
854   }
855 }
856
857 void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
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
872   gSystem->Load("libPWG0base");
873   .L AliHighMultiplicitySelector.cxx+g
874   x = new AliHighMultiplicitySelector();
875   x->ReadHistograms("highmult_hijing100k.root");
876   x->Ntrigger(kFALSE);
877
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;
892   Float_t rate = 250e3;
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
909     //Int_t nCuts = 6;
910     //Int_t cuts[] = { 0, 164, 178, 190, 204, 216 };
911     //Int_t nCuts = 4;
912     //Int_t cuts[] = { 0, 164, 190, 216 };
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 };
921
922     // desired trigger rate in Hz
923     //Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 };
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();
939
940     if (relative)
941       xSection->DrawCopy("HIST");
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
951       TH1* triggerEff = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
952
953       TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
954
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
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
991       if (relative)
992         proj2->Scale(1.0 / nTrigger * proj->Integral(1, 1000));
993
994       proj2->SetLineColor(colors[currentCut]);
995       proj2->SetMarkerStyle(markers[currentCut]);
996       proj2->SetMarkerColor(colors[currentCut]);
997
998       if (relative || currentCut > 0) {
999         proj2->DrawCopy("SAME P");
1000       } else
1001         proj2->DrawCopy(" P");
1002
1003       legend2->AddEntry(proj2, Form("%lld evts, FO > %d chips", nTrigger, cut));
1004
1005       TLine* line = new TLine(triggerLimit, proj2->GetMinimum(), triggerLimit, proj2->GetMaximum());
1006       line->SetLineColor(colors[currentCut]);
1007       line->Draw();
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
1017 void 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;
1106       //triggerRate2 *= collPerWindow * rate;
1107
1108       Printf("Rate for 2 collisions is %f Hz --> %.1f%%", triggerRate2, triggerRate2 / triggerRate * 100);
1109
1110       Double_t triggerRate3 = 0;
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
1185 void AliHighMultiplicitySelector::DrawHistograms()
1186 {
1187   // draws the histograms
1188
1189   /*
1190
1191   gSystem->Load("libPWG0base");
1192   .L AliHighMultiplicitySelector.cxx+
1193   x = new AliHighMultiplicitySelector();
1194   x->ReadHistograms("highmult_pythia.root");
1195   x->DrawHistograms();
1196
1197
1198   gSystem->Load("libPWG0base");
1199   .L AliHighMultiplicitySelector.cxx+
1200   x = new AliHighMultiplicitySelector();
1201   x->ReadHistograms("highmult_hijing.root");
1202   x->DrawHistograms();
1203
1204   gSystem->Load("libPWG0base");
1205   .L AliHighMultiplicitySelector.cxx+
1206   x = new AliHighMultiplicitySelector();
1207   x->ReadHistograms("highmult_central.root");
1208   x->DrawHistograms();
1209
1210   gSystem->Load("libPWG0base");
1211   .L AliHighMultiplicitySelector.cxx+
1212   x = new AliHighMultiplicitySelector();
1213   x->ReadHistograms("highmult_hijing100k.root");
1214   x->DrawHistograms();
1215
1216   */
1217
1218   /*TCanvas* canvas = new TCanvas("chips", "chips", 600, 400);
1219
1220   fChipsLayer2->SetLineColor(2);
1221   fChipsLayer2->SetStats(kFALSE);
1222   fChipsLayer1->SetStats(kFALSE);
1223   fChipsLayer2->SetTitle("");
1224   fChipsLayer2->DrawCopy();
1225   fChipsLayer1->DrawCopy("SAME");
1226   canvas->SaveAs("chips.gif");
1227
1228   canvas = new TCanvas("L1vsL2", "L1vsL2", 600, 400);
1229   fL1vsL2->SetStats(kFALSE);
1230   fL1vsL2->DrawCopy("COLZ");
1231   gPad->SetLogz();
1232   canvas->SaveAs("L1vsL2.gif");*/
1233
1234   TCanvas *canvas = new TCanvas("L1", "L1", 800, 600);
1235   canvas->SetTopMargin(0.05);
1236   canvas->SetRightMargin(0.12);
1237   fMvsL1->SetStats(kFALSE);
1238   fMvsL1->DrawCopy("COLZ");
1239   gPad->SetLogz();
1240
1241   canvas->SaveAs("L1NoCurve.gif");
1242   canvas->SaveAs("L1NoCurve.eps");
1243
1244   // draw corresponding theoretical curve
1245   TF1* func = new TF1("func", "[0]*(1-(1-1/[0])**x)", 1, 1000);
1246   func->SetParameter(0, 400-5*2);
1247   func->DrawCopy("SAME");
1248
1249   canvas->SaveAs("L1.gif");
1250
1251   canvas = new TCanvas("L2", "L2", 600, 400);
1252   //fMvsL2->GetYaxis()->SetRangeUser(0, 150);
1253   fMvsL2->SetStats(kFALSE);
1254   fMvsL2->DrawCopy("COLZ");
1255   gPad->SetLogz();
1256   func->SetParameter(0, 800-5*4);
1257   func->DrawCopy("SAME");
1258   canvas->SaveAs("L2.gif");
1259
1260   // get x-sections
1261   TFile* file = TFile::Open("crosssectionEx.root");
1262   if (file)
1263   {
1264     TH1* xSection2 = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1265     TH1* xSection15 = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1266
1267     MakeGraphs2("Layer1", xSection2, fMvsL1);
1268     return;
1269
1270     // 5 times the mean
1271     //MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 5)); //75 * 2 * 2);
1272     //MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 5)); //(Int_t) (75 * 1.5 * 2));
1273
1274     MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 8));
1275     MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 8));
1276     MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 9));
1277     MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 9));
1278     MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 10));
1279     MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 10));
1280
1281     file->Close();
1282   }
1283
1284   // make spread hists
1285   TGraph* spread = new TGraph;
1286   spread->SetTitle("Spread L1;true multiplicity;RMS");
1287
1288   for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
1289   {
1290     TH1* proj = fMvsL1->ProjectionY("proj", i, i);
1291     spread->SetPoint(spread->GetN(), i, proj->GetRMS());
1292   }
1293
1294   canvas = new TCanvas("SpreadL1", "SpreadL1", 600, 400);
1295   spread->Draw("A*");
1296   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1297
1298   TF1* log = new TF1("log", "[0]*log([1]*x)", 1, 150);
1299   log->SetParLimits(0, 0, 10);
1300   log->SetParLimits(1, 1e-5, 10);
1301
1302   spread->Fit(log, "", "", 1, 150);
1303   log->DrawCopy("SAME");
1304
1305   TGraph* spread2 = new TGraph;
1306   spread2->SetTitle("Spread L2;true multiplicity;RMS");
1307
1308   for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
1309   {
1310     TH1* proj = fMvsL2->ProjectionY("proj", i, i);
1311     spread2->SetPoint(spread2->GetN(), i, proj->GetRMS());
1312   }
1313
1314   canvas = new TCanvas("SpreadL2", "SpreadL2", 600, 400);
1315   spread2->Draw("A*");
1316   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1317
1318   spread2->Fit(log, "", "", 1, 150);
1319   log->DrawCopy("SAME");
1320
1321   canvas = new TCanvas("Clusters_L1", "Clusters_L1", 600, 400);
1322   fClvsL1->SetStats(kFALSE);
1323   fClvsL1->DrawCopy("COLZ");
1324   gPad->SetLogz();
1325
1326   func->SetParameter(0, 400-5*2);
1327   func->DrawCopy("SAME");
1328
1329   canvas->SaveAs("Clusters_L1.gif");
1330
1331   canvas = new TCanvas("Clusters_L2", "Clusters_L2", 600, 400);
1332   fClvsL2->SetStats(kFALSE);
1333   fClvsL2->DrawCopy("COLZ");
1334   gPad->SetLogz();
1335   func->SetParameter(0, 800-5*4);
1336   func->DrawCopy("SAME");
1337   canvas->SaveAs("Clusters_L2.gif");
1338
1339   canvas = new TCanvas("ChipsFired", "ChipsFired", 600, 400);
1340   //fChipsFired->GetYaxis()->SetRangeUser(0, 150);
1341   fChipsFired->SetStats(kFALSE);
1342   fChipsFired->DrawCopy("COLZ");
1343   canvas->SaveAs("ChipsFired.gif");
1344
1345   /*TH1F* tresholdHistL1 = new TH1F("tresholdHistL1", ";chip treshold;<n>", BINNING_LAYER1);
1346   TH1F* tresholdHistL2 = new TH1F("tresholdHistL2", ";chip treshold;<n>", BINNING_LAYER2);
1347
1348   for (Int_t treshold = 0; treshold < 800; treshold++)
1349   {
1350     if (fPrimaryL1->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
1351     {
1352       TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL1->GetHistogram());
1353       if (mult)
1354         tresholdHistL1->Fill(treshold, mult->GetMean());
1355     }
1356     if (fPrimaryL2->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
1357     {
1358       TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL2->GetHistogram());
1359       if (mult)
1360         tresholdHistL2->Fill(treshold, mult->GetMean());
1361     }
1362   }
1363
1364   canvas = new TCanvas("TresholdL1", "TresholdL1", 600, 400);
1365   tresholdHistL1->Draw();
1366   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1367
1368   canvas = new TCanvas("TresholdL2", "TresholdL2", 600, 400);
1369   tresholdHistL2->Draw();
1370   canvas->SaveAs(Form("%s.gif", canvas->GetName()));*/
1371
1372   fPrimaryL1->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff1", "", "prof goff");
1373   fPrimaryL2->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff2", "", "prof goff");
1374
1375   canvas = new TCanvas("Efficiency", "Efficiency", 600, 400);
1376   fPrimaryL1->GetHistogram()->SetStats(kFALSE);
1377   fPrimaryL1->GetHistogram()->Draw();
1378   fPrimaryL2->GetHistogram()->SetLineColor(2);
1379   fPrimaryL2->GetHistogram()->SetStats(kFALSE);
1380   fPrimaryL2->GetHistogram()->Draw("SAME");
1381   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1382
1383   canvas = new TCanvas("ClustersZL1", "ClustersZL1", 600, 400);
1384   fClusterZL1->Rebin(2);
1385   fClusterZL1->Draw();
1386   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1387
1388   canvas = new TCanvas("ClustersZL2", "ClustersZL2", 600, 400);
1389   fClusterZL2->Draw();
1390   fClusterZL2->Rebin(2);
1391   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1392 }
1393
1394 TGraph* AliHighMultiplicitySelector::IntFractRate()
1395 {
1396   // A plot which shows the fractional rate above any threshold
1397   // as function of threshold (i.e. the integral of dSigma/dN as function of
1398   // N, normalised to 1 for N=0)
1399
1400   /*
1401   gSystem->Load("libANALYSIS");
1402   gSystem->Load("libPWG0base");
1403   .L AliHighMultiplicitySelector.cxx+
1404   x = new AliHighMultiplicitySelector();
1405   x->ReadHistograms("highmult_hijing100k.root");
1406   x->IntFractRate();
1407   */
1408
1409   // get x-sections
1410   TFile* file = TFile::Open("crosssectionEx.root");
1411   if (!file)
1412     return 0;
1413
1414   TH1* xSection;
1415   xSection = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1416
1417   TGraph* result = new TGraph;
1418
1419   for (Int_t threshold = 0; threshold < 300; threshold += 2)
1420   {
1421     TH1* proj = GetXSectionCut(xSection, fMvsL1, threshold);
1422
1423     //new TCanvas; proj->DrawCopy();
1424
1425     Double_t integral = proj->Integral();
1426
1427     Printf("Cut at %d, integral is %e", threshold, integral);
1428
1429     result->SetPoint(result->GetN(), threshold, integral);
1430   }
1431
1432   TCanvas* canvas = new TCanvas("IntFractRate", "IntFractRate", 600, 400);
1433   gPad->SetLogy();
1434
1435   result->Draw("A*");
1436   result->GetXaxis()->SetTitle("threshold");
1437   result->GetYaxis()->SetTitle("integrated fractional rate above threshold");
1438
1439   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1440
1441   return result;
1442 }
1443
1444 void AliHighMultiplicitySelector::MBComparison()
1445 {
1446   //
1447   // finds the threshold from which onwards the number of found events above N times the mean
1448   // is higher using a high mult. trigger than just triggering with MB
1449   //
1450
1451   /*
1452
1453   gSystem->Load("libANALYSIS");
1454   gSystem->Load("libPWG0base");
1455   .L AliHighMultiplicitySelector.cxx+g
1456   x = new AliHighMultiplicitySelector();
1457   x->ReadHistograms("highmult_hijing100k.root");
1458   x->MBComparison();
1459
1460   */
1461
1462   // get x-sections
1463   TFile* file = TFile::Open("crosssectionEx.root");
1464   if (!file)
1465     return;
1466
1467   TH1* xSections[2];
1468   xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1469   xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1470
1471   // rate = L * sigma
1472   // sigma ~ 80 mb (Pythia 14 TeV)
1473   // 10^28 lum --> 8e2 Hz
1474   // 10^31 lum --> 8e5 Hz
1475   Int_t nRates = 4;
1476   Double_t rates[] = { 8e2, 8e3, 8e4, 8e5 };
1477
1478   // threshold in number of fired chips for corresponding rate
1479   //Int_t cuts[] = { 104, 134, 154, 170 }; // values for 20 Hz
1480   Int_t cuts[] = { 82, 124, 147, 164 };    // values for 50 Hz
1481
1482   // bandwidth, fractions (for MB, high mult.)
1483   Float_t bandwidth = 1e3;
1484   Float_t fractionMB = 0.5;
1485   Float_t fractionHM = 0.05;
1486
1487   // different limits to define "interesting events"
1488   Int_t nLimits = 9;
1489   Int_t limits[] = { 0, 1, 2, 4, 6, 7, 8, 9, 10 };
1490
1491   // put to 2 for second layer
1492   for (Int_t i=0; i<1; ++i)
1493   {
1494     if (!xSections[i])
1495       continue;
1496
1497     TH1* xSection = xSections[i];
1498     TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
1499
1500     // relative x-section, once we have a collision
1501     xSection->Scale(1.0 / xSection->Integral());
1502
1503     xSection->SetStats(kFALSE);
1504     xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
1505     xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5));
1506     xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 450 : 350);
1507     xSection->GetYaxis()->SetTitleOffset(1.2);
1508
1509     TCanvas* canvas = new TCanvas("MBComparison", "MBComparison", 1000, 800);
1510     canvas->Divide(3, 3);
1511
1512     for (Int_t currentLimit = 0; currentLimit<nLimits; currentLimit++)
1513     {
1514       // limit is N times the mean
1515       Int_t limit = (Int_t) (xSection->GetMean() * limits[currentLimit]);
1516       if (limit < 1)
1517         limit = 1;
1518
1519       TGraph* graphMB = new TGraph;
1520       graphMB->SetTitle(Form("Events with %d times above <n> (i.e. n >= %d)", limits[currentLimit], limit));
1521       graphMB->SetMarkerStyle(20);
1522
1523       TGraph* graphBoth = new TGraph;
1524       graphBoth->SetMarkerStyle(21);
1525
1526       Float_t min = bandwidth;
1527       Float_t max = 0;
1528
1529       for (Int_t current = 0; current<nRates; ++current)
1530       {
1531         Float_t rate = rates[current];
1532         Int_t cut = cuts[current];
1533
1534         TH1* triggerEff = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
1535         TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
1536
1537         Float_t downScaleMB1 = rate / bandwidth;
1538         if (downScaleMB1 < 1)
1539           downScaleMB1 = 1;
1540
1541         Float_t downScaleMB2 = rate / (bandwidth * fractionMB);
1542         if (downScaleMB2 < 1)
1543           downScaleMB2 = 1;
1544
1545         Float_t downScaleHM = rate * proj->Integral(1, xSection->GetNbinsX()) / (bandwidth * fractionHM);
1546         if (downScaleHM < 1)
1547           downScaleHM = 1;
1548
1549         Float_t rateMB1 = rate / downScaleMB1 * xSection->Integral(limit, xSection->GetNbinsX());
1550         Float_t rateMB2 = rate / downScaleMB2 * xSection->Integral(limit, xSection->GetNbinsX());
1551         Float_t rateHM = rate / downScaleHM * proj->Integral(limit, xSection->GetNbinsX());
1552         Float_t combinedRate = rateMB2 + rateHM;
1553
1554         graphMB->SetPoint(graphMB->GetN(), rate, rateMB1);
1555         graphBoth->SetPoint(graphBoth->GetN(), rate, combinedRate);
1556
1557         min = TMath::Min(min, TMath::Min(rateMB1, combinedRate));
1558         max = TMath::Max(min, TMath::Max(rateMB1, combinedRate));
1559
1560         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);
1561         Printf("   %.2e Hz in MB-only mode", rateMB1);
1562         Printf("   %.2e Hz = %.2e Hz + %.2e Hz in MB + high mult. mode", combinedRate, rateMB2, rateHM);
1563
1564         Printf("   The downscale factors are: %.2f %.2f %.2f", downScaleMB1, downScaleMB2, downScaleHM);
1565
1566         Int_t triggerLimit = 0;
1567         for (Int_t bin = 1; bin <= triggerEff->GetNbinsX(); bin++)
1568           if (triggerEff->GetBinContent(bin) < 0.5)
1569             triggerLimit = (Int_t) triggerEff->GetXaxis()->GetBinCenter(bin);
1570
1571         Printf("   Efficiency limit (50%%) is at multiplicity %d", triggerLimit);
1572
1573         if (triggerLimit > limit)
1574           Printf("   WARNING: interesting events also counted inside the trigger limit");
1575
1576         Printf("");
1577       }
1578
1579       canvas->cd(currentLimit+1)->SetLogx();
1580       canvas->cd(currentLimit+1)->SetLogy();
1581
1582       graphMB->Draw("AP");
1583       graphBoth->Draw("P SAME");
1584
1585       graphMB->GetYaxis()->SetRangeUser(0.5 * min, 2 * max);
1586       graphMB->GetXaxis()->SetTitle("Raw rate in Hz");
1587       graphMB->GetYaxis()->SetTitle("Event rate in Hz");
1588     }
1589   }
1590 }