+++ /dev/null
-/* $Id$ */
-
-#include "AliHighMultiplicitySelector.h"
-
-#include <TVector3.h>
-#include <TFile.h>
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TNtuple.h>
-#include <TProfile.h>
-#include <TParticle.h>
-#include <TCanvas.h>
-#include <TTree.h>
-#include <TGeoManager.h>
-#include <TF1.h>
-#include <TGraph.h>
-#include <TLegend.h>
-#include <TLine.h>
-#include <TMath.h>
-#include <TLegend.h>
-#include <TText.h>
-
-#include <AliLog.h>
-#include <AliESD.h>
-#include <AliRunLoader.h>
-#include <AliStack.h>
-
-#include <AliITSgeom.h>
-#include <AliITSLoader.h>
-#include <AliITSdigitSPD.h>
-#include <AliITSRecPoint.h>
-
-#include "AliPWG0Helper.h"
-
-//
-// Selector that produces plots needed for the high multiplicity analysis with the
-// pixel detector
-// Needs cluster file
-//
-
-ClassImp(AliHighMultiplicitySelector)
-
-AliHighMultiplicitySelector::AliHighMultiplicitySelector() :
- AliSelectorRL(),
- fChipsLayer1(0),
- fChipsLayer2(0),
- fL1vsL2(0),
- fMvsL1(0),
- fMvsL2(0),
- fChipsFired(0),
- fPrimaryL1(0),
- fPrimaryL2(0),
- fClusterZL1(0),
- fClusterZL2(0),
- fClvsL1(0),
- fClvsL2(0),
- fCentralRegion(kFALSE)
-{
- //
- // Constructor. Initialization of pointers
- //
-}
-
-AliHighMultiplicitySelector::~AliHighMultiplicitySelector()
-{
- //
- // Destructor
- //
-
- // histograms are in the output list and deleted when the output
- // list is deleted by the TSelector dtor
-}
-
-void AliHighMultiplicitySelector::SlaveBegin(TTree *tree)
-{
- // create histograms
-
- AliSelectorRL::SlaveBegin(tree);
-
- fChipsFired = new TH2F("fChipsFired", ";Module;Chip;Count", 240, -0.5, 239.5, 5, -0.5, 4.5);
-
- fPrimaryL1 = new TNtuple("fPrimaryL1", "", "multiplicity:firedChips:chipsByPrimaries:clusters");
- fPrimaryL2 = new TNtuple("fPrimaryL2", "", "multiplicity:firedChips:chipsByPrimaries:clusters");
-
- fClusterZL1 = new TH1F("fClusterZL1", ";z", 400, -20, 20);
- fClusterZL2 = new TH1F("fClusterZL2", ";z", 400, -20, 20);
-}
-
-void AliHighMultiplicitySelector::Init(TTree* tree)
-{
- // read the user objects
-
- AliSelectorRL::Init(tree);
-
- // enable only the needed branches
- if (tree)
- {
- tree->SetBranchStatus("*", 0);
- tree->SetBranchStatus("fTriggerMask", 1);
-
- /*if (fTree->GetCurrentFile())
- {
- TString fileName(fTree->GetCurrentFile()->GetName());
- fileName.ReplaceAll("AliESDs", "geometry");
-
- // load geometry
- TGeoManager::Import(fileName);
- }*/
- }
-}
-
-Bool_t AliHighMultiplicitySelector::Process(Long64_t entry)
-{
- //
- // processing
- //
-
- if (AliSelectorRL::Process(entry) == kFALSE)
- return kFALSE;
-
- // Check prerequisites
- if (!fESD)
- {
- AliDebug(AliLog::kError, "ESD branch not available");
- return kFALSE;
- }
-
- static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
- Bool_t eventTriggered = triggerAnalysis->IsTriggerBitFired(fESD->GetTriggerMask(), AliTriggerAnalysis::kMB1);
-
- if (!eventTriggered)
- {
- AliDebug(AliLog::kDebug, "Event not triggered. Skipping.");
- return kTRUE;
- }
-
- AliStack* stack = GetStack();
- if (!stack)
- {
- AliDebug(AliLog::kError, "Stack not available");
- return kFALSE;
- }
-
- Int_t nPrim = stack->GetNprimary();
- Int_t multiplicity21 = 0;
- Int_t multiplicity16 = 0;
-
- for (Int_t iMc = 0; iMc < nPrim; ++iMc)
- {
- AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
-
- TParticle* particle = stack->Particle(iMc);
-
- if (!particle)
- {
- AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
- continue;
- }
-
- if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
- continue;
-
- if (fCentralRegion)
- {
- if (TMath::Abs(particle->Eta()) < 1.05)
- multiplicity21++;
- if (TMath::Abs(particle->Eta()) < 0.8)
- multiplicity16++;
- }
- else
- {
- if (TMath::Abs(particle->Eta()) < 2.1)
- multiplicity21++;
- if (TMath::Abs(particle->Eta()) < 1.6)
- multiplicity16++;
- }
- }// end of mc particle
-
- AliRunLoader* runLoader = GetRunLoader();
- if (!runLoader)
- {
- AliDebug(AliLog::kError, "runloader not available");
- return kFALSE;
- }
-
- // TDirectory::TContext restores the current directory is restored when the scope ends.
- // This helps around ROOT bug #26025 and is good behaviour anyway
- TDirectory::TContext context(0);
- AliITSLoader* loader = (AliITSLoader*) runLoader->GetLoader( "ITSLoader" );
- loader->LoadDigits("READ");
- TTree* treeD = loader->TreeD();
- if (!treeD)
- {
- AliDebug(AliLog::kError, "Could not retrieve TreeD of ITS");
- return kFALSE;
- }
-
- treeD->SetBranchStatus("*", 0);
- treeD->SetBranchStatus("ITSDigitsSPD.fTracks*", 1);
- treeD->SetBranchStatus("ITSDigitsSPD.fCoord1", 1);
-
- TClonesArray* digits = 0;
- treeD->SetBranchAddress("ITSDigitsSPD", &digits);
- if (digits)
- digits->Clear();
-
- // each value for both layers
- Int_t totalNumberOfFO[2];
- Int_t chipsHitByPrimaries[2];
- //Int_t chipsHitBySecondaries[2];
-
- for (Int_t i=0; i<2; ++i)
- {
- totalNumberOfFO[i] = 0;
- chipsHitByPrimaries[i] = 0;
- //chipsHitBySecondaries[i] = 0;
- }
-
- Int_t startSPD = 0; //geom->GetStartSPD();
- Int_t lastSPD = 239; //geom->GetLastSPD();
-
- //printf("%d %d\n", startSPD, lastSPD);
-// for (Int_t l=0; l<240; ++l) { AliITSgeomTGeo::GetModuleId(l, i, j, k); printf("%d --> %d\n", l, i); }
-
- // loop over modules (ladders)
- for (Int_t moduleIndex=startSPD; moduleIndex<lastSPD+1; moduleIndex++)
- {
- if (fCentralRegion)
- {
- if ((moduleIndex % 4) == 0 || (moduleIndex % 4) == 3)
- continue;
- }
-
- Int_t currentLayer = 0;
- if (moduleIndex >= 80)
- currentLayer = 1;
-
- treeD->GetEvent(moduleIndex);
-
- // get number of digits in this module
- Int_t ndigitsInModule = digits->GetEntriesFast();
-
- // get number of digits in each chip of the module
- Int_t ndigitsInChip[5];
- Bool_t hitByPrimary[5];
- for( Int_t iChip=0; iChip<5; iChip++)
- {
- ndigitsInChip[iChip]=0;
- hitByPrimary[iChip] = kFALSE;
- }
-
- // loop over digits in this module
- for (Int_t iDig=0; iDig<ndigitsInModule; iDig++)
- {
- AliITSdigitSPD* dp = (AliITSdigitSPD*) digits->At(iDig);
- Int_t column = dp->GetCoord1();
- Int_t isChip = column / 32;
-
- //printf("Digit %d has column %d which translates to chip %d\n", iDig, column, isChip);
-
- fChipsFired->Fill(moduleIndex, isChip);
-
- ndigitsInChip[isChip]++;
-
- Bool_t debug = kFALSE;
-
- // find out which particle caused this chip to fire
- // if we find at least one primary we consider this chip being fired by a primary
- for (Int_t track=0; track<10; ++track)
- {
- Int_t label = dp->GetTrack(track);
-
- if (label < 0)
- continue;
-
- if (debug)
- printf("track %d contains label %d\n", track, label);
-
- TParticle* particle = stack->Particle(label);
-
- if (!particle)
- {
- AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (digit loop).", label));
- continue;
- }
-
- if (debug)
- {
- particle->Print();
- printf("statuscode = %d, p = %f, m = %d\n", particle->GetStatusCode(), particle->P(), particle->GetFirstMother());
- }
-
- // TODO delta electrons should be traced back to their mother. this is e.g. solved in AliITSClusterFinderV2::CheckLabels2
- while (particle->P() < 0.02 && particle->GetStatusCode() == 0 && particle->GetFirstMother() >= 0)
- {
- label = particle->GetFirstMother();
- particle = stack->Particle(label);
-
- if (!particle)
- break;
-
- if (debug)
- {
- printf("-->\n");
- printf("statuscode = %d, p = %f, m = %d\n", particle->GetStatusCode(), particle->P(), particle->GetFirstMother());
- particle->Print();
- }
- }
-
- if (!particle)
- {
- AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (digit loop, finding delta electrons).", label));
- continue;
- }
-
- if (label > nPrim)
- continue;
-
- if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
- continue;
-
- if (debug)
- printf("This was a primary (or delta-electron of a primary)!\n");
-
- hitByPrimary[isChip] = kTRUE;
- }
- }
-
- // get number of FOs in the module
- for (Int_t ifChip=0; ifChip<5; ifChip++)
- if( ndigitsInChip[ifChip] >= 1 )
- {
- totalNumberOfFO[currentLayer]++;
- if (hitByPrimary[ifChip])
- {
- chipsHitByPrimaries[currentLayer]++;
- }
- //else
- // chipsHitBySecondaries[currentLayer]++;
- }
- }
-
- //printf("Fired chips: %d %d\n", totalNumberOfFOLayer1, totalNumberOfFOLayer2);
-
- // now find clusters
- Int_t clustersLayer[2];
- clustersLayer[0] = 0;
- clustersLayer[1] = 0;
-
- loader->LoadRecPoints("READ");
- TTree* treeR = loader->TreeR();
- if (!treeR)
- {
- AliDebug(AliLog::kError, "Could not retrieve TreeR of ITS");
- return kFALSE;
- }
-
- // TODO branches!
- //treeR->SetBranchStatus("*", 0);
-
- TClonesArray* itsClusters = 0;
- treeR->SetBranchAddress("ITSRecPoints", &itsClusters);
-
- Int_t nTreeEntries = treeR->GetEntries();
- for (Int_t iEntry = 0; iEntry < nTreeEntries; ++iEntry)
- {
- if (!treeR->GetEvent(iEntry))
- continue;
-
- Int_t nClusters = itsClusters->GetEntriesFast();
-
- while(nClusters--)
- {
- AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters);
-
- if (cluster->GetLayer() == 0)
- {
- clustersLayer[0]++;
- fClusterZL1->Fill(cluster->GetZ());
- }
- else if (cluster->GetLayer() == 1)
- {
- clustersLayer[1]++;
- fClusterZL2->Fill(cluster->GetZ());
- }
- }
- }
-
- fPrimaryL1->Fill(multiplicity21, totalNumberOfFO[0], chipsHitByPrimaries[0], clustersLayer[0]);
- fPrimaryL2->Fill(multiplicity16, totalNumberOfFO[1], chipsHitByPrimaries[1], clustersLayer[1]);
-
- return kTRUE;
-}
-
-Bool_t AliHighMultiplicitySelector::Notify()
-{
- // get next ITS runloader
-
- AliRunLoader* runLoader = GetRunLoader();
-
- if (runLoader)
- {
- AliITSLoader* loader = (AliITSLoader* )runLoader->GetLoader( "ITSLoader" );
- if (loader)
- {
- loader->UnloadDigits();
- loader->UnloadRecPoints();
- }
- }
-
- return AliSelectorRL::Notify();
-}
-
-void AliHighMultiplicitySelector::SlaveTerminate()
-{
- // The SlaveTerminate() function is called after all entries or objects
- // have been processed. When running with PROOF SlaveTerminate() is called
- // on each slave server.
-
- AliSelectorRL::SlaveTerminate();
-
- // Add the histograms to the output on each slave server
- if (!fOutput)
- {
- AliDebug(AliLog::kError, "ERROR: Output list not initialized.");
- return;
- }
-
- fOutput->Add(fChipsFired);
- fOutput->Add(fPrimaryL1);
- fOutput->Add(fPrimaryL2);
- fOutput->Add(fClusterZL1);
- fOutput->Add(fClusterZL2);
-}
-
-void AliHighMultiplicitySelector::Terminate()
-{
- // The Terminate() function is the last function to be called during
- // a query. It always runs on the client, it can be used to present
- // the results graphically or save the results to file.
-
- AliSelectorRL::Terminate();
-
- fChipsFired = dynamic_cast<TH2F*> (fOutput->FindObject("fChipsFired"));
- fPrimaryL1 = dynamic_cast<TNtuple*> (fOutput->FindObject("fPrimaryL1"));
- fPrimaryL2 = dynamic_cast<TNtuple*> (fOutput->FindObject("fPrimaryL2"));
- fClusterZL1 = dynamic_cast<TH1F*> (fOutput->FindObject("fClusterZL1"));
- fClusterZL2 = dynamic_cast<TH1F*> (fOutput->FindObject("fClusterZL2"));
-
- if (!fClusterZL1 || !fClusterZL2 || !fChipsFired || !fPrimaryL1 || !fPrimaryL2)
- {
- AliError("Histograms not available");
- return;
- }
-
- WriteHistograms();
-}
-
-void AliHighMultiplicitySelector::WriteHistograms(const char* filename)
-{
- // write the histograms to a file
-
- TFile* file = TFile::Open(filename, "RECREATE");
-
- fChipsFired->Write();
- fPrimaryL1->Write();
- fPrimaryL2->Write();
- fClusterZL1->Write();
- fClusterZL2->Write();
-
- file->Close();
-}
-
-void AliHighMultiplicitySelector::ReadHistograms(const char* filename)
-{
- // read the data from a file and fill histograms
-
- TFile* file = TFile::Open(filename);
-
- if (!file)
- return;
-
- fPrimaryL1 = dynamic_cast<TNtuple*> (file->Get("fPrimaryL1"));
- fPrimaryL2 = dynamic_cast<TNtuple*> (file->Get("fPrimaryL2"));
- 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
- #define BINNING_LAYER2 801, -0.5, 800.5
-
- fChipsLayer1 = new TH1F("fChipsLayer1", "Layer 1;Fired Chips;Count", BINNING_LAYER1);
- fChipsLayer2 = new TH1F("fChipsLayer2", "Layer 2;Fired Chips;Count", BINNING_LAYER2);
-
- fL1vsL2 = new TH2F("fL1vsL2", ";Fired Chips Layer 1;Fired Chips Layer 2", BINNING_LAYER1, BINNING_LAYER2);
- fMvsL1 = new TH2F("fMvsL1", ";true multiplicity;fired chips layer1", MULT, BINNING_LAYER1);
- fMvsL2 = new TH2F("fMvsL2", ";true multiplicity;fired chips layer2", MULT, BINNING_LAYER2);
-
- fClvsL1 = new TH2F("fClvsL1", ";clusters layer1;fired chips layer1", MULT, BINNING_LAYER1);
- fClvsL2 = new TH2F("fClvsL2", ";clusters layer2;fired chips layer2", MULT, BINNING_LAYER2);
-
- for (Int_t i = 0; i < fPrimaryL1->GetEntries(); i++)
- {
- fPrimaryL1->GetEvent(i);
- fPrimaryL2->GetEvent(i);
-
- Int_t multiplicity21 = (Int_t) fPrimaryL1->GetArgs()[0];
- Int_t multiplicity16 = (Int_t) fPrimaryL2->GetArgs()[0];
-
- Int_t totalNumberOfFO[2];
- totalNumberOfFO[0] = (Int_t) fPrimaryL1->GetArgs()[1];
- totalNumberOfFO[1] = (Int_t) fPrimaryL2->GetArgs()[1];
-
- Int_t chipsHitByPrimaries[2];
- chipsHitByPrimaries[0] = (Int_t) fPrimaryL1->GetArgs()[2];
- chipsHitByPrimaries[1] = (Int_t) fPrimaryL2->GetArgs()[2];
-
- Int_t clustersLayer[2];
- clustersLayer[0] = (Int_t) fPrimaryL1->GetArgs()[3];
- clustersLayer[1] = (Int_t) fPrimaryL2->GetArgs()[3];
-
- fChipsLayer1->Fill(totalNumberOfFO[0]);
- fChipsLayer2->Fill(totalNumberOfFO[1]);
-
- fL1vsL2->Fill(totalNumberOfFO[0], totalNumberOfFO[1]);
-
- fMvsL1->Fill(multiplicity21, totalNumberOfFO[0]);
- fMvsL2->Fill(multiplicity16, totalNumberOfFO[1]);
-
- fClvsL1->Fill(clustersLayer[0], totalNumberOfFO[0]);
- fClvsL2->Fill(clustersLayer[1], totalNumberOfFO[1]);
- }
-}
-
-TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut, Int_t upperCut)
-{
- //
- // returns the trigger efficiency as function of multiplicity with a given cut
- //
-
- //cut and multiply with x-section
- TH1* allEvents = multVsLayer->ProjectionX("fMvsL_x_total", 1, 1001);
- //allEvents->Sumw2();
-
- //cut and multiply with x-section
- TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, upperCut);
- //proj->Sumw2();
-
- //new TCanvas; allEvents->DrawCopy(); gPad->SetLogy();
- //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
-
- // make probability distribution out of it
- // TODO binomial errors do not work??? weird...
- proj->Divide(proj, allEvents, 1, 1, "B");
-
- return proj;
-}
-
-TH1* AliHighMultiplicitySelector::GetXSectionCut(TH1* xSection, TH2* multVsLayer, Int_t cut)
-{
- // returns the rel. cross section of the true spectrum that is measured when a cut at <cut> is performed
-
- TH1* proj = GetTriggerEfficiency(multVsLayer, cut);
-
- //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
-
- for (Int_t i=1; i<=proj->GetNbinsX(); i++)
- {
- if (i <= xSection->GetNbinsX())
- {
- Double_t value = proj->GetBinContent(i) * xSection->GetBinContent(i);
- Double_t error = 0;
-
- if (value != 0)
- error = value * (proj->GetBinError(i) / proj->GetBinContent(i) + xSection->GetBinError(i) / xSection->GetBinContent(i));
-
- proj->SetBinContent(i, value);
- proj->SetBinError(i, error);
- }
- else
- {
- proj->SetBinContent(i, 0);
- proj->SetBinError(i, 0);
- }
- }
-
- //new TCanvas; proj->DrawCopy(); gPad->SetLogy();
-
- return proj;
-}
-
-void AliHighMultiplicitySelector::MakeGraphs2(const char* title, TH1* xSection, TH2* fMvsL)
-{
- // creates graphs
-
- TGraph* effGraph = new TGraph;
- effGraph->SetTitle(Form("%s;Cut on fired chips;mult. where eff. >50%%", title));
- TGraph* ratioGraph = new TGraph;
- ratioGraph->SetTitle(Form("%s;Cut on fired chips;x-section_(>=eff. limit) / x-section_(total)", title));
- TGraph* totalGraph = new TGraph;
- totalGraph->SetTitle(Form("%s;Cut on fired chips;rel x-section_(>=eff. limit)", title));
-
- for (Int_t cut = 0; cut <= 300; cut+=50)
- {
- TH1* proj = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("clone");
-
- //proj->Rebin(3);
- //proj->Scale(1.0 / 3);
-
- new TCanvas; proj->DrawCopy();
-
- Int_t limitBin = proj->GetNbinsX()+1;
- while (limitBin > 1 && proj->GetBinContent(limitBin-1) > 0.5)
- limitBin--;
-
- Float_t limit = proj->GetXaxis()->GetBinCenter(limitBin);
-
- effGraph->SetPoint(effGraph->GetN(), cut, limit);
-
- proj = GetXSectionCut(xSection, fMvsL, cut);
-
- Double_t ratio = 0;
- Double_t total = 0;
- if (proj->Integral(1, 1001) > 0)
- {
- ratio = proj->Integral(proj->FindBin(limit), 1001) / proj->Integral(1, 1001);
- total = proj->Integral(proj->FindBin(limit), 1001);
- }
-
- ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio);
- totalGraph->SetPoint(totalGraph->GetN(), cut, total);
-
- 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);
- }
-
- TCanvas* canvas = new TCanvas(Form("%s_Efficiency", title), Form("%s_Efficiency", title), 1200, 800);
- canvas->Divide(2, 2);
-
- canvas->cd(1);
- effGraph->Draw("A*");
-
- for (Int_t i=8; i<=10; ++i)
- {
- TLine* line = new TLine(0, xSection->GetMean() * i, 300, xSection->GetMean() * i);
- line->Draw();
- }
-
- canvas->cd(2); ratioGraph->Draw("A*");
- canvas->cd(3); gPad->SetLogy(); totalGraph->Draw("A*");
-
- canvas->SaveAs(Form("%s.gif", canvas->GetName()));
-}
-
-void AliHighMultiplicitySelector::MakeGraphs(const char* title, TH1* xSection, TH2* fMvsL, Int_t limit)
-{
- // relative x-section, once we have a collision
-
- xSection->Scale(1.0 / xSection->Integral());
-
- TGraph* ratioGraph = new TGraph;
- ratioGraph->SetTitle(Form("%s;Cut on fired chips;x-section_(>=%d) / x-section_(total)", title, limit));
- TGraph* totalGraph = new TGraph;
- totalGraph->SetTitle(Form("%s;Cut on fired chips;rel x-section_(>=%d)", title, limit));
-
- Double_t max = 0;
- Int_t bestCut = -1;
- Double_t bestRatio = -1;
- Double_t bestTotal = -1;
- Int_t fullCut = -1;
- Double_t fullRatio = -1;
- Double_t fullTotal = -1;
-
- fMvsL->Sumw2();
-
- for (Int_t cut = 50; cut <= 300; cut+=2)
- {
- TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
-
- Double_t ratio = 0;
- Double_t total = 0;
- if (proj->Integral(1, 1000) > 0)
- {
- ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000);
- total = proj->Integral(limit, 1000);
- }
-
- max = TMath::Max(max, total);
-
- //printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio);
-
- if (total < max * 0.9 && bestCut == -1)
- {
- bestCut = cut;
- bestRatio = ratio;
- bestTotal = total;
- }
-
- if (ratio == 1 && fullCut == -1)
- {
- fullCut = cut;
- fullRatio = ratio;
- fullTotal = total;
- }
-
- ratioGraph->SetPoint(ratioGraph->GetN(), cut, ratio);
- totalGraph->SetPoint(totalGraph->GetN(), cut, total);
- }
-
- if (bestCut != -1)
- printf("Best cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", bestCut, limit, bestTotal, limit, bestRatio);
- if (fullCut != -1)
- printf("100%% cut at %d: rel. x-section_(>=%d) = %e %%; x-section_(>=%d) / x-section_(total) = %f %%\n", fullCut, limit, fullTotal, limit, fullRatio);
-
- TCanvas* canvas = new TCanvas(Form("%s_RatioXSection_%d", title, limit), Form("%s_RatioXSection_%d", title, limit), 600, 400);
- ratioGraph->Draw("A*");
- canvas->SaveAs(Form("%s.gif", canvas->GetName()));
-
- canvas = new TCanvas(Form("%s_TotalXSection_%d", title, limit), Form("%s_TotalXSection_%d", title, limit), 600, 400);
- totalGraph->Draw("A*");
- canvas->SaveAs(Form("%s.gif", canvas->GetName()));
-}
-
-void AliHighMultiplicitySelector::JPRPlots()
-{
- // more plots
-
- /*
-
- gSystem->Load("libPWG0base");
- .L AliHighMultiplicitySelector.cxx+
- x = new AliHighMultiplicitySelector();
- x->ReadHistograms("highmult_hijing100k.root");
- x->JPRPlots();
-
- */
-
- // 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"));
-
- for (Int_t i=0; i<2; ++i)
- {
- if (!xSections[i])
- continue;
-
- TH1* xSection = xSections[i];
- TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
- //Int_t cut = (i == 0) ? 164 : 150; // 8 times the mean
- //Int_t cut = (i == 0) ? 178 : 166; // 9 times the mean
- Int_t cut = (i == 0) ? 190 : 182; // 10 times the mean
-
- // limit is N times the mean
- Int_t limit = (Int_t) (xSection->GetMean() * 10);
-
- // 10^28 lum --> 1.2 kHz
- // 10^31 lum --> 1200 kHz
- Float_t rate = 1200e3;
-
- // time in s
- Float_t lengthRun = 1e6;
-
- 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) ? 400 : 350);
- //xSection->GetYaxis()->SetTitle("relative cross section");
- xSection->GetYaxis()->SetTitleOffset(1.2);
-
- // relative x-section, once we have a collision
- xSection->Scale(1.0 / xSection->Integral());
-
- TH1* proj = GetXSectionCut(xSection, fMvsL, cut);
-
- Double_t ratio = 0;
- Double_t total = 0;
- if (proj->Integral(1, 1000) > 0)
- {
- ratio = proj->Integral(limit, 1000) / proj->Integral(1, 1000);
- total = proj->Integral(limit, 1000);
- }
-
- printf("Cut at %d: rel. x-section_(>=%d) = %e; x-section_(>=%d) / x-section_(total) = %f\n", cut, limit, total, limit, ratio);
-
- TCanvas* canvas = new TCanvas(Form("HMPlots_%d", i), Form("HMPlots_%d", i), 800, 600);
- canvas->SetLogy();
- xSection->DrawCopy();
- proj->SetLineColor(2);
- proj->SetStats(kFALSE);
- proj->DrawCopy("SAME");
-
- TLegend* legend = new TLegend(0.15, 0.15, 0.45, 0.3);
- legend->SetFillColor(0);
- legend->AddEntry(xSection, "no trigger");
- legend->AddEntry(proj, Form("FO trigger > %d chips", cut));
- legend->Draw();
-
- TLine* line = new TLine(limit, xSection->GetMinimum() * 0.5, limit, xSection->GetMaximum() * 2);
- line->SetLineWidth(2);
- line->Draw();
-
- canvas->SaveAs(Form("%s.gif", canvas->GetName()));
-
- TCanvas* canvas2 = new TCanvas(Form("HMPlots_%d_Random", i), Form("HMPlots_%d_Random", i), 800, 600);
- //canvas2->SetTopMargin(0.05);
- //canvas2->SetRightMargin(0.05);
- canvas2->SetLogy();
- xSection->DrawCopy("HIST");
-
- TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.3);
- legend2->SetFillColor(0);
- legend2->AddEntry(xSection, "no trigger");
-
- TH1* proj2 = (TH1*) proj->Clone("random");
- proj2->Reset();
- // MB lengthRun s 100 Hz
- Int_t nTrigger = (Int_t) (100 * lengthRun * proj->Integral(1, 1000));
- proj2->FillRandom(proj, nTrigger);
-
- TH1* proj3 = (TH1*) proj2->Clone("random_clone");
- proj3->Divide(proj);
- proj3->Fit("pol0", "0", "");
- proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0));
-
- /*
- proj2->DrawCopy("SAME");
- legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio)));
- */
-
- proj2 = (TH1*) proj->Clone("random2");
- proj2->Reset();
- // 10^31 lum --> 1200 kHz; lengthRun s
- nTrigger = (Int_t) (rate * proj->Integral(1, 1000) * lengthRun);
- proj2->FillRandom(proj, nTrigger);
-
- proj3 = (TH1*) proj2->Clone("random_clone2");
- proj3->Divide(proj);
- proj3->Fit("pol0", "0", "");
- proj2->Scale(1.0 / proj3->GetFunction("pol0")->GetParameter(0));
-
- proj2->SetLineColor(4);
- proj2->SetMarkerStyle(7);
- proj2->SetMarkerColor(4);
- proj2->DrawCopy("SAME P");
- //legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio)));
- legend2->AddEntry(proj2, Form("FO trigger > %d chips", cut));
-
- legend2->Draw();
- line->Draw();
-
- canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
- canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
- }
-}
-
-void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
-{
- //
- // 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->Ntrigger();
-
- gSystem->Load("libPWG0base");
- .L AliHighMultiplicitySelector.cxx+g
- x = new AliHighMultiplicitySelector();
- x->ReadHistograms("highmult_hijing100k.root");
- x->Ntrigger(kFALSE);
-
- */
-
- // 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"));
-
- // 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 = 1e5;
-
- 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;
-
- TH1* xSection = xSections[i];
- TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
-
- //Int_t nCuts = 6;
- //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 };
- //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 };
-
- 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()->SetTitle("relative cross section");
- xSection->GetYaxis()->SetTitleOffset(1.2);
-
- // relative x-section, once we have a collision
- xSection->Scale(1.0 / xSection->Integral());
-
- TCanvas* canvas2 = new TCanvas(Form("HMPlots2_%d_Random", i), Form("HMPlots2_%d_Random", i), 800, 600);
- canvas2->SetTopMargin(0.05);
- canvas2->SetRightMargin(0.05);
- canvas2->SetLogy();
-
- if (relative)
- xSection->DrawCopy("HIST");
-
- TLegend* legend2 = new TLegend(0.15, 0.15, 0.6, 0.4);
- legend2->SetFillColor(0);
- legend2->AddEntry(xSection, "cross-section");
-
- for (Int_t currentCut = 0; currentCut<nCuts; ++currentCut)
- {
- 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);
-
- printf("Cut at %d: rel. x-section = %e\n", cut, total);
-
- TH1* proj2 = (TH1*) proj->Clone("random2");
- proj2->Reset();
-
- // calculate downscale factor
- Float_t normalRate = rate * proj->Integral(1, 1000);
- Float_t downScale = normalRate / ratePerTrigger[currentCut];
- 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;
- for (Int_t bin=1; bin<proj2->GetNbinsX(); ++bin)
- if (proj2->GetBinContent(bin) < 5)
- {
- removed += (Int_t) proj2->GetBinContent(bin);
- proj2->SetBinContent(bin, 0);
- }
-
- Printf("Removed %d events", removed);
-
- if (relative)
- proj2->Scale(1.0 / nTrigger * proj->Integral(1, 1000));
-
- proj2->SetLineColor(colors[currentCut]);
- proj2->SetMarkerStyle(markers[currentCut]);
- proj2->SetMarkerColor(colors[currentCut]);
-
- 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();
-
- canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
- canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
- }
-}
-
-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
-
- /*
-
- gSystem->Load("libPWG0base");
- .L AliHighMultiplicitySelector.cxx+
- x = new AliHighMultiplicitySelector();
- x->ReadHistograms("highmult_pythia.root");
- x->DrawHistograms();
-
-
- gSystem->Load("libPWG0base");
- .L AliHighMultiplicitySelector.cxx+
- x = new AliHighMultiplicitySelector();
- x->ReadHistograms("highmult_hijing.root");
- x->DrawHistograms();
-
- gSystem->Load("libPWG0base");
- .L AliHighMultiplicitySelector.cxx+
- x = new AliHighMultiplicitySelector();
- x->ReadHistograms("highmult_central.root");
- x->DrawHistograms();
-
- gSystem->Load("libANALYSIS");
- gSystem->Load("libPWG0base");
- .L AliHighMultiplicitySelector.cxx++
- x = new AliHighMultiplicitySelector();
- x->ReadHistograms("highmult_hijing100k.root");
- x->DrawHistograms();
-
- */
-
- /*TCanvas* canvas = new TCanvas("chips", "chips", 600, 400);
-
- fChipsLayer2->SetLineColor(2);
- fChipsLayer2->SetStats(kFALSE);
- fChipsLayer1->SetStats(kFALSE);
- fChipsLayer2->SetTitle("");
- fChipsLayer2->DrawCopy();
- fChipsLayer1->DrawCopy("SAME");
- canvas->SaveAs("chips.gif");
-
- canvas = new TCanvas("L1vsL2", "L1vsL2", 600, 400);
- fL1vsL2->SetStats(kFALSE);
- fL1vsL2->DrawCopy("COLZ");
- gPad->SetLogz();
- canvas->SaveAs("L1vsL2.gif");*/
-
- TCanvas *canvas = new TCanvas("L1", "L1", 800, 600);
- canvas->SetTopMargin(0.05);
- canvas->SetRightMargin(0.12);
- fMvsL1->SetStats(kFALSE);
- fMvsL1->DrawCopy("COLZ");
- gPad->SetLogz();
-
- 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);
- func->DrawCopy("SAME");
-
- canvas->SaveAs("L1.gif");
-
- canvas = new TCanvas("L2", "L2", 600, 400);
- //fMvsL2->GetYaxis()->SetRangeUser(0, 150);
- fMvsL2->SetStats(kFALSE);
- fMvsL2->DrawCopy("COLZ");
- gPad->SetLogz();
- func->SetParameter(0, 800-5*4);
- func->DrawCopy("SAME");
- canvas->SaveAs("L2.gif");
-
- // get x-sections
- TFile* file = TFile::Open("crosssectionEx.root");
- if (file)
- {
- TH1* xSection2 = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
- TH1* xSection15 = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
-
- MakeGraphs2("Layer1", xSection2, fMvsL1);
- return;
-
- // 5 times the mean
- //MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 5)); //75 * 2 * 2);
- //MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 5)); //(Int_t) (75 * 1.5 * 2));
-
- MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 8));
- MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 8));
- MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 9));
- MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 9));
- MakeGraphs("Layer1", xSection2, fMvsL1, (Int_t) (xSection2->GetMean() * 10));
- MakeGraphs("Layer2", xSection15, fMvsL2, (Int_t) (xSection15->GetMean() * 10));
-
- file->Close();
- }
-
- // make spread hists
- TGraph* spread = new TGraph;
- spread->SetTitle("Spread L1;true multiplicity;RMS");
-
- for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
- {
- TH1* proj = fMvsL1->ProjectionY("proj", i, i);
- spread->SetPoint(spread->GetN(), i, proj->GetRMS());
- }
-
- canvas = new TCanvas("SpreadL1", "SpreadL1", 600, 400);
- spread->Draw("A*");
- canvas->SaveAs(Form("%s.gif", canvas->GetName()));
-
- TF1* log = new TF1("log", "[0]*log([1]*x)", 1, 150);
- log->SetParLimits(0, 0, 10);
- log->SetParLimits(1, 1e-5, 10);
-
- spread->Fit(log, "", "", 1, 150);
- log->DrawCopy("SAME");
-
- TGraph* spread2 = new TGraph;
- spread2->SetTitle("Spread L2;true multiplicity;RMS");
-
- for (Int_t i=1; i<=fMvsL1->GetNbinsX(); ++i)
- {
- TH1* proj = fMvsL2->ProjectionY("proj", i, i);
- spread2->SetPoint(spread2->GetN(), i, proj->GetRMS());
- }
-
- canvas = new TCanvas("SpreadL2", "SpreadL2", 600, 400);
- spread2->Draw("A*");
- canvas->SaveAs(Form("%s.gif", canvas->GetName()));
-
- spread2->Fit(log, "", "", 1, 150);
- log->DrawCopy("SAME");
-
- canvas = new TCanvas("Clusters_L1", "Clusters_L1", 600, 400);
- fClvsL1->SetStats(kFALSE);
- fClvsL1->DrawCopy("COLZ");
- gPad->SetLogz();
-
- func->SetParameter(0, 400-5*2);
- func->DrawCopy("SAME");
-
- canvas->SaveAs("Clusters_L1.gif");
-
- canvas = new TCanvas("Clusters_L2", "Clusters_L2", 600, 400);
- fClvsL2->SetStats(kFALSE);
- fClvsL2->DrawCopy("COLZ");
- gPad->SetLogz();
- func->SetParameter(0, 800-5*4);
- func->DrawCopy("SAME");
- canvas->SaveAs("Clusters_L2.gif");
-
- canvas = new TCanvas("ChipsFired", "ChipsFired", 600, 400);
- //fChipsFired->GetYaxis()->SetRangeUser(0, 150);
- fChipsFired->SetStats(kFALSE);
- fChipsFired->DrawCopy("COLZ");
- canvas->SaveAs("ChipsFired.gif");
-
- /*TH1F* tresholdHistL1 = new TH1F("tresholdHistL1", ";chip treshold;<n>", BINNING_LAYER1);
- TH1F* tresholdHistL2 = new TH1F("tresholdHistL2", ";chip treshold;<n>", BINNING_LAYER2);
-
- for (Int_t treshold = 0; treshold < 800; treshold++)
- {
- if (fPrimaryL1->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
- {
- TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL1->GetHistogram());
- if (mult)
- tresholdHistL1->Fill(treshold, mult->GetMean());
- }
- if (fPrimaryL2->Draw("multiplicity>>mult", Form("firedChips>%d", treshold), "goff") > 0)
- {
- TH1F* mult = dynamic_cast<TH1F*> (fPrimaryL2->GetHistogram());
- if (mult)
- tresholdHistL2->Fill(treshold, mult->GetMean());
- }
- }
-
- canvas = new TCanvas("TresholdL1", "TresholdL1", 600, 400);
- tresholdHistL1->Draw();
- canvas->SaveAs(Form("%s.gif", canvas->GetName()));
-
- canvas = new TCanvas("TresholdL2", "TresholdL2", 600, 400);
- tresholdHistL2->Draw();
- canvas->SaveAs(Form("%s.gif", canvas->GetName()));*/
-
- fPrimaryL1->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff1", "", "prof goff");
- fPrimaryL2->Draw("(chipsByPrimaries/firedChips):multiplicity>>eff2", "", "prof goff");
-
- canvas = new TCanvas("Efficiency", "Efficiency", 600, 400);
- fPrimaryL1->GetHistogram()->SetStats(kFALSE);
- fPrimaryL1->GetHistogram()->Draw();
- fPrimaryL2->GetHistogram()->SetLineColor(2);
- fPrimaryL2->GetHistogram()->SetStats(kFALSE);
- fPrimaryL2->GetHistogram()->Draw("SAME");
- canvas->SaveAs(Form("%s.gif", canvas->GetName()));
-
- canvas = new TCanvas("ClustersZL1", "ClustersZL1", 600, 400);
- fClusterZL1->Rebin(2);
- fClusterZL1->Draw();
- canvas->SaveAs(Form("%s.gif", canvas->GetName()));
-
- canvas = new TCanvas("ClustersZL2", "ClustersZL2", 600, 400);
- fClusterZL2->Draw();
- fClusterZL2->Rebin(2);
- canvas->SaveAs(Form("%s.gif", canvas->GetName()));
-}
-
-TGraph* AliHighMultiplicitySelector::IntFractRate()
-{
- // A plot which shows the fractional rate above any threshold
- // as function of threshold (i.e. the integral of dSigma/dN as function of
- // N, normalised to 1 for N=0)
-
- /*
- gSystem->Load("libANALYSIS");
- gSystem->Load("libPWG0base");
- .L AliHighMultiplicitySelector.cxx+
- x = new AliHighMultiplicitySelector();
- x->ReadHistograms("highmult_hijing100k.root");
- x->IntFractRate();
- */
-
- // get x-sections
- TFile* file = TFile::Open("crosssectionEx.root");
- if (!file)
- return 0;
-
- TH1* xSection;
- xSection = static_cast<TH1*> (gFile->Get("xSection2Ex"));
-
- TGraph* result = new TGraph;
-
- for (Int_t threshold = 0; threshold < 300; threshold += 2)
- {
- TH1* proj = GetXSectionCut(xSection, fMvsL1, threshold);
-
- //new TCanvas; proj->DrawCopy();
-
- Double_t integral = proj->Integral();
-
- Printf("Cut at %d, integral is %e", threshold, integral);
-
- result->SetPoint(result->GetN(), threshold, integral);
- }
-
- TCanvas* canvas = new TCanvas("IntFractRate", "IntFractRate", 600, 400);
- gPad->SetLogy();
-
- result->Draw("A*");
- result->GetXaxis()->SetTitle("threshold");
- result->GetYaxis()->SetTitle("integrated fractional rate above threshold");
-
- canvas->SaveAs(Form("%s.gif", canvas->GetName()));
-
- 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 / commitdiff