3 #include "AliHighMultiplicitySelector.h"
11 #include <TParticle.h>
14 #include <TGeoManager.h>
22 #include <AliRunLoader.h>
25 #include <AliITSgeom.h>
26 #include <AliITSLoader.h>
27 #include <AliITSdigitSPD.h>
28 #include <AliITSRecPoint.h>
30 #include "AliPWG0Helper.h"
33 // Selector that produces plots needed for the high multiplicity analysis with the
38 ClassImp(AliHighMultiplicitySelector)
40 AliHighMultiplicitySelector::AliHighMultiplicitySelector() :
54 fCentralRegion(kFALSE)
57 // Constructor. Initialization of pointers
61 AliHighMultiplicitySelector::~AliHighMultiplicitySelector()
67 // histograms are in the output list and deleted when the output
68 // list is deleted by the TSelector dtor
71 void AliHighMultiplicitySelector::SlaveBegin(TTree *tree)
75 AliSelectorRL::SlaveBegin(tree);
77 fChipsFired = new TH2F("fChipsFired", ";Module;Chip;Count", 240, -0.5, 239.5, 5, -0.5, 4.5);
79 fPrimaryL1 = new TNtuple("fPrimaryL1", "", "multiplicity:firedChips:chipsByPrimaries:clusters");
80 fPrimaryL2 = new TNtuple("fPrimaryL2", "", "multiplicity:firedChips:chipsByPrimaries:clusters");
82 fClusterZL1 = new TH1F("fClusterZL1", ";z", 400, -20, 20);
83 fClusterZL2 = new TH1F("fClusterZL2", ";z", 400, -20, 20);
86 void AliHighMultiplicitySelector::Init(TTree* tree)
88 // read the user objects
90 AliSelectorRL::Init(tree);
92 // enable only the needed branches
95 tree->SetBranchStatus("*", 0);
96 tree->SetBranchStatus("fTriggerMask", 1);
98 /*if (fTree->GetCurrentFile())
100 TString fileName(fTree->GetCurrentFile()->GetName());
101 fileName.ReplaceAll("AliESDs", "geometry");
104 TGeoManager::Import(fileName);
109 Bool_t AliHighMultiplicitySelector::Process(Long64_t entry)
115 if (AliSelectorRL::Process(entry) == kFALSE)
118 // Check prerequisites
121 AliDebug(AliLog::kError, "ESD branch not available");
125 Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, AliPWG0Helper::kMB1);
129 AliDebug(AliLog::kDebug, "Event not triggered. Skipping.");
133 AliStack* stack = GetStack();
136 AliDebug(AliLog::kError, "Stack not available");
140 Int_t nPrim = stack->GetNprimary();
141 Int_t multiplicity21 = 0;
142 Int_t multiplicity16 = 0;
144 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
146 AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
148 TParticle* particle = stack->Particle(iMc);
152 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
156 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
161 if (TMath::Abs(particle->Eta()) < 1.05)
163 if (TMath::Abs(particle->Eta()) < 0.8)
168 if (TMath::Abs(particle->Eta()) < 2.1)
170 if (TMath::Abs(particle->Eta()) < 1.6)
173 }// end of mc particle
175 AliRunLoader* runLoader = GetRunLoader();
178 AliDebug(AliLog::kError, "runloader not available");
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();
190 AliDebug(AliLog::kError, "Could not retrieve TreeD of ITS");
194 treeD->SetBranchStatus("*", 0);
195 treeD->SetBranchStatus("ITSDigitsSPD.fTracks*", 1);
196 treeD->SetBranchStatus("ITSDigitsSPD.fCoord1", 1);
198 TClonesArray* digits = 0;
199 treeD->SetBranchAddress("ITSDigitsSPD", &digits);
203 // each value for both layers
204 Int_t totalNumberOfFO[2];
205 Int_t chipsHitByPrimaries[2];
206 //Int_t chipsHitBySecondaries[2];
208 for (Int_t i=0; i<2; ++i)
210 totalNumberOfFO[i] = 0;
211 chipsHitByPrimaries[i] = 0;
212 //chipsHitBySecondaries[i] = 0;
215 Int_t startSPD = 0; //geom->GetStartSPD();
216 Int_t lastSPD = 239; //geom->GetLastSPD();
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); }
221 // loop over modules (ladders)
222 for (Int_t moduleIndex=startSPD; moduleIndex<lastSPD+1; moduleIndex++)
226 if ((moduleIndex % 4) == 0 || (moduleIndex % 4) == 3)
230 Int_t currentLayer = 0;
231 if (moduleIndex >= 80)
234 treeD->GetEvent(moduleIndex);
236 // get number of digits in this module
237 Int_t ndigitsInModule = digits->GetEntriesFast();
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++)
244 ndigitsInChip[iChip]=0;
245 hitByPrimary[iChip] = kFALSE;
248 // loop over digits in this module
249 for (Int_t iDig=0; iDig<ndigitsInModule; iDig++)
251 AliITSdigitSPD* dp = (AliITSdigitSPD*) digits->At(iDig);
252 Int_t column = dp->GetCoord1();
253 Int_t isChip = column / 32;
255 //printf("Digit %d has column %d which translates to chip %d\n", iDig, column, isChip);
257 fChipsFired->Fill(moduleIndex, isChip);
259 ndigitsInChip[isChip]++;
261 Bool_t debug = kFALSE;
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)
267 Int_t label = dp->GetTrack(track);
273 printf("track %d contains label %d\n", track, label);
275 TParticle* particle = stack->Particle(label);
279 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (digit loop).", label));
286 printf("statuscode = %d, p = %f, m = %d\n", particle->GetStatusCode(), particle->P(), particle->GetFirstMother());
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)
292 label = particle->GetFirstMother();
293 particle = stack->Particle(label);
301 printf("statuscode = %d, p = %f, m = %d\n", particle->GetStatusCode(), particle->P(), particle->GetFirstMother());
308 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (digit loop, finding delta electrons).", label));
315 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
319 printf("This was a primary (or delta-electron of a primary)!\n");
321 hitByPrimary[isChip] = kTRUE;
325 // get number of FOs in the module
326 for (Int_t ifChip=0; ifChip<5; ifChip++)
327 if( ndigitsInChip[ifChip] >= 1 )
329 totalNumberOfFO[currentLayer]++;
330 if (hitByPrimary[ifChip])
332 chipsHitByPrimaries[currentLayer]++;
335 // chipsHitBySecondaries[currentLayer]++;
339 //printf("Fired chips: %d %d\n", totalNumberOfFOLayer1, totalNumberOfFOLayer2);
342 Int_t clustersLayer[2];
343 clustersLayer[0] = 0;
344 clustersLayer[1] = 0;
346 loader->LoadRecPoints("READ");
347 TTree* treeR = loader->TreeR();
350 AliDebug(AliLog::kError, "Could not retrieve TreeR of ITS");
355 //treeR->SetBranchStatus("*", 0);
357 TClonesArray* itsClusters = 0;
358 treeR->SetBranchAddress("ITSRecPoints", &itsClusters);
360 Int_t nTreeEntries = treeR->GetEntries();
361 for (Int_t iEntry = 0; iEntry < nTreeEntries; ++iEntry)
363 if (!treeR->GetEvent(iEntry))
366 Int_t nClusters = itsClusters->GetEntriesFast();
370 AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters);
372 if (cluster->GetLayer() == 0)
375 fClusterZL1->Fill(cluster->GetZ());
377 else if (cluster->GetLayer() == 1)
380 fClusterZL2->Fill(cluster->GetZ());
385 fPrimaryL1->Fill(multiplicity21, totalNumberOfFO[0], chipsHitByPrimaries[0], clustersLayer[0]);
386 fPrimaryL2->Fill(multiplicity16, totalNumberOfFO[1], chipsHitByPrimaries[1], clustersLayer[1]);
391 Bool_t AliHighMultiplicitySelector::Notify()
393 // get next ITS runloader
395 AliRunLoader* runLoader = GetRunLoader();
399 AliITSLoader* loader = (AliITSLoader* )runLoader->GetLoader( "ITSLoader" );
402 loader->UnloadDigits();
403 loader->UnloadRecPoints();
407 return AliSelectorRL::Notify();
410 void AliHighMultiplicitySelector::SlaveTerminate()
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.
416 AliSelectorRL::SlaveTerminate();
418 // Add the histograms to the output on each slave server
421 AliDebug(AliLog::kError, "ERROR: Output list not initialized.");
425 fOutput->Add(fChipsFired);
426 fOutput->Add(fPrimaryL1);
427 fOutput->Add(fPrimaryL2);
428 fOutput->Add(fClusterZL1);
429 fOutput->Add(fClusterZL2);
432 void AliHighMultiplicitySelector::Terminate()
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.
438 AliSelectorRL::Terminate();
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"));
446 if (!fClusterZL1 || !fClusterZL2 || !fChipsFired || !fPrimaryL1 || !fPrimaryL2)
448 AliError("Histograms not available");
455 void AliHighMultiplicitySelector::WriteHistograms(const char* filename)
457 // write the histograms to a file
459 TFile* file = TFile::Open(filename, "RECREATE");
461 fChipsFired->Write();
464 fClusterZL1->Write();
465 fClusterZL2->Write();
470 void AliHighMultiplicitySelector::ReadHistograms(const char* filename)
472 // read the data from a file and fill histograms
474 TFile* file = TFile::Open(filename);
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"));
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
489 fChipsLayer1 = new TH1F("fChipsLayer1", "Layer 1;Fired Chips;Count", BINNING_LAYER1);
490 fChipsLayer2 = new TH1F("fChipsLayer2", "Layer 2;Fired Chips;Count", BINNING_LAYER2);
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);
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);
499 for (Int_t i = 0; i < fPrimaryL1->GetEntries(); i++)
501 fPrimaryL1->GetEvent(i);
502 fPrimaryL2->GetEvent(i);
504 Int_t multiplicity21 = (Int_t) fPrimaryL1->GetArgs()[0];
505 Int_t multiplicity16 = (Int_t) fPrimaryL2->GetArgs()[0];
507 Int_t totalNumberOfFO[2];
508 totalNumberOfFO[0] = (Int_t) fPrimaryL1->GetArgs()[1];
509 totalNumberOfFO[1] = (Int_t) fPrimaryL2->GetArgs()[1];
511 Int_t chipsHitByPrimaries[2];
512 chipsHitByPrimaries[0] = (Int_t) fPrimaryL1->GetArgs()[2];
513 chipsHitByPrimaries[1] = (Int_t) fPrimaryL2->GetArgs()[2];
515 Int_t clustersLayer[2];
516 clustersLayer[0] = (Int_t) fPrimaryL1->GetArgs()[3];
517 clustersLayer[1] = (Int_t) fPrimaryL2->GetArgs()[3];
519 fChipsLayer1->Fill(totalNumberOfFO[0]);
520 fChipsLayer2->Fill(totalNumberOfFO[1]);
522 fL1vsL2->Fill(totalNumberOfFO[0], totalNumberOfFO[1]);
524 fMvsL1->Fill(multiplicity21, totalNumberOfFO[0]);
525 fMvsL2->Fill(multiplicity16, totalNumberOfFO[1]);
527 fClvsL1->Fill(clustersLayer[0], totalNumberOfFO[0]);
528 fClvsL2->Fill(clustersLayer[1], totalNumberOfFO[1]);
532 TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut)
535 // returns the trigger efficiency as function of multiplicity with a given cut
538 //cut and multiply with x-section
539 TH1* allEvents = multVsLayer->ProjectionX("fMvsL_x_total", 1, 1001);
540 //allEvents->Sumw2();
542 //cut and multiply with x-section
543 TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, 1001);
546 //new TCanvas; allEvents->DrawCopy(); gPad->SetLogy();
547 //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
549 // make probability distribution out of it
550 // TODO binomial errors do not work??? weird...
551 proj->Divide(proj, allEvents, 1, 1, "B");
556 TH1* AliHighMultiplicitySelector::GetXSectionCut(TH1* xSection, TH2* multVsLayer, Int_t cut)
558 // returns the rel. cross section of the true spectrum that is measured when a cut at <cut> is performed
560 TH1* proj = GetTriggerEfficiency(multVsLayer, cut);
562 //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
564 for (Int_t i=1; i<=proj->GetNbinsX(); i++)
566 if (i <= xSection->GetNbinsX())
568 Double_t value = proj->GetBinContent(i) * xSection->GetBinContent(i);
572 error = value * (proj->GetBinError(i) / proj->GetBinContent(i) + xSection->GetBinError(i) / xSection->GetBinContent(i));
574 proj->SetBinContent(i, value);
575 proj->SetBinError(i, error);
579 proj->SetBinContent(i, 0);
580 proj->SetBinError(i, 0);
584 //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
589 void AliHighMultiplicitySelector::MakeGraphs2(const char* title, TH1* xSection, TH2* fMvsL)
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));
600 for (Int_t cut = 0; cut <= 300; cut+=50)
602 TH1* proj = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("clone");
605 //proj->Scale(1.0 / 3);
607 new TCanvas; proj->DrawCopy();
609 Int_t limitBin = proj->GetNbinsX()+1;
610 while (limitBin > 1 && proj->GetBinContent(limitBin-1) > 0.5)
613 Float_t limit = proj->GetXaxis()->GetBinCenter(limitBin);
615 effGraph->SetPoint(effGraph->GetN(), cut, limit);
617 proj = GetXSectionCut(xSection, fMvsL, cut);
621 if (proj->Integral(1, 1001) > 0)
623 ratio = proj->Integral(proj->FindBin(limit), 1001) / proj->Integral(1, 1001);
624 total = proj->Integral(proj->FindBin(limit), 1001);
627 ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio);
628 totalGraph->SetPoint(totalGraph->GetN(), cut, total);
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);
633 TCanvas* canvas = new TCanvas(Form("%s_Efficiency", title), Form("%s_Efficiency", title), 1200, 800);
634 canvas->Divide(2, 2);
637 effGraph->Draw("A*");
639 for (Int_t i=8; i<=10; ++i)
641 TLine* line = new TLine(0, xSection->GetMean() * i, 300, xSection->GetMean() * i);
645 canvas->cd(2); ratioGraph->Draw("A*");
646 canvas->cd(3); gPad->SetLogy(); totalGraph->Draw("A*");
648 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
651 void AliHighMultiplicitySelector::MakeGraphs(const char* title, TH1* xSection, TH2* fMvsL, Int_t limit)
653 // relative x-section, once we have a collision
655 xSection->Scale(1.0 / xSection->Integral());
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));
664 Double_t bestRatio = -1;
665 Double_t bestTotal = -1;
667 Double_t fullRatio = -1;
668 Double_t fullTotal = -1;
672 for (Int_t cut = 50; cut <= 300; cut+=2)
674 TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
678 if (proj->Integral(1, 1000) > 0)
680 ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000);
681 total = proj->Integral(limit, 1000);
684 max = TMath::Max(max, total);
686 //printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio);
688 if (total < max * 0.9 && bestCut == -1)
695 if (ratio == 1 && fullCut == -1)
702 ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio);
703 totalGraph->SetPoint(totalGraph->GetN(), cut, total);
707 printf("Best cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", bestCut, limit, bestTotal, limit, bestRatio);
709 printf("100%% cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", fullCut, limit, fullTotal, limit, fullRatio);
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()));
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()));
720 void AliHighMultiplicitySelector::JPRPlots()
726 gSystem->Load("libPWG0base");
727 .L AliHighMultiplicitySelector.cxx+
728 x = new AliHighMultiplicitySelector();
729 x->ReadHistograms("highmult_hijing100k.root");
735 TFile* file = TFile::Open("crosssectionEx.root");
740 xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
741 xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
743 for (Int_t i=0; i<2; ++i)
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
754 // limit is N times the mean
755 Int_t limit = (Int_t) (xSection->GetMean() * 10);
757 // 10^28 lum --> 1.2 kHz
758 // 10^31 lum --> 1200 kHz
759 Float_t rate = 1200e3;
762 Float_t lengthRun = 1e6;
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);
771 // relative x-section, once we have a collision
772 xSection->Scale(1.0 / xSection->Integral());
774 TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
778 if (proj->Integral(1, 1000) > 0)
780 ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000);
781 total = proj->Integral(limit, 1000);
784 printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio);
786 TCanvas* canvas = new TCanvas(Form("HMPlots_%d", i), Form("HMPlots_%d", i), 800, 600);
788 xSection->DrawCopy();
789 proj->SetLineColor(2);
790 proj->SetStats(kFALSE);
791 proj->DrawCopy("SAME");
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));
799 TLine* line = new TLine(limit, xSection->GetMinimum() * 0.5, limit, xSection->GetMaximum() * 2);
800 line->SetLineWidth(2);
803 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
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);
809 xSection->DrawCopy("HIST");
811 TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.3);
812 legend2->SetFillColor(0);
813 legend2->AddEntry(xSection, "no trigger");
815 TH1* proj2 = (TH1*) proj->Clone("random");
817 // MB lengthRun s 100 Hz
818 Int_t nTrigger = (Int_t) (100 * lengthRun * proj->Integral(1, 1000));
819 proj2->FillRandom(proj, nTrigger);
821 TH1* proj3 = (TH1*) proj2->Clone("random_clone");
823 proj3->Fit("pol0", "0", "");
824 proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0));
827 proj2->DrawCopy("SAME");
828 legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio)));
831 proj2 = (TH1*) proj->Clone("random2");
833 // 10^31 lum --> 1200 kHz; lengthRun s
834 nTrigger = (Int_t) (rate * proj->Integral(1, 1000) * lengthRun);
835 proj2->FillRandom(proj, nTrigger);
837 proj3 = (TH1*) proj2->Clone("random_clone2");
839 proj3->Fit("pol0", "0", "");
840 proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0));
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));
852 canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
853 canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
857 void AliHighMultiplicitySelector::Ntrigger()
860 // produces a spectrum created with N triggers
861 // number of triggers and thresholds for the moment fixed
866 gSystem->Load("libPWG0base");
867 .L AliHighMultiplicitySelector.cxx+g
868 x = new AliHighMultiplicitySelector();
869 x->ReadHistograms("highmult_hijing100k.root");
875 TFile* file = TFile::Open("crosssectionEx.root");
880 xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
881 xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
883 // 10^28 lum --> 1.2 kHz
884 // 10^31 lum --> 1200 kHz
885 //Float_t rate = 1200e3;
886 Float_t rate = 250e3;
889 Float_t lengthRun = 1e6;
891 Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
892 Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
894 // put to 2 for second layer
895 for (Int_t i=0; i<1; ++i)
900 TH1* xSection = xSections[i];
901 TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
904 //Int_t cuts[] = { 0, 164, 178, 190, 204, 216 };
906 //Int_t cuts[] = { 0, 164, 190, 216 };
908 Int_t cuts[] = { 0, 114, 145, 165 };
910 // desired trigger rate in Hz
911 //Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 };
912 Float_t ratePerTrigger[] = { 60, 13.3, 13.3, 13.3 };
914 xSection->SetStats(kFALSE);
915 xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
916 xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5));
917 xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 450 : 350);
918 //xSection->GetYaxis()->SetTitle("relative cross section");
919 xSection->GetYaxis()->SetTitleOffset(1.2);
921 // relative x-section, once we have a collision
922 xSection->Scale(1.0 / xSection->Integral());
924 TCanvas* canvas2 = new TCanvas(Form("HMPlots2_%d_Random", i), Form("HMPlots2_%d_Random", i), 800, 600);
925 canvas2->SetTopMargin(0.05);
926 canvas2->SetRightMargin(0.05);
928 xSection->DrawCopy("HIST");
930 TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.4);
931 legend2->SetFillColor(0);
932 legend2->AddEntry(xSection, "cross-section");
934 for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut)
936 Int_t cut = cuts[currentCut];
938 TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
941 if (proj->Integral(1, 1000) > 0)
942 total = proj->Integral(1, 1000);
944 printf("Cut at %d: rel. x-section = %e\n", cut, total);
946 TH1* proj2 = (TH1*) proj->Clone("random2");
949 // calculate downscale factor
950 Float_t normalRate = rate * proj->Integral(1, 1000);
951 Float_t downScale = normalRate / ratePerTrigger[currentCut];
954 Long64_t nTrigger = (Long64_t) (normalRate / downScale * lengthRun);
956 Printf("Normal rate is %f, downscale: %f, Simulating %lld triggers", normalRate, downScale, nTrigger);
957 proj2->FillRandom(proj, nTrigger);
960 for (Int_t bin=1; bin<proj2->GetNbinsX(); ++bin)
961 if (proj2->GetBinContent(bin) < 5)
963 removed += (Int_t) proj2->GetBinContent(bin);
964 proj2->SetBinContent(bin, 0);
967 Printf("Removed %d events", removed);
969 proj2->Scale(1.0 / nTrigger * proj->Integral(1, 1000));
971 proj2->SetLineColor(colors[currentCut]);
972 proj2->SetMarkerStyle(markers[currentCut]);
973 proj2->SetMarkerColor(colors[currentCut]);
974 proj2->DrawCopy("SAME P");
975 legend2->AddEntry(proj2, Form("%lld evts, FO > %d chips", nTrigger, cut));
980 canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
981 canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
985 void AliHighMultiplicitySelector::DrawHistograms()
987 // draws the histograms
991 gSystem->Load("libPWG0base");
992 .L AliHighMultiplicitySelector.cxx+
993 x = new AliHighMultiplicitySelector();
994 x->ReadHistograms("highmult_pythia.root");
998 gSystem->Load("libPWG0base");
999 .L AliHighMultiplicitySelector.cxx+
1000 x = new AliHighMultiplicitySelector();
1001 x->ReadHistograms("highmult_hijing.root");
1002 x->DrawHistograms();
1004 gSystem->Load("libPWG0base");
1005 .L AliHighMultiplicitySelector.cxx+
1006 x = new AliHighMultiplicitySelector();
1007 x->ReadHistograms("highmult_central.root");
1008 x->DrawHistograms();
1010 gSystem->Load("libPWG0base");
1011 .L AliHighMultiplicitySelector.cxx+
1012 x = new AliHighMultiplicitySelector();
1013 x->ReadHistograms("highmult_hijing100k.root");
1014 x->DrawHistograms();
1018 /*TCanvas* canvas = new TCanvas("chips", "chips", 600, 400);
1020 fChipsLayer2->SetLineColor(2);
1021 fChipsLayer2->SetStats(kFALSE);
1022 fChipsLayer1->SetStats(kFALSE);
1023 fChipsLayer2->SetTitle("");
1024 fChipsLayer2->DrawCopy();
1025 fChipsLayer1->DrawCopy("SAME");
1026 canvas->SaveAs("chips.gif");
1028 canvas = new TCanvas("L1vsL2", "L1vsL2", 600, 400);
1029 fL1vsL2->SetStats(kFALSE);
1030 fL1vsL2->DrawCopy("COLZ");
1032 canvas->SaveAs("L1vsL2.gif");*/
1034 TCanvas *canvas = new TCanvas("L1", "L1", 800, 600);
1035 canvas->SetTopMargin(0.05);
1036 canvas->SetRightMargin(0.12);
1037 fMvsL1->SetStats(kFALSE);
1038 fMvsL1->DrawCopy("COLZ");
1041 canvas->SaveAs("L1NoCurve.gif");
1042 canvas->SaveAs("L1NoCurve.eps");
1044 // draw corresponding theoretical curve
1045 TF1* func = new TF1("func", "[0]*(1-(1-1/[0])**x)", 1, 1000);
1046 func->SetParameter(0, 400-5*2);
1047 func->DrawCopy("SAME");
1049 canvas->SaveAs("L1.gif");
1051 canvas = new TCanvas("L2", "L2", 600, 400);
1052 //fMvsL2->GetYaxis()->SetRangeUser(0, 150);
1053 fMvsL2->SetStats(kFALSE);
1054 fMvsL2->DrawCopy("COLZ");
1056 func->SetParameter(0, 800-5*4);
1057 func->DrawCopy("SAME");
1058 canvas->SaveAs("L2.gif");
1061 TFile* file = TFile::Open("crosssectionEx.root");
1064 TH1* xSection2 = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1065 TH1* xSection15 = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
1067 MakeGraphs2("Layer1", xSection2, fMvsL1);
1071 //MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 5)); //75 * 2 * 2);
1072 //MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 5)); //(Int_t) (75 * 1.5 * 2));
1074 MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 8));
1075 MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 8));
1076 MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 9));
1077 MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 9));
1078 MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 10));
1079 MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 10));
1084 // make spread hists
1085 TGraph* spread = new TGraph;
1086 spread->SetTitle("Spread L1;true multiplicity;RMS");
1088 for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
1090 TH1* proj = fMvsL1->ProjectionY("proj", i, i);
1091 spread->SetPoint(spread->GetN(), i, proj->GetRMS());
1094 canvas = new TCanvas("SpreadL1", "SpreadL1", 600, 400);
1096 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1098 TF1* log = new TF1("log", "[0]*log([1]*x)", 1, 150);
1099 log->SetParLimits(0, 0, 10);
1100 log->SetParLimits(1, 1e-5, 10);
1102 spread->Fit(log, "", "", 1, 150);
1103 log->DrawCopy("SAME");
1105 TGraph* spread2 = new TGraph;
1106 spread2->SetTitle("Spread L2;true multiplicity;RMS");
1108 for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
1110 TH1* proj = fMvsL2->ProjectionY("proj", i, i);
1111 spread2->SetPoint(spread2->GetN(), i, proj->GetRMS());
1114 canvas = new TCanvas("SpreadL2", "SpreadL2", 600, 400);
1115 spread2->Draw("A*");
1116 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1118 spread2->Fit(log, "", "", 1, 150);
1119 log->DrawCopy("SAME");
1121 canvas = new TCanvas("Clusters_L1", "Clusters_L1", 600, 400);
1122 fClvsL1->SetStats(kFALSE);
1123 fClvsL1->DrawCopy("COLZ");
1126 func->SetParameter(0, 400-5*2);
1127 func->DrawCopy("SAME");
1129 canvas->SaveAs("Clusters_L1.gif");
1131 canvas = new TCanvas("Clusters_L2", "Clusters_L2", 600, 400);
1132 fClvsL2->SetStats(kFALSE);
1133 fClvsL2->DrawCopy("COLZ");
1135 func->SetParameter(0, 800-5*4);
1136 func->DrawCopy("SAME");
1137 canvas->SaveAs("Clusters_L2.gif");
1139 canvas = new TCanvas("ChipsFired", "ChipsFired", 600, 400);
1140 //fChipsFired->GetYaxis()->SetRangeUser(0, 150);
1141 fChipsFired->SetStats(kFALSE);
1142 fChipsFired->DrawCopy("COLZ");
1143 canvas->SaveAs("ChipsFired.gif");
1145 /*TH1F* tresholdHistL1 = new TH1F("tresholdHistL1", ";chip treshold;<n>", BINNING_LAYER1);
1146 TH1F* tresholdHistL2 = new TH1F("tresholdHistL2", ";chip treshold;<n>", BINNING_LAYER2);
1148 for (Int_t treshold = 0; treshold < 800; treshold++)
1150 if (fPrimaryL1->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
1152 TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL1->GetHistogram());
1154 tresholdHistL1->Fill(treshold, mult->GetMean());
1156 if (fPrimaryL2->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
1158 TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL2->GetHistogram());
1160 tresholdHistL2->Fill(treshold, mult->GetMean());
1164 canvas = new TCanvas("TresholdL1", "TresholdL1", 600, 400);
1165 tresholdHistL1->Draw();
1166 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1168 canvas = new TCanvas("TresholdL2", "TresholdL2", 600, 400);
1169 tresholdHistL2->Draw();
1170 canvas->SaveAs(Form("%s.gif", canvas->GetName()));*/
1172 fPrimaryL1->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff1", "", "prof goff");
1173 fPrimaryL2->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff2", "", "prof goff");
1175 canvas = new TCanvas("Efficiency", "Efficiency", 600, 400);
1176 fPrimaryL1->GetHistogram()->SetStats(kFALSE);
1177 fPrimaryL1->GetHistogram()->Draw();
1178 fPrimaryL2->GetHistogram()->SetLineColor(2);
1179 fPrimaryL2->GetHistogram()->SetStats(kFALSE);
1180 fPrimaryL2->GetHistogram()->Draw("SAME");
1181 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1183 canvas = new TCanvas("ClustersZL1", "ClustersZL1", 600, 400);
1184 fClusterZL1->Rebin(2);
1185 fClusterZL1->Draw();
1186 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1188 canvas = new TCanvas("ClustersZL2", "ClustersZL2", 600, 400);
1189 fClusterZL2->Draw();
1190 fClusterZL2->Rebin(2);
1191 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1194 TGraph* AliHighMultiplicitySelector::IntFractRate()
1196 // A plot which shows the fractional rate above any threshold
1197 // as function of threshold (i.e. the integral of dSigma/dN as function of
1198 // N, normalised to 1 for N=0)
1201 gSystem->Load("libPWG0base");
1202 .L AliHighMultiplicitySelector.cxx+
1203 x = new AliHighMultiplicitySelector();
1204 x->ReadHistograms("highmult_hijing100k.root");
1209 TFile* file = TFile::Open("crosssectionEx.root");
1214 xSection = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
1216 TGraph* result = new TGraph;
1218 for (Int_t threshold = 0; threshold < 300; threshold += 2)
1220 TH1* proj = GetXSectionCut(xSection, fMvsL1, threshold);
1222 //new TCanvas; proj->DrawCopy();
1224 Double_t integral = proj->Integral();
1226 Printf("Cut at %d, intregral is %e", threshold, integral);
1228 result->SetPoint(result->GetN(), threshold, integral);
1231 TCanvas* canvas = new TCanvas("IntFractRate", "IntFractRate", 600, 400);
1235 result->GetXaxis()->SetTitle("threshold");
1236 result->GetYaxis()->SetTitle("integrated fractional rate above threshold");
1238 canvas->SaveAs(Form("%s.gif", canvas->GetName()));