#include <TGraph.h>
#include <TLegend.h>
#include <TLine.h>
+#include <TMath.h>
+#include <TLegend.h>
+#include <TText.h>
#include <AliLog.h>
#include <AliESD.h>
return kFALSE;
}
- Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, AliPWG0Helper::kMB1);
+ static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
+ Bool_t eventTriggered = triggerAnalysis->IsTriggerBitFired(fESD->GetTriggerMask(), AliTriggerAnalysis::kMB1);
if (!eventTriggered)
{
TClonesArray* digits = 0;
treeD->SetBranchAddress("ITSDigitsSPD", &digits);
- if (digits);
+ if (digits)
digits->Clear();
// each value for both layers
fChipsFired = dynamic_cast<TH2F*> (file->Get("fChipsFired"));
fClusterZL1 = dynamic_cast<TH1F*> (file->Get("fClusterZL1"));
fClusterZL2 = dynamic_cast<TH1F*> (file->Get("fClusterZL2"));
+
+ if (!fPrimaryL1 || !fPrimaryL2 || !fChipsFired || !fClusterZL1 || !fClusterZL2)
+ return;
#define MULT 1001, -0.5, 1000.5
#define BINNING_LAYER1 401, -0.5, 400.5
}
}
-TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut)
+TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut, Int_t upperCut)
{
//
// returns the trigger efficiency as function of multiplicity with a given cut
//allEvents->Sumw2();
//cut and multiply with x-section
- TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, 1001);
+ TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, upperCut);
//proj->Sumw2();
//new TCanvas; allEvents->DrawCopy(); gPad->SetLogy();
}
}
-void AliHighMultiplicitySelector::Ntrigger()
+void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
{
//
// produces a spectrum created with N triggers
/*
+ gSystem->Load("libANALYSIS");
gSystem->Load("libPWG0base");
.L AliHighMultiplicitySelector.cxx+g
x = new AliHighMultiplicitySelector();
x->ReadHistograms("highmult_hijing100k.root");
x->Ntrigger();
+ gSystem->Load("libPWG0base");
+ .L AliHighMultiplicitySelector.cxx+g
+ x = new AliHighMultiplicitySelector();
+ x->ReadHistograms("highmult_hijing100k.root");
+ x->Ntrigger(kFALSE);
+
*/
// get x-sections
xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
- // 10^28 lum --> 1.2 kHz
- // 10^31 lum --> 1200 kHz
- //Float_t rate = 1200e3;
- Float_t rate = 250e3;
+ // 10^28 lum --> 1.4 kHz
+ // 10^31 lum --> 1400 kHz
+ //Float_t rate = 1400e3;
+ Float_t rate = 1.4e3;
// time in s
- Float_t lengthRun = 1e6;
+ Float_t lengthRun = 1e5;
Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
//Int_t cuts[] = { 0, 164, 178, 190, 204, 216 };
//Int_t nCuts = 4;
//Int_t cuts[] = { 0, 164, 190, 216 };
- Int_t nCuts = 4;
- Int_t cuts[] = { 0, 114, 145, 165 };
+
+ //Int_t nCuts = 4;
+ //Int_t cuts[] = { 0, 114, 145, 165 };
+ //Float_t ratePerTrigger[] = { 60, 13.3, 13.3, 13.3 };
+
+ Int_t nCuts = 3;
+ Int_t cuts[] = { 0, 114, 148 };
+
+ //Int_t nCuts = 3;
+ //Int_t cuts[] = { 0, 126, 162 };
+ //Float_t ratePerTrigger[] = { 60, 20.0, 20.0 };
// desired trigger rate in Hz
- //Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 };
- Float_t ratePerTrigger[] = { 60, 13.3, 13.3, 13.3 };
+ Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 };
xSection->SetStats(kFALSE);
xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
canvas2->SetTopMargin(0.05);
canvas2->SetRightMargin(0.05);
canvas2->SetLogy();
- xSection->DrawCopy("HIST");
+
+ if (relative)
+ xSection->DrawCopy("HIST");
TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.4);
legend2->SetFillColor(0);
{
Int_t cut = cuts[currentCut];
+ TH1* triggerEff = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
+
TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
+ Float_t triggerLimit = 0;
+ for (Int_t bin = 1; bin <= triggerEff->GetNbinsX(); bin++)
+ if (triggerEff->GetBinContent(bin) < 0.5)
+ triggerLimit = triggerEff->GetXaxis()->GetBinCenter(bin);
+
+ Printf("Efficiency limit (50%%) is at multiplicity %f", triggerLimit);
+
Double_t total = 0;
if (proj->Integral(1, 1000) > 0)
total = proj->Integral(1, 1000);
if (downScale < 1)
downScale = 1;
Long64_t nTrigger = (Long64_t) (normalRate / downScale * lengthRun);
+ nTrigger = TMath::Nint(((Float_t) nTrigger) / 1000) * 1000;
Printf("Normal rate is %f, downscale: %f, Simulating %lld triggers", normalRate, downScale, nTrigger);
+ if (nTrigger == 0)
+ continue;
+
proj2->FillRandom(proj, nTrigger);
Int_t removed = 0;
Printf("Removed %d events", removed);
- proj2->Scale(1.0 / nTrigger * proj->Integral(1, 1000));
+ if (relative)
+ proj2->Scale(1.0 / nTrigger * proj->Integral(1, 1000));
proj2->SetLineColor(colors[currentCut]);
proj2->SetMarkerStyle(markers[currentCut]);
proj2->SetMarkerColor(colors[currentCut]);
- proj2->DrawCopy("SAME P");
- legend2->AddEntry(proj2, Form("%lld evts, FO > %d chips", nTrigger, cut));
+
+ if (relative || currentCut > 0) {
+ proj2->DrawCopy("SAME P");
+ } else
+ proj2->DrawCopy(" P");
+
+ TString eventStr;
+ if (nTrigger > 1e6)
+ {
+ eventStr.Form("%lld M", nTrigger / 1000 / 1000);
+ }
+ else if (nTrigger > 1e3)
+ {
+ eventStr.Form("%lld K", nTrigger / 1000);
+ }
+ else
+ eventStr.Form("%lld", nTrigger);
+
+ TString triggerStr;
+ if (cut == 0)
+ {
+ triggerStr = "minimum bias";
+ }
+ else
+ triggerStr.Form("FO > %d chips", cut);
+
+ legend2->AddEntry(proj2, Form("%s evts, %s", eventStr.Data(), triggerStr.Data()));
+
+ if (triggerLimit > 1)
+ {
+ TLine* line = new TLine(triggerLimit, proj2->GetMinimum(), triggerLimit, proj2->GetMaximum());
+ line->SetLineColor(colors[currentCut]);
+ line->Draw();
+ }
}
legend2->Draw();
}
}
+void AliHighMultiplicitySelector::Contamination()
+{
+ //
+ //
+
+ /*
+
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libPWG0base");
+ .L AliHighMultiplicitySelector.cxx+g
+ x = new AliHighMultiplicitySelector();
+ x->ReadHistograms("highmult_hijing100k.root");
+ x->Contamination();
+
+ */
+
+ // get x-sections
+ TFile* file = TFile::Open("crosssectionEx.root");
+ if (!file)
+ return;
+
+ TH1* xSections[2];
+ xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
+ xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
+
+ // rate = L * sigma
+ // sigma ~ 80 mb (Pythia 14 TeV)
+ // 10^28 lum --> 8e2 Hz
+ // 10^31 lum --> 8e5 Hz
+ Double_t rates[] = { 8e2, 8e3, 8e4, 8e5 };
+
+ Int_t nCuts = 4;
+ Int_t cuts[] = { 104, 134, 154, 170 };
+
+ // put to 2 for second layer
+ for (Int_t i=0; i<1; ++i)
+ {
+ if (!xSections[i])
+ continue;
+
+ // relative x-section, once we have a collision
+ xSections[i]->Scale(1.0 / xSections[i]->Integral());
+
+ Int_t max = xSections[i]->GetNbinsX();
+ max = 500;
+
+ Float_t* xSection = new Float_t[max];
+ for (Int_t mult = 0; mult < max; mult++)
+ xSection[mult] = xSections[i]->GetBinContent(mult+1);
+
+ TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
+
+ TGraph* graph = new TGraph;
+
+ for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut)
+ {
+ Int_t cut = cuts[currentCut];
+ Double_t rate = rates[currentCut];
+ //Double_t rate = rates[3];
+
+ // coll. in 100 ns window
+ Double_t windowSize = 100e-9;
+ //Double_t windowSize = 25e-9;
+ Double_t collPerWindow = windowSize * rate;
+ Printf("coll/window = %f", collPerWindow);
+ Double_t windowsPerSecond = 1.0 / windowSize;
+
+ TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
+ Float_t* triggerEff = new Float_t[max*2];
+ for (Int_t mult = 0; mult < max; mult++)
+ triggerEff[mult] = triggerEffHist->GetBinContent(mult+1);
+
+ Double_t triggerRate = 0;
+ for (Int_t mult = 0; mult < max; mult++)
+ triggerRate += xSection[mult] * triggerEff[mult];
+
+ triggerRate *= TMath::Poisson(1, collPerWindow) * windowsPerSecond;
+
+ Printf("Rate for 1 collision is %f Hz", triggerRate);
+
+ Double_t triggerRate2 = 0;
+ for (Int_t mult = 0; mult < max; mult++)
+ for (Int_t mult2 = mult; mult2 < max; mult2++)
+ if (mult+mult2 < max)
+ triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2];
+
+ triggerRate2 *= TMath::Poisson(2, collPerWindow) * windowsPerSecond;
+
+ Printf("Rate for 2 collisions is %f Hz --> %.1f%%", triggerRate2, triggerRate2 / triggerRate * 100);
+
+ Double_t triggerRate3 = 0;
+
+ for (Int_t mult = 0; mult < max; mult++)
+ for (Int_t mult2 = mult; mult2 < max-mult; mult2++)
+ for (Int_t mult3 = 0; mult3 < max-mult-mult2; mult3++)
+ //if (mult+mult2+mult3 < max)
+ triggerRate3 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * xSection[mult3] * triggerEff[mult+mult2+mult3];
+
+ triggerRate3 *= TMath::Poisson(3, collPerWindow) * windowsPerSecond;
+ //triggerRate3 *= collPerWindow * collPerWindow * rate;
+
+ Printf("Rate for 3 collisions is %f Hz --> %.1f%%", triggerRate3, triggerRate3 / triggerRate * 100);
+
+ Float_t totalContamination = (triggerRate2 + triggerRate3) / (triggerRate + triggerRate2 + triggerRate3);
+
+ Printf("Total contamination is %.1f%%", totalContamination * 100);
+
+ graph->SetPoint(graph->GetN(), cut, totalContamination);
+
+ continue;
+
+ Double_t triggerRate4 = 0;
+ for (Int_t mult = 0; mult < max; mult++)
+ for (Int_t mult2 = mult; mult2 < max-mult; mult2++)
+ for (Int_t mult3 = 0; mult3 < max-mult-mult2; mult3++)
+ for (Int_t mult4 = 0; mult4 < max-mult-mult2-mult3; mult4++)
+ //if (mult+mult2+mult3+mult4 < max)
+ triggerRate4 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * xSection[mult3] * xSection[mult4] * triggerEff[mult+mult2+mult3+mult4];
+
+ //triggerRate4 *= collPerWindow * collPerWindow * collPerWindow * rate;
+ triggerRate4 *= TMath::Poisson(4, collPerWindow) * windowsPerSecond;
+
+ Printf("Rate for 4 collisions is %f Hz --> %.1f%%", triggerRate4, triggerRate4 / triggerRate * 100);
+
+ // general code for n collisions follows, however much slower...
+ /*
+ const Int_t maxdepth = 4;
+ for (Int_t depth = 1; depth <= maxdepth; depth++) {
+ Double_t triggerRate = 0;
+
+ Int_t m[maxdepth];
+ for (Int_t d=0; d<maxdepth; d++)
+ m[d] = 0;
+
+ while (m[0] < max) {
+ Double_t value = 1;
+ Int_t sum = 0;
+ for (Int_t d=0; d<depth; d++) {
+ value *= xSection[m[d]];
+ sum += m[d];
+ }
+
+ if (sum < max) {
+ value *= triggerEff[sum];
+ triggerRate += value;
+ }
+
+ Int_t increase = depth-1;
+ ++m[increase];
+ while (m[increase] == max && increase > 0) {
+ m[increase] = 0;
+ --increase;
+ ++m[increase];
+ }
+ }
+
+ triggerRate *= rate * TMath::Power(collPerWindow, depth - 1);
+
+ Printf("Rate for %d collisions is %f Hz", depth, triggerRate);
+ }*/
+ }
+
+ new TCanvas; graph->Draw("AP*");
+ }
+}
+
+void AliHighMultiplicitySelector::Contamination2()
+{
+ //
+ // produces a spectrum created with N triggers
+ // number of triggers and thresholds for the moment fixed
+ //
+
+ /*
+
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libPWG0base");
+ .L AliHighMultiplicitySelector.cxx+g
+ x = new AliHighMultiplicitySelector();
+ x->ReadHistograms("highmult_hijing100k.root");
+ x->Contamination2();
+
+ */
+
+ // get x-sections
+ TFile* file = TFile::Open("crosssectionEx.root");
+ if (!file)
+ return;
+
+ TH1* xSections[2];
+ xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
+ xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
+
+ Int_t nCuts = 4;
+ Int_t cuts[] = { 104, 134, 154, 170 };
+
+ new TCanvas;
+
+ Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
+ Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
+
+ // put to 2 for second layer
+ for (Int_t i=0; i<1; ++i)
+ {
+ if (!xSections[i])
+ continue;
+
+ // relative x-section, once we have a collision
+ xSections[i]->Scale(1.0 / xSections[i]->Integral());
+
+ Int_t max = xSections[i]->GetNbinsX();
+ max = 500;
+
+ Float_t* xSection = new Float_t[max];
+ for (Int_t mult = 0; mult < max; mult++)
+ xSection[mult] = xSections[i]->GetBinContent(mult+1);
+
+ TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
+
+ for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut)
+ {
+ TGraph* graph = new TGraph;
+ graph->SetMarkerColor(colors[currentCut]);
+ graph->SetMarkerStyle(markers[currentCut]);
+
+ Int_t cut = cuts[currentCut];
+
+ TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
+ Float_t* triggerEff = new Float_t[max*2];
+ for (Int_t mult = 0; mult < max; mult++)
+ triggerEff[mult] = triggerEffHist->GetBinContent(mult+1);
+
+ Double_t triggerRate = 0;
+ for (Int_t mult = 0; mult < max; mult++)
+ triggerRate += xSection[mult] * triggerEff[mult];
+
+ Printf("Raw value for 1 collision is %e", triggerRate);
+
+ Double_t triggerRate2 = 0;
+ for (Int_t mult = 0; mult < max; mult++)
+ for (Int_t mult2 = mult; mult2 < max; mult2++)
+ if (mult+mult2 < max)
+ triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2];
+
+ Printf("Raw value for 2 collisions is %e", triggerRate2);
+
+ for (Double_t doubleRate = 0; doubleRate <= 0.3; doubleRate += 0.005)
+ {
+ Float_t totalContamination = (triggerRate2 * doubleRate) / (triggerRate + triggerRate2 * doubleRate);
+
+ //Printf("Total contamination is %.1f%%", totalContamination * 100);
+
+ graph->SetPoint(graph->GetN(), doubleRate, totalContamination);
+ }
+
+ graph->Draw((currentCut == 0) ? "A*" : "* SAME");
+ graph->GetXaxis()->SetRangeUser(0, 1);
+ }
+ }
+}
+
+void AliHighMultiplicitySelector::Contamination3()
+{
+ //
+ // draws the contamination as function of treshold depending on a number a set of input MC and rate parameters
+ //
+
+ /*
+
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libPWG0base");
+ .L AliHighMultiplicitySelector.cxx+g
+ x = new AliHighMultiplicitySelector();
+ x->ReadHistograms("highmult_hijing100k.root");
+ x->Contamination3();
+
+ */
+
+ // output file
+ TFile* output = TFile::Open("contamination3.root", "RECREATE");
+
+ // get x-sections
+ TFile* file = TFile::Open("crosssectionEx_10TeV.root");
+ if (!file)
+ return;
+
+ TCanvas* c = new TCanvas;
+ c->SetGridx();
+ c->SetGridy();
+
+ TLegend* legend = new TLegend(0.7, 0.2, 1, 0.5);
+ legend->SetNColumns(2);
+
+ TH2* dummy = new TH2F("dummy", ";Layer 1 Threshold;Contamination", 100, 95, 255, 100, 0, 1);
+ dummy->SetStats(kFALSE);
+ dummy->Draw();
+
+ for (Int_t mc = 0; mc < 6; mc++)
+ {
+ TH1* xSections[2];
+ TString str;
+ str.Form("xSection2Ex_%d_%d", mc/3, mc%3);
+ Printf("%s", str.Data());
+ file->cd();
+ xSections[0] = dynamic_cast<TH1*> (gFile->Get(str));
+ //xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
+
+ // prob for a collision in a bunch crossing
+ Int_t nRates = 1;
+ //Float_t rates[] = {0.02, 0.05, 0.1, 0.15, 0.2};
+ Float_t rates[] = {0.0636};
+
+ // bunch crossing rate in Hz
+ Float_t bunchCrossingRate = 24. * 11245.5;
+
+ Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
+ Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
+
+ // put to 2 for second layer
+ for (Int_t i=0; i<1; ++i)
+ {
+ if (!xSections[i])
+ continue;
+
+ // relative x-section, once we have a collision
+ xSections[i]->Scale(1.0 / xSections[i]->Integral());
+
+ Int_t max = xSections[i]->GetNbinsX();
+ max = 500;
+
+ Float_t* xSection = new Float_t[max];
+ for (Int_t mult = 0; mult < max; mult++)
+ xSection[mult] = xSections[i]->GetBinContent(mult+1);
+
+ TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
+
+ for (Int_t currentRate = 0; currentRate<nRates; ++currentRate)
+ {
+ TGraph* graph = new TGraph;
+ graph->SetMarkerColor(colors[currentRate]);
+ graph->SetMarkerStyle(markers[currentRate]);
+
+ TGraph* graph2 = new TGraph;
+
+ Float_t rate = rates[currentRate];
+
+ Double_t singleRate = TMath::Poisson(1, rate);
+ Double_t doubleRate = TMath::Poisson(2, rate);
+ Double_t tripleRate = TMath::Poisson(3, rate);
+
+ Printf("single = %f, double = %f, triple = %f", singleRate, doubleRate, tripleRate);
+
+ for (Int_t cut = 100; cut <= 251; cut += 10)
+ {
+ Printf("Cut at %d", cut);
+
+ TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
+ Float_t* triggerEff = new Float_t[max*2];
+ for (Int_t mult = 0; mult < max; mult++)
+ triggerEff[mult] = triggerEffHist->GetBinContent(mult+1);
+
+ Double_t triggerRate = 0;
+ for (Int_t mult = 0; mult < max; mult++)
+ triggerRate += xSection[mult] * triggerEff[mult];
+
+ //Printf(" Raw value for 1 collision is %e; Rate: %.1f Hz", triggerRate, triggerRate * singleRate * bunchCrossingRate);
+
+ Double_t triggerRate2 = 0;
+ for (Int_t mult = 0; mult < max; mult++)
+ for (Int_t mult2 = mult; mult2 < max; mult2++)
+ if (mult+mult2 < max)
+ triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2];
+
+ //Printf(" Raw value for 2 collisions is %e; Rate: %.1f Hz", triggerRate2, triggerRate2 * doubleRate * bunchCrossingRate);
+
+ Double_t triggerRate3 = 0;
+ for (Int_t mult = 0; mult < max; mult++)
+ for (Int_t mult2 = 0; mult2 < max; mult2++)
+ for (Int_t mult3 = 0; mult3 < max; mult3++)
+ if (mult+mult2+mult3 < max)
+ triggerRate3 += xSection[mult] * xSection[mult2] * xSection[mult3] * triggerEff[mult+mult2+mult3];
+
+ //Printf(" Raw value for 3 collisions is %e; Rate: %.1f Hz", triggerRate3, triggerRate3 * tripleRate * bunchCrossingRate);
+
+ Printf(" Rates: %.1f Hz; %.1f Hz; %.1f Hz", triggerRate * singleRate * bunchCrossingRate, triggerRate2 * doubleRate * bunchCrossingRate, triggerRate3 * tripleRate * bunchCrossingRate);
+
+ Float_t totalTrigger = (triggerRate * singleRate + triggerRate2 * doubleRate + triggerRate3 * tripleRate);
+
+ Printf(" Total trigger rate: %.1f Hz", totalTrigger * bunchCrossingRate);
+
+ //if (totalTrigger * bunchCrossingRate > 200)
+ // continue;
+
+ Float_t totalContamination = (triggerRate2 * doubleRate + triggerRate3 * tripleRate) / totalTrigger;
+ //if (totalContamination > 0.99)
+ // break;
+
+ Printf(" Total contamination is %.1f%%", totalContamination * 100);
+
+ graph->SetPoint(graph->GetN(), cut, totalContamination);
+ graph2->SetPoint(graph->GetN(), cut, totalTrigger * bunchCrossingRate);
+ }
+
+ graph->SetMarkerStyle(mc+20);
+ graph->SetMarkerColor(currentRate+1);
+ graph->Draw("P SAME");
+ graph->GetXaxis()->SetTitle("Layer 1 threshold");
+ graph->GetYaxis()->SetTitle("Contamination");
+ graph->GetYaxis()->SetRangeUser(0, 1);
+
+ if (currentRate == 0)
+ {
+ const char* legendLabel[] = { "Pythia Slope 1", "Pythia Slope 2", "Pythia Slope 3", "Phojet Slope 1", "Phojet Slope 2", "Phojet Slope 3" };
+ legend->AddEntry(graph, legendLabel[mc], "P");
+ }
+
+ output->cd();
+ graph->Write(Form("%s_%d_cont", str.Data(), currentRate));
+ graph2->Write(Form("%s_%d_rate", str.Data(), currentRate));
+ }
+ }
+ }
+
+ output->Close();
+
+ legend->Draw();
+}
+
+void AliHighMultiplicitySelector::Contamination_Reach()
+{
+ // plot the multiplicity reach based on the output from Contamination3()
+ // for each rate case and each MC, a certain number of events is required starting counting from the highest multiplicity
+ // note that the reach of different MC cannot be compared with each other
+
+ /*
+
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libPWG0base");
+ .L AliHighMultiplicitySelector.cxx+g
+ x = new AliHighMultiplicitySelector();
+ x->ReadHistograms("highmult_hijing100k.root");
+ x->Contamination_Reach();
+
+ */
+
+ TCanvas* c = new TCanvas("c", "c", 800, 600);
+ c->Divide(2, 3);
+
+ // prob for a collision in a bunch crossing
+ Int_t nRates = 1;
+ //Float_t rates[] = {0.02, 0.05, 0.1, 0.15, 0.2};
+ Float_t rates[] = {0.0636};
+
+ // bunch crossing rate in Hz
+ Float_t bunchCrossingRate = 24. * 11245.5;
+
+ TH2* dummy = new TH2F("dummy", ";Coll/bunch crossing;multiplicity reach", 100, 0, 0.3, 100, 50, 350);
+ //TH2* dummy = new TH2F("dummy", ";Coll/bunch crossing;fractional cross-section", 100, 0, 0.3, 1000, 1e-6, 0.1);
+ dummy->SetStats(kFALSE);
+
+ const char* legendLabel[] = { "Pythia Slope 1", "Pythia Slope 2", "Pythia Slope 3", "Phojet Slope 1", "Phojet Slope 2", "Phojet Slope 3" };
+
+ TFile* mcFile = TFile::Open("crosssectionEx_10TeV.root");
+ TFile* contFile = TFile::Open("contamination3.root");
+
+ // for comparison: how many MB events can one take at the same time
+ Int_t mbEvents = 2000000 * 500;
+
+ for (Int_t mc=0; mc<6; mc++)
+ {
+ mcFile->cd();
+ TH1* mcHist = (TH1*) gFile->Get(Form("xSection2Ex_%d_%d", mc/3, mc%3));
+ mcHist->Scale(1.0 / mcHist->Integral());
+
+ c->cd(mc+1);//->SetLogy();
+ c->SetGridx();
+ c->SetGridy();
+ dummy->Draw();
+
+ Int_t color = 0;
+ for (Int_t requiredEvents = 300; requiredEvents <= 3000000; requiredEvents *= 10)
+ {
+ TGraph* reach = new TGraph;
+
+ color++;
+ if (color == 5)
+ color++;
+
+ Float_t requiredRate = (Float_t) requiredEvents / 1e6;
+ Printf("Required rate is %f", requiredRate);
+
+ // find reach without trigger
+ Int_t mbReach = 1000;
+ while (mcHist->Integral(mcHist->FindBin(mbReach), mcHist->GetNbinsX()) < (Float_t) requiredEvents / mbEvents && mbReach > 1)
+ mbReach--;
+ Printf("MB reach is %d with %f events", mbReach, mcHist->Integral(mcHist->FindBin(mbReach), mcHist->GetNbinsX()) * mbEvents);
+
+ for (Int_t rate=0; rate<nRates; rate++)
+ {
+ contFile->cd();
+ TGraph* cont = (TGraph*) gFile->Get(Form("xSection2Ex_%d_%d_%d_cont", mc/3, mc%3, rate));
+ TGraph* rateh = (TGraph*) gFile->Get(Form("xSection2Ex_%d_%d_%d_rate", mc/3, mc%3, rate));
+
+ Double_t singleRate = TMath::Poisson(1, rates[rate]);
+ Double_t totalCollRate = singleRate * bunchCrossingRate;
+ Printf("collisions/bc: %f; coll. rate: %f", singleRate, totalCollRate);
+
+ // find 200 Hz limit
+ Int_t low = 100;
+ while (rateh->Eval(low) > 200)
+ low++;
+
+ // find contamination limit
+ Int_t high = 100;
+ while (cont->Eval(high) < 0.9 && high < 250)
+ high++;
+
+ Printf("MC %d Rate %d: Acceptable threshold region is %d to %d", mc, rate, low, high);
+ // find reachable multiplicity; include contamination in rate calculation
+ // TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
+
+ // trigger efficiency as function of multiplicity in range mult <= ... <= high
+ //new TCanvas; fMvsL1->Draw("COLZ");
+ TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL1, low, high);
+
+ //new TCanvas; triggerEffHist->DrawCopy();
+ triggerEffHist->Multiply(mcHist);
+
+ Float_t fractionXSection = triggerEffHist->Integral();
+ Printf("The fraction of the cross-section is %f", fractionXSection);
+
+ //new TCanvas; triggerEffHist->DrawCopy();
+ triggerEffHist->Scale(totalCollRate);
+ //new TCanvas; triggerEffHist->DrawCopy(); gPad->SetLogy();
+
+ Float_t achievedRate = 0;
+ Int_t mult = 1000;
+ while (1)
+ {
+ achievedRate = triggerEffHist->Integral(triggerEffHist->FindBin(mult), triggerEffHist->GetNbinsX());
+
+ if (achievedRate >= requiredRate)
+ break;
+
+ if (mult == 1)
+ break;
+
+ mult--;
+ }
+
+ Printf("Achieved rate %f above multiplicity %d", achievedRate, mult);
+
+ if (achievedRate < requiredRate)
+ {
+ Printf("Achieved rate too low");
+ continue;
+ }
+
+ reach->SetPoint(reach->GetN(), rates[rate], mult);
+ //reach->SetPoint(reach->GetN(), rates[rate], fractionXSection);
+
+
+ //return;
+ }
+
+ reach->SetMarkerColor(color);
+ reach->Draw("SAME*");
+
+ TLine* line = new TLine(0, mbReach, 0.3, mbReach);
+ line->SetLineColor(color);
+ //line->Draw();
+ }
+
+ TText* text = new TText;
+ text->DrawText(0.2, 325, legendLabel[mc]);
+ //return;
+ }
+}
+
void AliHighMultiplicitySelector::DrawHistograms()
{
// draws the histograms
x->ReadHistograms("highmult_central.root");
x->DrawHistograms();
+ gSystem->Load("libANALYSIS");
gSystem->Load("libPWG0base");
- .L AliHighMultiplicitySelector.cxx+
+ .L AliHighMultiplicitySelector.cxx++
x = new AliHighMultiplicitySelector();
x->ReadHistograms("highmult_hijing100k.root");
x->DrawHistograms();
canvas->SaveAs("L1NoCurve.gif");
canvas->SaveAs("L1NoCurve.eps");
+ TLine* line = new TLine(fMvsL1->GetXaxis()->GetXmin(), 150, fMvsL1->GetXaxis()->GetXmax(), 150);
+ line->SetLineWidth(2);
+ line->SetLineColor(kRed);
+ line->Draw();
+
+ canvas->SaveAs("L1NoCurveCut.gif");
+ canvas->SaveAs("L1NoCurveCut.eps");
+
+ return;
+
// draw corresponding theoretical curve
TF1* func = new TF1("func", "[0]*(1-(1-1/[0])**x)", 1, 1000);
func->SetParameter(0, 400-5*2);
// N, normalised to 1 for N=0)
/*
+ gSystem->Load("libANALYSIS");
gSystem->Load("libPWG0base");
.L AliHighMultiplicitySelector.cxx+
x = new AliHighMultiplicitySelector();
return 0;
TH1* xSection;
- xSection = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
+ xSection = static_cast<TH1*> (gFile->Get("xSection2Ex"));
TGraph* result = new TGraph;
Double_t integral = proj->Integral();
- Printf("Cut at %d, intregral is %e", threshold, integral);
+ Printf("Cut at %d, integral is %e", threshold, integral);
result->SetPoint(result->GetN(), threshold, integral);
}
return result;
}
+
+void AliHighMultiplicitySelector::MBComparison()
+{
+ //
+ // finds the threshold from which onwards the number of found events above N times the mean
+ // is higher using a high mult. trigger than just triggering with MB
+ //
+
+ /*
+
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libPWG0base");
+ .L AliHighMultiplicitySelector.cxx+g
+ x = new AliHighMultiplicitySelector();
+ x->ReadHistograms("highmult_hijing100k.root");
+ x->MBComparison();
+
+ */
+
+ // get x-sections
+ TFile* file = TFile::Open("crosssectionEx.root");
+ if (!file)
+ return;
+
+ TH1* xSections[2];
+ xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
+ xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
+
+ // rate = L * sigma
+ // sigma ~ 80 mb (Pythia 14 TeV)
+ // 10^28 lum --> 8e2 Hz
+ // 10^31 lum --> 8e5 Hz
+ Int_t nRates = 4;
+ Double_t rates[] = { 8e2, 8e3, 8e4, 8e5 };
+
+ // threshold in number of fired chips for corresponding rate
+ //Int_t cuts[] = { 104, 134, 154, 170 }; // values for 20 Hz
+ Int_t cuts[] = { 82, 124, 147, 164 }; // values for 50 Hz
+
+ // bandwidth, fractions (for MB, high mult.)
+ Float_t bandwidth = 1e3;
+ Float_t fractionMB = 0.5;
+ Float_t fractionHM = 0.05;
+
+ // different limits to define "interesting events"
+ Int_t nLimits = 9;
+ Int_t limits[] = { 0, 1, 2, 4, 6, 7, 8, 9, 10 };
+
+ // put to 2 for second layer
+ for (Int_t i=0; i<1; ++i)
+ {
+ if (!xSections[i])
+ continue;
+
+ TH1* xSection = xSections[i];
+ TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
+
+ // relative x-section, once we have a collision
+ xSection->Scale(1.0 / xSection->Integral());
+
+ xSection->SetStats(kFALSE);
+ xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
+ xSection->GetXaxis()->SetTitle(Form("true multiplicity in |#eta| < %.1f", (i == 0) ? 2.0 : 1.5));
+ xSection->GetXaxis()->SetRangeUser(0, (i == 0) ? 450 : 350);
+ xSection->GetYaxis()->SetTitleOffset(1.2);
+
+ TCanvas* canvas = new TCanvas("MBComparison", "MBComparison", 1000, 800);
+ canvas->Divide(3, 3);
+
+ for (Int_t currentLimit = 0; currentLimit<nLimits; currentLimit++)
+ {
+ // limit is N times the mean
+ Int_t limit = (Int_t) (xSection->GetMean() * limits[currentLimit]);
+ if (limit < 1)
+ limit = 1;
+
+ TGraph* graphMB = new TGraph;
+ graphMB->SetTitle(Form("Events with %d times above <n> (i.e. n >= %d)", limits[currentLimit], limit));
+ graphMB->SetMarkerStyle(20);
+
+ TGraph* graphBoth = new TGraph;
+ graphBoth->SetMarkerStyle(21);
+
+ Float_t min = bandwidth;
+ Float_t max = 0;
+
+ for (Int_t current = 0; current<nRates; ++current)
+ {
+ Float_t rate = rates[current];
+ Int_t cut = cuts[current];
+
+ TH1* triggerEff = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
+ TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
+
+ Float_t downScaleMB1 = rate / bandwidth;
+ if (downScaleMB1 < 1)
+ downScaleMB1 = 1;
+
+ Float_t downScaleMB2 = rate / (bandwidth * fractionMB);
+ if (downScaleMB2 < 1)
+ downScaleMB2 = 1;
+
+ Float_t downScaleHM = rate * proj->Integral(1, xSection->GetNbinsX()) / (bandwidth * fractionHM);
+ if (downScaleHM < 1)
+ downScaleHM = 1;
+
+ Float_t rateMB1 = rate / downScaleMB1 * xSection->Integral(limit, xSection->GetNbinsX());
+ Float_t rateMB2 = rate / downScaleMB2 * xSection->Integral(limit, xSection->GetNbinsX());
+ Float_t rateHM = rate / downScaleHM * proj->Integral(limit, xSection->GetNbinsX());
+ Float_t combinedRate = rateMB2 + rateHM;
+
+ graphMB->SetPoint(graphMB->GetN(), rate, rateMB1);
+ graphBoth->SetPoint(graphBoth->GetN(), rate, combinedRate);
+
+ min = TMath::Min(min, TMath::Min(rateMB1, combinedRate));
+ max = TMath::Max(min, TMath::Max(rateMB1, combinedRate));
+
+ 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);
+ Printf(" %.2e Hz in MB-only mode", rateMB1);
+ Printf(" %.2e Hz = %.2e Hz + %.2e Hz in MB + high mult. mode", combinedRate, rateMB2, rateHM);
+
+ Printf(" The downscale factors are: %.2f %.2f %.2f", downScaleMB1, downScaleMB2, downScaleHM);
+
+ Int_t triggerLimit = 0;
+ for (Int_t bin = 1; bin <= triggerEff->GetNbinsX(); bin++)
+ if (triggerEff->GetBinContent(bin) < 0.5)
+ triggerLimit = (Int_t) triggerEff->GetXaxis()->GetBinCenter(bin);
+
+ Printf(" Efficiency limit (50%%) is at multiplicity %d", triggerLimit);
+ Float_t fractionGood = proj->Integral(triggerLimit, proj->GetNbinsX()) / proj->Integral();
+ Printf(" %.2f %% of the events are above the trigger limit", fractionGood * 100);
+
+ if (triggerLimit > limit)
+ Printf(" WARNING: interesting events also counted inside the trigger limit");
+
+ Printf(" ");
+ }
+
+ canvas->cd(currentLimit+1)->SetLogx();
+ canvas->cd(currentLimit+1)->SetLogy();
+
+ graphMB->Draw("AP");
+ graphBoth->Draw("P SAME");
+
+ graphMB->GetYaxis()->SetRangeUser(0.5 * min, 2 * max);
+ graphMB->GetXaxis()->SetTitle("Raw rate in Hz");
+ graphMB->GetYaxis()->SetTitle("Event rate in Hz");
+ }
+ }
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff