c6a4c3ab69d3774844ae959d3ccc64359eccb5a7
[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 //
34
35 ClassImp(AliHighMultiplicitySelector)
36
37 AliHighMultiplicitySelector::AliHighMultiplicitySelector() :
38   AliSelectorRL(),
39   fChipsLayer1(0),
40   fChipsLayer2(0),
41   fL1vsL2(0),
42   fMvsL1(0),
43   fMvsL2(0),
44   fChipsFired(0),
45   fPrimaryL1(0),
46   fPrimaryL2(0),
47   fClusterZL1(0),
48   fClusterZL2(0),
49   fClvsL1(0),
50   fClvsL2(0),
51   centralRegion(kFALSE)
52 {
53   //
54   // Constructor. Initialization of pointers
55   //
56 }
57
58 AliHighMultiplicitySelector::~AliHighMultiplicitySelector()
59 {
60   //
61   // Destructor
62   //
63
64   // histograms are in the output list and deleted when the output
65   // list is deleted by the TSelector dtor
66 }
67
68 void AliHighMultiplicitySelector::SlaveBegin(TTree *tree)
69 {
70   AliSelectorRL::SlaveBegin(tree);
71
72   fChipsFired = new TH2F("fChipsFired", ";Module;Chip;Count", 240, -0.5, 239.5, 5, -0.5, 4.5);
73
74   fPrimaryL1 = new TNtuple("fPrimaryL1", "", "multiplicity:firedChips:chipsByPrimaries:clusters");
75   fPrimaryL2 = new TNtuple("fPrimaryL2", "", "multiplicity:firedChips:chipsByPrimaries:clusters");
76
77   fClusterZL1 = new TH1F("fClusterZL1", ";z", 400, -20, 20);
78   fClusterZL2 = new TH1F("fClusterZL2", ";z", 400, -20, 20);
79 }
80
81 void AliHighMultiplicitySelector::Init(TTree* tree)
82 {
83   // read the user objects
84
85   AliSelectorRL::Init(tree);
86
87   // enable only the needed branches
88   if (tree)
89   {
90     tree->SetBranchStatus("*", 0);
91     tree->SetBranchStatus("fTriggerMask", 1);
92
93     /*if (fTree->GetCurrentFile())
94     {
95       TString fileName(fTree->GetCurrentFile()->GetName());
96       fileName.ReplaceAll("AliESDs", "geometry");
97
98       // load geometry
99       TGeoManager::Import(fileName);
100       }*/
101   }
102 }
103
104 Bool_t AliHighMultiplicitySelector::Process(Long64_t entry)
105 {
106   //
107   // processing
108   //
109
110   if (AliSelectorRL::Process(entry) == kFALSE)
111     return kFALSE;
112
113   // Check prerequisites
114   if (!fESD)
115   {
116     AliDebug(AliLog::kError, "ESD branch not available");
117     return kFALSE;
118   }
119
120   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, AliPWG0Helper::kMB1);
121
122   if (!eventTriggered)
123   {
124     AliDebug(AliLog::kDebug, "Event not triggered. Skipping.");
125     return kTRUE;
126   }
127
128   AliStack* stack = GetStack();
129   if (!stack)
130   {
131     AliDebug(AliLog::kError, "Stack not available");
132     return kFALSE;
133   }
134
135   Int_t nPrim  = stack->GetNprimary();
136   Int_t multiplicity21 = 0;
137   Int_t multiplicity16 = 0;
138
139   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
140   {
141     AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
142
143     TParticle* particle = stack->Particle(iMc);
144
145     if (!particle)
146     {
147       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
148       continue;
149     }
150
151     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
152       continue;
153
154     if (centralRegion)
155     {
156       if (TMath::Abs(particle->Eta()) < 1.05)
157         multiplicity21++;
158       if (TMath::Abs(particle->Eta()) < 0.8)
159         multiplicity16++;
160     }
161     else
162     {
163       if (TMath::Abs(particle->Eta()) < 2.1)
164         multiplicity21++;
165       if (TMath::Abs(particle->Eta()) < 1.6)
166         multiplicity16++;
167     }
168   }// end of mc particle
169
170   AliRunLoader* runLoader = GetRunLoader();
171   if (!runLoader)
172   {
173     AliDebug(AliLog::kError, "runloader not available");
174     return kFALSE;
175   }
176
177   // TDirectory::TContext restores the current directory is restored when the scope ends.
178   // This helps around ROOT bug #26025 and is good behaviour anyway
179   TDirectory::TContext context(0);
180   AliITSLoader* loader = (AliITSLoader*) runLoader->GetLoader( "ITSLoader" );
181   loader->LoadDigits("READ");
182   TTree* treeD = loader->TreeD();
183   if (!treeD)
184   {
185     AliDebug(AliLog::kError, "Could not retrieve TreeD of ITS");
186     return kFALSE;
187   }
188
189   treeD->SetBranchStatus("*", 0);
190   treeD->SetBranchStatus("ITSDigitsSPD.fTracks*", 1);
191   treeD->SetBranchStatus("ITSDigitsSPD.fCoord1", 1);
192
193   TClonesArray* digits = 0;
194   treeD->SetBranchAddress("ITSDigitsSPD", &digits);
195   if (digits);
196     digits->Clear();
197
198   // each value for both layers
199   Int_t totalNumberOfFO[2];
200   Int_t chipsHitByPrimaries[2];
201   //Int_t chipsHitBySecondaries[2];
202
203   for (Int_t i=0; i<2; ++i)
204   {
205     totalNumberOfFO[i] = 0;
206     chipsHitByPrimaries[i] = 0;
207     //chipsHitBySecondaries[i] = 0;
208   }
209
210   Int_t startSPD = 0; //geom->GetStartSPD();
211   Int_t lastSPD  = 239; //geom->GetLastSPD();
212
213   //printf("%d %d\n", startSPD, lastSPD);
214 // for (Int_t l=0; l<240; ++l) { AliITSgeomTGeo::GetModuleId(l, i, j, k); printf("%d --> %d\n", l, i); }
215
216   // loop over modules (ladders)
217   for (Int_t moduleIndex=startSPD; moduleIndex<lastSPD+1; moduleIndex++)
218   {
219     if (centralRegion)
220     {
221       if ((moduleIndex % 4) == 0 || (moduleIndex % 4) == 3)
222         continue;
223     }
224
225     Int_t currentLayer = 0;
226     if (moduleIndex >= 80)
227       currentLayer = 1;
228
229     treeD->GetEvent(moduleIndex);
230
231     // get number of digits in this module
232     Int_t ndigitsInModule = digits->GetEntriesFast();
233
234     // get number of digits in each chip of the module
235     Int_t ndigitsInChip[5];
236     Bool_t hitByPrimary[5];
237     for( Int_t iChip=0; iChip<5; iChip++)
238     {
239       ndigitsInChip[iChip]=0;
240       hitByPrimary[iChip] = kFALSE;
241     }
242
243     // loop over digits in this module
244     for (Int_t iDig=0; iDig<ndigitsInModule; iDig++)
245     {
246       AliITSdigitSPD* dp = (AliITSdigitSPD*) digits->At(iDig);
247       Int_t column = dp->GetCoord1();
248       Int_t isChip = column / 32;
249
250       //printf("Digit %d has column %d which translates to chip %d\n", iDig, column, isChip);
251
252       fChipsFired->Fill(moduleIndex, isChip);
253
254       ndigitsInChip[isChip]++;
255
256       Bool_t debug = kFALSE;
257
258       // find out which particle caused this chip to fire
259       // if we find at least one primary we consider this chip being fired by a primary
260       for (Int_t track=0; track<10; ++track)
261       {
262         Int_t label = dp->GetTrack(track);
263
264         if (label < 0)
265           continue;
266
267         if (debug)
268           printf("track %d contains label %d\n", track, label);
269
270         TParticle* particle = stack->Particle(label);
271
272         if (!particle)
273         {
274           AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (digit loop).", label));
275           continue;
276         }
277
278         if (debug)
279         {
280           particle->Print();
281           printf("statuscode = %d, p = %f, m = %d\n", particle->GetStatusCode(), particle->P(), particle->GetFirstMother());
282         }
283
284         // TODO delta electrons should be traced back to their mother. this is e.g. solved in AliITSClusterFinderV2::CheckLabels2
285         while (particle->P() < 0.02 && particle->GetStatusCode() == 0 && particle->GetFirstMother() >= 0)
286         {
287           label = particle->GetFirstMother();
288           particle = stack->Particle(label);
289
290           if (!particle)
291             break;
292
293           if (debug)
294           {
295             printf("-->\n");
296             printf("statuscode = %d, p = %f, m = %d\n", particle->GetStatusCode(), particle->P(), particle->GetFirstMother());
297             particle->Print();
298           }
299         }
300
301         if (!particle)
302         {
303           AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (digit loop, finding delta electrons).", label));
304           continue;
305         }
306
307         if (label > nPrim)
308           continue;
309
310         if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
311           continue;
312
313         if (debug)
314           printf("This was a primary (or delta-electron of a primary)!\n");
315
316         hitByPrimary[isChip] = kTRUE;
317       }
318     }
319
320     // get number of FOs in the module
321     for (Int_t ifChip=0; ifChip<5; ifChip++)
322       if( ndigitsInChip[ifChip] >= 1 )
323       {
324         totalNumberOfFO[currentLayer]++;
325         if (hitByPrimary[ifChip])
326         {
327           chipsHitByPrimaries[currentLayer]++;
328         }
329         //else
330         //  chipsHitBySecondaries[currentLayer]++;
331       }
332   }
333
334   //printf("Fired chips: %d %d\n", totalNumberOfFOLayer1, totalNumberOfFOLayer2);
335
336   // now find clusters
337   Int_t clustersLayer[2];
338   clustersLayer[0] = 0;
339   clustersLayer[1] = 0;
340
341   loader->LoadRecPoints("READ");
342   TTree* treeR = loader->TreeR();
343   if (!treeR)
344   {
345     AliDebug(AliLog::kError, "Could not retrieve TreeR of ITS");
346     return kFALSE;
347   }
348
349   // TODO branches!
350   //treeR->SetBranchStatus("*", 0);
351
352   TClonesArray* itsClusters = 0;
353   treeR->SetBranchAddress("ITSRecPoints", &itsClusters);
354
355   Int_t nTreeEntries = treeR->GetEntries();
356   for (Int_t iEntry = 0; iEntry < nTreeEntries; ++iEntry)
357   {
358     if (!treeR->GetEvent(iEntry))
359       continue;
360
361     Int_t nClusters = itsClusters->GetEntriesFast();
362
363     while(nClusters--)
364     {
365       AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters);
366
367       if (cluster->GetLayer() == 0)
368       {
369         clustersLayer[0]++;
370         fClusterZL1->Fill(cluster->GetZ());
371       }
372       else if (cluster->GetLayer() == 1)
373       {
374         clustersLayer[1]++;
375         fClusterZL2->Fill(cluster->GetZ());
376       }
377     }
378   }
379
380   fPrimaryL1->Fill(multiplicity21, totalNumberOfFO[0], chipsHitByPrimaries[0], clustersLayer[0]);
381   fPrimaryL2->Fill(multiplicity16, totalNumberOfFO[1], chipsHitByPrimaries[1], clustersLayer[1]);
382
383   return kTRUE;
384 }
385
386 Bool_t AliHighMultiplicitySelector::Notify()
387 {
388   AliRunLoader* runLoader = GetRunLoader();
389
390   if (runLoader)
391   {
392     AliITSLoader* loader = (AliITSLoader* )runLoader->GetLoader( "ITSLoader" );
393     if (loader)
394     {
395       loader->UnloadDigits();
396       loader->UnloadRecPoints();
397     }
398   }
399
400   return AliSelectorRL::Notify();
401 }
402
403 void AliHighMultiplicitySelector::SlaveTerminate()
404 {
405   // The SlaveTerminate() function is called after all entries or objects
406   // have been processed. When running with PROOF SlaveTerminate() is called
407   // on each slave server.
408
409   AliSelectorRL::SlaveTerminate();
410
411   // Add the histograms to the output on each slave server
412   if (!fOutput)
413   {
414     AliDebug(AliLog::kError, "ERROR: Output list not initialized.");
415     return;
416   }
417
418   fOutput->Add(fChipsFired);
419   fOutput->Add(fPrimaryL1);
420   fOutput->Add(fPrimaryL2);
421   fOutput->Add(fClusterZL1);
422   fOutput->Add(fClusterZL2);
423 }
424
425 void AliHighMultiplicitySelector::Terminate()
426 {
427   // The Terminate() function is the last function to be called during
428   // a query. It always runs on the client, it can be used to present
429   // the results graphically or save the results to file.
430
431   AliSelectorRL::Terminate();
432
433   fChipsFired = dynamic_cast<TH2F*> (fOutput->FindObject("fChipsFired"));
434   fPrimaryL1 = dynamic_cast<TNtuple*> (fOutput->FindObject("fPrimaryL1"));
435   fPrimaryL2 = dynamic_cast<TNtuple*> (fOutput->FindObject("fPrimaryL2"));
436   fClusterZL1 = dynamic_cast<TH1F*> (fOutput->FindObject("fClusterZL1"));
437   fClusterZL2 = dynamic_cast<TH1F*> (fOutput->FindObject("fClusterZL2"));
438
439   if (!fClusterZL1 || !fClusterZL2 || !fChipsFired || !fPrimaryL1 || !fPrimaryL2)
440   {
441     AliError("Histograms not available");
442     return;
443   }
444
445   WriteHistograms();
446 }
447
448 void AliHighMultiplicitySelector::WriteHistograms(const char* filename)
449 {
450   TFile* file = TFile::Open(filename, "RECREATE");
451
452   fChipsFired->Write();
453   fPrimaryL1->Write();
454   fPrimaryL2->Write();
455   fClusterZL1->Write();
456   fClusterZL2->Write();
457
458   file->Close();
459 }
460
461 void AliHighMultiplicitySelector::ReadHistograms(const char* filename)
462 {
463   TFile* file = TFile::Open(filename);
464
465   if (!file)
466     return;
467
468   fPrimaryL1  = dynamic_cast<TNtuple*> (file->Get("fPrimaryL1"));
469   fPrimaryL2  = dynamic_cast<TNtuple*> (file->Get("fPrimaryL2"));
470   fChipsFired  = dynamic_cast<TH2F*> (file->Get("fChipsFired"));
471   fClusterZL1  = dynamic_cast<TH1F*> (file->Get("fClusterZL1"));
472   fClusterZL2  = dynamic_cast<TH1F*> (file->Get("fClusterZL2"));
473
474   #define MULT   1001, -0.5, 1000.5
475   #define BINNING_LAYER1 401, -0.5, 400.5
476   #define BINNING_LAYER2 801, -0.5, 800.5
477
478   fChipsLayer1 = new TH1F("fChipsLayer1", "Layer 1;Fired Chips;Count", BINNING_LAYER1);
479   fChipsLayer2 = new TH1F("fChipsLayer2", "Layer 2;Fired Chips;Count", BINNING_LAYER2);
480
481   fL1vsL2 = new TH2F("fL1vsL2", ";Fired Chips Layer 1;Fired Chips Layer 2", BINNING_LAYER1, BINNING_LAYER2);
482   fMvsL1 = new TH2F("fMvsL1", ";true multiplicity;fired chips layer1", MULT, BINNING_LAYER1);
483   fMvsL2 = new TH2F("fMvsL2", ";true multiplicity;fired chips layer2", MULT, BINNING_LAYER2);
484
485   fClvsL1 = new TH2F("fClvsL1", ";clusters layer1;fired chips layer1", MULT, BINNING_LAYER1);
486   fClvsL2 = new TH2F("fClvsL2", ";clusters layer2;fired chips layer2", MULT, BINNING_LAYER2);
487
488   for (Int_t i = 0; i < fPrimaryL1->GetEntries(); i++)
489   {
490     fPrimaryL1->GetEvent(i);
491     fPrimaryL2->GetEvent(i);
492
493     Int_t multiplicity21 = (Int_t) fPrimaryL1->GetArgs()[0];
494     Int_t multiplicity16 = (Int_t) fPrimaryL2->GetArgs()[0];
495
496     Int_t totalNumberOfFO[2];
497     totalNumberOfFO[0] = (Int_t) fPrimaryL1->GetArgs()[1];
498     totalNumberOfFO[1] = (Int_t) fPrimaryL2->GetArgs()[1];
499
500     Int_t chipsHitByPrimaries[2];
501     chipsHitByPrimaries[0] = (Int_t) fPrimaryL1->GetArgs()[2];
502     chipsHitByPrimaries[1] = (Int_t) fPrimaryL2->GetArgs()[2];
503
504     Int_t clustersLayer[2];
505     clustersLayer[0] = (Int_t) fPrimaryL1->GetArgs()[3];
506     clustersLayer[1] = (Int_t) fPrimaryL2->GetArgs()[3];
507
508     fChipsLayer1->Fill(totalNumberOfFO[0]);
509     fChipsLayer2->Fill(totalNumberOfFO[1]);
510
511     fL1vsL2->Fill(totalNumberOfFO[0], totalNumberOfFO[1]);
512
513     fMvsL1->Fill(multiplicity21, totalNumberOfFO[0]);
514     fMvsL2->Fill(multiplicity16, totalNumberOfFO[1]);
515
516     fClvsL1->Fill(clustersLayer[0], totalNumberOfFO[0]);
517     fClvsL2->Fill(clustersLayer[1], totalNumberOfFO[1]);
518   }
519 }
520
521 TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut)
522 {
523   //
524   // returns the trigger efficiency as function of multiplicity with a given cut
525   //
526
527   //cut and multiply with x-section
528   TH1* allEvents = multVsLayer->ProjectionX("fMvsL_x_total", 1, 1001);
529   //allEvents->Sumw2();
530
531   //cut and multiply with x-section
532   TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, 1001);
533   //proj->Sumw2();
534
535   //new TCanvas; allEvents->DrawCopy(); gPad->SetLogy();
536   //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
537
538   // make probability distribution out of it
539   // TODO binomial errors do not work??? weird...
540   proj->Divide(proj, allEvents, 1, 1, "B");
541
542   return proj;
543 }
544
545 TH1* AliHighMultiplicitySelector::GetXSectionCut(TH1* xSection, TH2* multVsLayer, Int_t cut)
546 {
547   // returns the rel. cross section of the true spectrum that is measured when a cut at <cut> is performed
548
549   TH1* proj = GetTriggerEfficiency(multVsLayer, cut);
550
551   //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
552
553   for (Int_t i=1; i<=proj->GetNbinsX(); i++)
554   {
555     if (i <= xSection->GetNbinsX())
556     {
557       Double_t value = proj->GetBinContent(i) * xSection->GetBinContent(i);
558       Double_t error = 0;
559
560       if (value != 0)
561         error = value * (proj->GetBinError(i) / proj->GetBinContent(i) + xSection->GetBinError(i) / xSection->GetBinContent(i));
562
563       proj->SetBinContent(i, value);
564       proj->SetBinError(i, error);
565     }
566     else
567     {
568       proj->SetBinContent(i, 0);
569       proj->SetBinError(i, 0);
570     }
571   }
572
573   //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
574
575   return proj;
576 }
577
578 void AliHighMultiplicitySelector::MakeGraphs2(const char* title, TH1* xSection, TH2* fMvsL)
579 {
580   TGraph* effGraph = new TGraph;
581   effGraph->SetTitle(Form("%s;Cut on fired chips;mult. where eff. >50%%", title));
582   TGraph* ratioGraph = new TGraph;
583   ratioGraph->SetTitle(Form("%s;Cut on fired chips;x-section_(>=eff. limit) / x-section_(total)", title));
584   TGraph* totalGraph = new TGraph;
585   totalGraph->SetTitle(Form("%s;Cut on fired chips;rel x-section_(>=eff. limit)", title));
586
587   for (Int_t cut = 0; cut <= 300; cut+=50)
588   {
589     TH1* proj = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("clone");
590
591     //proj->Rebin(3);
592     //proj->Scale(1.0 / 3);
593
594     new TCanvas; proj->DrawCopy();
595
596     Int_t limitBin = proj->GetNbinsX()+1;
597     while (limitBin > 1 && proj->GetBinContent(limitBin-1) > 0.5)
598       limitBin--;
599
600     Float_t limit = proj->GetXaxis()->GetBinCenter(limitBin);
601
602     effGraph->SetPoint(effGraph->GetN(), cut, limit);
603
604     proj = GetXSectionCut(xSection, fMvsL, cut);
605
606     Double_t ratio = 0;
607     Double_t total = 0;
608     if (proj->Integral(1, 1001) > 0)
609     {
610       ratio = proj->Integral(proj->FindBin(limit), 1001) / proj->Integral(1, 1001);
611       total = proj->Integral(proj->FindBin(limit), 1001);
612     }
613
614     ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio);
615     totalGraph->SetPoint(totalGraph->GetN(), cut, total);
616
617     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);
618   }
619
620   TCanvas* canvas = new TCanvas(Form("%s_Efficiency", title), Form("%s_Efficiency", title), 1200, 800);
621   canvas->Divide(2, 2);
622
623   canvas->cd(1);
624   effGraph->Draw("A*");
625
626   for (Int_t i=8; i<=10; ++i)
627   {
628     TLine* line = new TLine(0, xSection->GetMean() * i, 300, xSection->GetMean() * i);
629     line->Draw();
630   }
631
632   canvas->cd(2);  ratioGraph->Draw("A*");
633   canvas->cd(3);  gPad->SetLogy(); totalGraph->Draw("A*");
634
635   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
636 }
637
638 void AliHighMultiplicitySelector::MakeGraphs(const char* title, TH1* xSection, TH2* fMvsL, Int_t limit)
639 {
640   // relative x-section, once we have a collision
641   xSection->Scale(1.0 / xSection->Integral());
642
643   TGraph* ratioGraph = new TGraph;
644   ratioGraph->SetTitle(Form("%s;Cut on fired chips;x-section_(>=%d) / x-section_(total)", title, limit));
645   TGraph* totalGraph = new TGraph;
646   totalGraph->SetTitle(Form("%s;Cut on fired chips;rel x-section_(>=%d)", title, limit));
647
648   Double_t max = 0;
649   Int_t bestCut = -1;
650   Double_t bestRatio = -1;
651   Double_t bestTotal = -1;
652   Int_t fullCut = -1;
653   Double_t fullRatio = -1;
654   Double_t fullTotal = -1;
655
656   fMvsL->Sumw2();
657
658   for (Int_t cut = 50; cut <= 300; cut+=2)
659   {
660     TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
661
662     Double_t ratio = 0;
663     Double_t total = 0;
664     if (proj->Integral(1, 1000) > 0)
665     {
666       ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000);
667       total = proj->Integral(limit, 1000);
668     }
669
670     max = TMath::Max(max, total);
671
672     //printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio);
673
674     if (total < max * 0.9 && bestCut == -1)
675     {
676       bestCut = cut;
677       bestRatio = ratio;
678       bestTotal = total;
679     }
680
681     if (ratio == 1 && fullCut == -1)
682     {
683       fullCut = cut;
684       fullRatio = ratio;
685       fullTotal = total;
686     }
687
688     ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio);
689     totalGraph->SetPoint(totalGraph->GetN(), cut, total);
690   }
691
692   if (bestCut != -1)
693     printf("Best cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", bestCut, limit, bestTotal, limit, bestRatio);
694   if (fullCut != -1)
695     printf("100%% cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", fullCut, limit, fullTotal, limit, fullRatio);
696
697   TCanvas* canvas = new TCanvas(Form("%s_RatioXSection_%d", title, limit), Form("%s_RatioXSection_%d", title, limit), 600, 400);
698   ratioGraph->Draw("A*");
699   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
700
701   canvas = new TCanvas(Form("%s_TotalXSection_%d", title, limit), Form("%s_TotalXSection_%d", title, limit), 600, 400);
702   totalGraph->Draw("A*");
703   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
704 }
705
706 void AliHighMultiplicitySelector::JPRPlots()
707 {
708   /*
709
710   gSystem->Load("libPWG0base");
711   .L AliHighMultiplicitySelector.cxx+
712   x = new AliHighMultiplicitySelector();
713   x->ReadHistograms("highmult_hijing100k.root");
714   x->JPRPlots();
715
716   */
717
718   // get x-sections
719   TFile* file = TFile::Open("crosssectionEx.root");
720   if (!file)
721     return;
722
723   TH1* xSections[2];
724   xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
725   xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
726
727   for (Int_t i=0; i<2; ++i)
728   {
729     if (!xSections[i])
730       continue;
731
732     TH1* xSection = xSections[i];
733     TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
734     //Int_t cut = (i == 0) ? 164 : 150; // 8 times the mean
735     //Int_t cut = (i == 0) ? 178 : 166; // 9 times the mean
736     Int_t cut = (i == 0) ? 190 : 182; // 10 times the mean
737
738     // limit is N times the mean
739     Int_t limit = (Int_t) (xSection->GetMean() * 10);
740
741     // 10^28 lum --> 1.2 kHz
742     // 10^31 lum --> 1200 kHz
743     Float_t rate = 1200e3;
744
745     // time in s
746     Float_t lengthRun = 1e6;
747
748     xSection->SetStats(kFALSE);
749     xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
750     xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5));
751     xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 400 : 350);
752     //xSection->GetYaxis()->SetTitle("relative cross section");
753     xSection->GetYaxis()->SetTitleOffset(1.2);
754
755     // relative x-section, once we have a collision
756     xSection->Scale(1.0 / xSection->Integral());
757
758     TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
759
760     Double_t ratio = 0;
761     Double_t total = 0;
762     if (proj->Integral(1, 1000) > 0)
763     {
764       ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000);
765       total = proj->Integral(limit, 1000);
766     }
767
768     printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio);
769
770     TCanvas* canvas = new TCanvas(Form("HMPlots_%d", i), Form("HMPlots_%d", i), 800, 600);
771     canvas->SetLogy();
772     xSection->DrawCopy();
773     proj->SetLineColor(2);
774     proj->SetStats(kFALSE);
775     proj->DrawCopy("SAME");
776
777     TLegend* legend = new TLegend(0.15, 0.15, 0.45, 0.3);
778     legend->SetFillColor(0);
779     legend->AddEntry(xSection, "no trigger");
780     legend->AddEntry(proj, Form("FO trigger > %d chips", cut));
781     legend->Draw();
782
783     TLine* line = new TLine(limit, xSection->GetMinimum() * 0.5, limit, xSection->GetMaximum() * 2);
784     line->SetLineWidth(2);
785     line->Draw();
786
787     canvas->SaveAs(Form("%s.gif", canvas->GetName()));
788
789     TCanvas* canvas2 = new TCanvas(Form("HMPlots_%d_Random", i), Form("HMPlots_%d_Random", i), 800, 600);
790     //canvas2->SetTopMargin(0.05);
791     //canvas2->SetRightMargin(0.05);
792     canvas2->SetLogy();
793     xSection->DrawCopy("HIST");
794
795     TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.3);
796     legend2->SetFillColor(0);
797     legend2->AddEntry(xSection, "no trigger");
798
799     TH1* proj2 = (TH1*) proj->Clone("random");
800     proj2->Reset();
801     // MB lengthRun s 100 Hz
802     Int_t nTrigger = (Int_t) (100 * lengthRun * proj->Integral(1, 1000));
803     proj2->FillRandom(proj, nTrigger);
804
805     TH1* proj3 = (TH1*) proj2->Clone("random_clone");
806     proj3->Divide(proj);
807     proj3->Fit("pol0", "0", "");
808     proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0));
809
810     /*
811     proj2->DrawCopy("SAME");
812     legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio)));
813     */
814
815     proj2 = (TH1*) proj->Clone("random2");
816     proj2->Reset();
817     // 10^31 lum --> 1200 kHz; lengthRun s
818     nTrigger = (Int_t) (rate * proj->Integral(1, 1000) * lengthRun);
819     proj2->FillRandom(proj, nTrigger);
820
821     proj3 = (TH1*) proj2->Clone("random_clone2");
822     proj3->Divide(proj);
823     proj3->Fit("pol0", "0", "");
824     proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0));
825
826     proj2->SetLineColor(4);
827     proj2->SetMarkerStyle(7);
828     proj2->SetMarkerColor(4);
829     proj2->DrawCopy("SAME P");
830     //legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio)));
831     legend2->AddEntry(proj2, Form("FO trigger > %d chips", cut));
832
833     legend2->Draw();
834     line->Draw();
835
836     canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
837     canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
838   }
839 }
840
841 void AliHighMultiplicitySelector::Ntrigger()
842 {
843   //
844   // produces a spectrum created with N triggers
845   // number of triggers and thresholds for the moment fixed
846   //
847
848   /*
849
850   gSystem->Load("libPWG0base");
851   .L AliHighMultiplicitySelector.cxx+g
852   x = new AliHighMultiplicitySelector();
853   x->ReadHistograms("highmult_hijing100k.root");
854   x->Ntrigger();
855
856   */
857
858   // get x-sections
859   TFile* file = TFile::Open("crosssectionEx.root");
860   if (!file)
861     return;
862
863   TH1* xSections[2];
864   xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
865   xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
866
867   // 10^28 lum --> 1.2 kHz
868   // 10^31 lum --> 1200 kHz
869   //Float_t rate = 1200e3;
870   Float_t rate = 1200e3;
871
872   // time in s
873   Float_t lengthRun = 1e6;
874
875   Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
876   Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
877
878   // put to 2 for second layer
879   for (Int_t i=0; i<1; ++i)
880   {
881     if (!xSections[i])
882       continue;
883
884     TH1* xSection = xSections[i];
885     TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
886
887     Int_t nCuts = 6;
888     Int_t cuts[] = { 0, 164, 178, 190, 204, 216 };
889     //Int_t nCuts = 4;
890     //Int_t cuts[] = { 0, 164, 190, 216 };
891     // desired trigger rate in Hz
892     Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 };
893
894     xSection->SetStats(kFALSE);
895     xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
896     xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5));
897     xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 450 : 350);
898     //xSection->GetYaxis()->SetTitle("relative cross section");
899     xSection->GetYaxis()->SetTitleOffset(1.2);
900
901     // relative x-section, once we have a collision
902     xSection->Scale(1.0 / xSection->Integral());
903
904     TCanvas* canvas2 = new TCanvas(Form("HMPlots2_%d_Random", i), Form("HMPlots2_%d_Random", i), 800, 600);
905     canvas2->SetTopMargin(0.05);
906     canvas2->SetRightMargin(0.05);
907     canvas2->SetLogy();
908     xSection->DrawCopy("HIST");
909
910     TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.4);
911     legend2->SetFillColor(0);
912     legend2->AddEntry(xSection, "cross-section");
913
914     for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut)
915     {
916       Int_t cut = cuts[currentCut];
917
918       TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
919
920       Double_t total = 0;
921       if (proj->Integral(1, 1000) > 0)
922         total = proj->Integral(1, 1000);
923
924       printf("Cut at %d: rel. x-section = %e\n", cut, total);
925
926       TH1* proj2 = (TH1*) proj->Clone("random2");
927       proj2->Reset();
928
929       // calculate downscale factor
930       Float_t normalRate = rate * proj->Integral(1, 1000);
931       Float_t downScale = normalRate / ratePerTrigger[currentCut];
932       if (downScale < 1)
933         downScale = 1;
934       Long64_t nTrigger = (Long64_t) (normalRate / downScale * lengthRun);
935
936       Printf("Normal rate is %f, downscale: %f, Simulating %lld triggers", normalRate, downScale, nTrigger);
937       proj2->FillRandom(proj, nTrigger);
938
939       Int_t removed = 0;
940       for (Int_t bin=1; bin<proj2->GetNbinsX(); ++bin)
941         if (proj2->GetBinContent(bin) < 5)
942         {
943           removed += (Int_t) proj2->GetBinContent(bin);
944           proj2->SetBinContent(bin, 0);
945         }
946
947       Printf("Removed %d events", removed);
948
949       proj2->Scale(1.0 / nTrigger * proj->Integral(1, 1000));
950
951       proj2->SetLineColor(colors[currentCut]);
952       proj2->SetMarkerStyle(markers[currentCut]);
953       proj2->SetMarkerColor(colors[currentCut]);
954       proj2->DrawCopy("SAME P");
955       legend2->AddEntry(proj2, Form("%lld evts, FO > %d chips", nTrigger, cut));
956     }
957
958     legend2->Draw();
959
960     canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
961     canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
962   }
963 }
964
965 void AliHighMultiplicitySelector::DrawHistograms()
966 {
967   /*
968
969   gSystem->Load("libPWG0base");
970   .L AliHighMultiplicitySelector.cxx+
971   x = new AliHighMultiplicitySelector();
972   x->ReadHistograms("highmult_pythia.root");
973   x->DrawHistograms();
974
975
976   gSystem->Load("libPWG0base");
977   .L AliHighMultiplicitySelector.cxx+
978   x = new AliHighMultiplicitySelector();
979   x->ReadHistograms("highmult_hijing.root");
980   x->DrawHistograms();
981
982   gSystem->Load("libPWG0base");
983   .L AliHighMultiplicitySelector.cxx+
984   x = new AliHighMultiplicitySelector();
985   x->ReadHistograms("highmult_central.root");
986   x->DrawHistograms();
987
988   gSystem->Load("libPWG0base");
989   .L AliHighMultiplicitySelector.cxx+
990   x = new AliHighMultiplicitySelector();
991   x->ReadHistograms("highmult_hijing100k.root");
992   x->DrawHistograms();
993
994   */
995
996   /*TCanvas* canvas = new TCanvas("chips", "chips", 600, 400);
997
998   fChipsLayer2->SetLineColor(2);
999   fChipsLayer2->SetStats(kFALSE);
1000   fChipsLayer1->SetStats(kFALSE);
1001   fChipsLayer2->SetTitle("");
1002   fChipsLayer2->DrawCopy();
1003   fChipsLayer1->DrawCopy("SAME");
1004   canvas->SaveAs("chips.gif");
1005
1006   canvas = new TCanvas("L1vsL2", "L1vsL2", 600, 400);
1007   fL1vsL2->SetStats(kFALSE);
1008   fL1vsL2->DrawCopy("COLZ");
1009   gPad->SetLogz();
1010   canvas->SaveAs("L1vsL2.gif");*/
1011
1012   TCanvas *canvas = new TCanvas("L1", "L1", 800, 600);
1013   canvas->SetTopMargin(0.05);
1014   canvas->SetRightMargin(0.12);
1015   fMvsL1->SetStats(kFALSE);
1016   fMvsL1->DrawCopy("COLZ");
1017   gPad->SetLogz();
1018
1019   canvas->SaveAs("L1NoCurve.gif");
1020   canvas->SaveAs("L1NoCurve.eps");
1021
1022   // draw corresponding theoretical curve
1023   TF1* func = new TF1("func", "[0]*(1-(1-1/[0])**x)", 1, 1000);
1024   func->SetParameter(0, 400-5*2);
1025   func->DrawCopy("SAME");
1026
1027   canvas->SaveAs("L1.gif");
1028
1029   canvas = new TCanvas("L2", "L2", 600, 400);
1030   //fMvsL2->GetYaxis()->SetRangeUser(0, 150);
1031   fMvsL2->SetStats(kFALSE);
1032   fMvsL2->DrawCopy("COLZ");
1033   gPad->SetLogz();
1034   func->SetParameter(0, 800-5*4);
1035   func->DrawCopy("SAME");
1036   canvas->SaveAs("L2.gif");
1037
1038   // get x-sections
1039   TFile* file = TFile::Open("crosssectionEx.root");
1040   if (file)
1041   {
1042     TH1* xSection2 = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1043     TH1* xSection15 = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1044
1045     MakeGraphs2("Layer1", xSection2, fMvsL1);
1046     return;
1047
1048     // 5 times the mean
1049     //MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 5)); //75 * 2 * 2);
1050     //MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 5)); //(Int_t) (75 * 1.5 * 2));
1051
1052     MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 8));
1053     MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 8));
1054     MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 9));
1055     MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 9));
1056     MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 10));
1057     MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 10));
1058
1059     file->Close();
1060   }
1061
1062   // make spread hists
1063   TGraph* spread = new TGraph;
1064   spread->SetTitle("Spread L1;true multiplicity;RMS");
1065
1066   for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
1067   {
1068     TH1* proj = fMvsL1->ProjectionY("proj", i, i);
1069     spread->SetPoint(spread->GetN(), i, proj->GetRMS());
1070   }
1071
1072   canvas = new TCanvas("SpreadL1", "SpreadL1", 600, 400);
1073   spread->Draw("A*");
1074   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1075
1076   TF1* log = new TF1("log", "[0]*log([1]*x)", 1, 150);
1077   log->SetParLimits(0, 0, 10);
1078   log->SetParLimits(1, 1e-5, 10);
1079
1080   spread->Fit(log, "", "", 1, 150);
1081   log->DrawCopy("SAME");
1082
1083   TGraph* spread2 = new TGraph;
1084   spread2->SetTitle("Spread L2;true multiplicity;RMS");
1085
1086   for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
1087   {
1088     TH1* proj = fMvsL2->ProjectionY("proj", i, i);
1089     spread2->SetPoint(spread2->GetN(), i, proj->GetRMS());
1090   }
1091
1092   canvas = new TCanvas("SpreadL2", "SpreadL2", 600, 400);
1093   spread2->Draw("A*");
1094   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1095
1096   spread2->Fit(log, "", "", 1, 150);
1097   log->DrawCopy("SAME");
1098
1099   canvas = new TCanvas("Clusters_L1", "Clusters_L1", 600, 400);
1100   fClvsL1->SetStats(kFALSE);
1101   fClvsL1->DrawCopy("COLZ");
1102   gPad->SetLogz();
1103
1104   func->SetParameter(0, 400-5*2);
1105   func->DrawCopy("SAME");
1106
1107   canvas->SaveAs("Clusters_L1.gif");
1108
1109   canvas = new TCanvas("Clusters_L2", "Clusters_L2", 600, 400);
1110   fClvsL2->SetStats(kFALSE);
1111   fClvsL2->DrawCopy("COLZ");
1112   gPad->SetLogz();
1113   func->SetParameter(0, 800-5*4);
1114   func->DrawCopy("SAME");
1115   canvas->SaveAs("Clusters_L2.gif");
1116
1117   canvas = new TCanvas("ChipsFired", "ChipsFired", 600, 400);
1118   //fChipsFired->GetYaxis()->SetRangeUser(0, 150);
1119   fChipsFired->SetStats(kFALSE);
1120   fChipsFired->DrawCopy("COLZ");
1121   canvas->SaveAs("ChipsFired.gif");
1122
1123   /*TH1F* tresholdHistL1 = new TH1F("tresholdHistL1", ";chip treshold;<n>", BINNING_LAYER1);
1124   TH1F* tresholdHistL2 = new TH1F("tresholdHistL2", ";chip treshold;<n>", BINNING_LAYER2);
1125
1126   for (Int_t treshold = 0; treshold < 800; treshold++)
1127   {
1128     if (fPrimaryL1->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
1129     {
1130       TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL1->GetHistogram());
1131       if (mult)
1132         tresholdHistL1->Fill(treshold, mult->GetMean());
1133     }
1134     if (fPrimaryL2->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
1135     {
1136       TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL2->GetHistogram());
1137       if (mult)
1138         tresholdHistL2->Fill(treshold, mult->GetMean());
1139     }
1140   }
1141
1142   canvas = new TCanvas("TresholdL1", "TresholdL1", 600, 400);
1143   tresholdHistL1->Draw();
1144   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1145
1146   canvas = new TCanvas("TresholdL2", "TresholdL2", 600, 400);
1147   tresholdHistL2->Draw();
1148   canvas->SaveAs(Form("%s.gif", canvas->GetName()));*/
1149
1150   fPrimaryL1->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff1", "", "prof goff");
1151   fPrimaryL2->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff2", "", "prof goff");
1152
1153   canvas = new TCanvas("Efficiency", "Efficiency", 600, 400);
1154   fPrimaryL1->GetHistogram()->SetStats(kFALSE);
1155   fPrimaryL1->GetHistogram()->Draw();
1156   fPrimaryL2->GetHistogram()->SetLineColor(2);
1157   fPrimaryL2->GetHistogram()->SetStats(kFALSE);
1158   fPrimaryL2->GetHistogram()->Draw("SAME");
1159   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1160
1161   canvas = new TCanvas("ClustersZL1", "ClustersZL1", 600, 400);
1162   fClusterZL1->Rebin(2);
1163   fClusterZL1->Draw();
1164   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1165
1166   canvas = new TCanvas("ClustersZL2", "ClustersZL2", 600, 400);
1167   fClusterZL2->Draw();
1168   fClusterZL2->Rebin(2);
1169   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1170 }