adding selector for high multiplicity trigger study
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 18 May 2007 11:49:19 +0000 (11:49 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 18 May 2007 11:49:19 +0000 (11:49 +0000)
PWG0/highMultiplicity/AliHighMultiplicitySelector.cxx [new file with mode: 0644]
PWG0/highMultiplicity/AliHighMultiplicitySelector.h [new file with mode: 0644]
PWG0/highMultiplicity/extrapolateXSection.C [new file with mode: 0644]
PWG0/highMultiplicity/runHighMultiplicitySelector.C [new file with mode: 0644]

diff --git a/PWG0/highMultiplicity/AliHighMultiplicitySelector.cxx b/PWG0/highMultiplicity/AliHighMultiplicitySelector.cxx
new file mode 100644 (file)
index 0000000..8c45e12
--- /dev/null
@@ -0,0 +1,957 @@
+/* $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 <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"
+
+//
+//
+
+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),
+  centralRegion(kTRUE)
+{
+  //
+  // 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)
+{
+  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;
+  }
+
+  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, AliPWG0Helper::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 (centralRegion)
+    {
+      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 (centralRegion)
+    {
+      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()
+{
+  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)
+{
+  TFile* file = TFile::Open(filename, "RECREATE");
+
+  fChipsFired->Write();
+  fPrimaryL1->Write();
+  fPrimaryL2->Write();
+  fClusterZL1->Write();
+  fClusterZL2->Write();
+
+  file->Close();
+}
+
+void AliHighMultiplicitySelector::ReadHistograms(const char* filename)
+{
+  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"));
+
+  #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::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
+
+  //cut and multiply with x-section
+  TH1* allEvents = multVsLayer->ProjectionX("fMvsL_x_total", 1, 1000);
+  //allEvents->Sumw2();
+
+  //cut and multiply with x-section
+  TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, 1000);
+  //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");
+
+  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::MakeGraphs(const char* title, TH1* xSection, TH2* fMvsL, Int_t limit)
+{
+  if (!xSection)
+    return;
+
+  // 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()
+{
+  /*
+
+  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
+
+    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);
+
+    // limit is N times the mean
+    Int_t limit = (Int_t) (xSection->GetMean() * 10);
+
+    // 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->SetLogy();
+    xSection->DrawCopy();
+
+    TLegend* legend2 = new TLegend(0.15, 0.15, 0.5, 0.3);
+    legend2->SetFillColor(0);
+    legend2->AddEntry(xSection, "no trigger");
+
+    TH1* proj2 = (TH1*) proj->Clone("random");
+    proj2->Reset();
+    // MB 10^5 s 100 Hz
+    Int_t nTrigger = (Int_t) (100 * 1e5 * 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; 1e5 s
+    nTrigger = (Int_t) (12e5 * proj->Integral(1, 1000) * 1e5);
+    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->DrawCopy("SAME");
+    legend2->AddEntry(proj2, Form("%d evts, FO > %d chips (%d evts)", nTrigger, cut, (Int_t) (nTrigger * ratio)));
+
+    legend2->Draw();
+    line->Draw();
+
+    canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
+  }
+}
+
+void AliHighMultiplicitySelector::DrawHistograms()
+{
+  /*
+
+  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("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", 600, 400);
+  fMvsL1->SetStats(kFALSE);
+  fMvsL1->DrawCopy("COLZ");
+  gPad->SetLogz();
+
+  canvas->SaveAs("L1NoCurve.gif");
+
+  // 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");
+
+  // 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");
+
+  // 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"));
+
+    // 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();
+  }
+
+  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()));
+}
diff --git a/PWG0/highMultiplicity/AliHighMultiplicitySelector.h b/PWG0/highMultiplicity/AliHighMultiplicitySelector.h
new file mode 100644 (file)
index 0000000..85c4394
--- /dev/null
@@ -0,0 +1,62 @@
+/* $Id$ */
+
+#ifndef AliHighMultiplicitySelector_H
+#define AliHighMultiplicitySelector_H
+
+#include "AliSelectorRL.h"
+
+class TH1F;
+class TH2F;
+class TH1;
+class TH2;
+class TNtuple;
+
+class AliHighMultiplicitySelector : public AliSelectorRL {
+  public:
+    AliHighMultiplicitySelector();
+    virtual ~AliHighMultiplicitySelector();
+
+    virtual void    Init(TTree *tree);
+    virtual void    SlaveBegin(TTree *tree);
+    virtual Bool_t  Process(Long64_t entry);
+    virtual Bool_t  Notify();
+    virtual void    SlaveTerminate();
+    virtual void    Terminate();
+
+    void WriteHistograms(const char* filename = "highmult.root");
+    void ReadHistograms(const char* filename = "highmult.root");
+    void DrawHistograms();
+    void JPRPlots();
+
+ protected:
+    void MakeGraphs(const char* title, TH1* xSection, TH2* fMvsL, Int_t limit);
+    TH1* GetXSectionCut(TH1* xSection, TH2* multVsLayer, Int_t cut);
+
+    TH1F* fChipsLayer1;
+    TH1F* fChipsLayer2;
+
+    TH2F* fL1vsL2;
+    TH2F* fMvsL1;
+    TH2F* fMvsL2;
+
+    TH2F* fChipsFired;
+
+    TNtuple* fPrimaryL1;  // multiplicity vs. number of fired chips vs. number of chips fired by primaries
+    TNtuple* fPrimaryL2;  // multiplicity vs. number of fired chips vs. number of chips fired by primaries
+
+    TH1F*    fClusterZL1;
+    TH1F*    fClusterZL2;
+
+    TH2F* fClvsL1;
+    TH2F* fClvsL2;
+
+    Bool_t centralRegion;
+
+ private:
+    AliHighMultiplicitySelector(const AliHighMultiplicitySelector&);
+    AliHighMultiplicitySelector& operator=(const AliHighMultiplicitySelector&);
+
+  ClassDef(AliHighMultiplicitySelector, 0);
+};
+
+#endif
diff --git a/PWG0/highMultiplicity/extrapolateXSection.C b/PWG0/highMultiplicity/extrapolateXSection.C
new file mode 100644 (file)
index 0000000..312b772
--- /dev/null
@@ -0,0 +1,65 @@
+void extrapolateXSection()
+{
+  TFile* file = TFile::Open("crosssection.root");
+  TH1* xSection2 = dynamic_cast<TH1*> (gFile->Get("xSection2"));
+  TH1* xSection15 = dynamic_cast<TH1*> (gFile->Get("xSection15"));
+
+  gROOT->cd();
+
+  TH1F* xSection2Ex = new TH1F("xSection2Ex", ";Npart", 1001, -0.5, 1000.5);
+  TH1F* xSection15Ex = new TH1F("xSection15Ex", ";Npart", 1001, -0.5, 1000.5);
+
+  new TCanvas;
+  xSection2->Draw();
+  xSection2->Fit("expo", "", "", 200, 250);
+  gPad->SetLogy();
+
+  for (Int_t i=1; i<=1000; ++i)
+  {
+    if (i < 250)
+    {
+      xSection2Ex->SetBinContent(i, xSection2->GetBinContent(i));
+      xSection2Ex->SetBinError(i, xSection2->GetBinError(i));
+    }
+    else
+      xSection2Ex->SetBinContent(i, xSection2->GetFunction("expo")->Eval(i));
+  }
+
+  new TCanvas;
+  xSection2Ex->Draw();
+  xSection2->SetLineColor(2);
+  xSection2->Draw("SAME");
+  gPad->SetLogy();
+
+  new TCanvas;
+  xSection15->Draw();
+  xSection15->Fit("expo", "", "", 145, 250);
+  gPad->SetLogy();
+
+  for (Int_t i=1; i<=1000; ++i)
+  {
+    if (i < 145)
+    {
+      xSection15Ex->SetBinContent(i,xSection15->GetBinContent(i));
+      xSection15Ex->SetBinError(i, xSection15->GetBinError(i));
+    }
+    else
+      xSection15Ex->SetBinContent(i, xSection15->GetFunction("expo")->Eval(i));
+  }
+
+  new TCanvas;
+  xSection15Ex->Draw();
+  xSection15->SetLineColor(2);
+  xSection15->Draw("SAME");
+  gPad->SetLogy();
+
+  TFile* file2 = TFile::Open("crosssectionEx.root", "RECREATE");
+
+  xSection2Ex->Write();
+  xSection15Ex->Write();
+
+  file2->Close();
+}
+
+
+
diff --git a/PWG0/highMultiplicity/runHighMultiplicitySelector.C b/PWG0/highMultiplicity/runHighMultiplicitySelector.C
new file mode 100644 (file)
index 0000000..7a7bd4a
--- /dev/null
@@ -0,0 +1,44 @@
+/* $Id$ */
+
+//
+// Script to run the AliFirstPlotsSelector
+//
+
+#include "../CreateESDChain.C"
+#include "../PWG0Helper.C"
+
+void runHighMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE,  const char* option = "", const char* proofServer = "jgrosseo@lxb6046")
+{
+  if (aProof)
+  {
+    connectProof(proofServer);
+    gProof->AddInput(new TParameter<long>("PROOF_MaxSlavesPerNode", (long)2));
+    gProof->AddInput(new TNamed("PROOF_Packetizer", "TAdaptivePacketizer"));
+  }
+
+  //gProof->SetParallel(1);
+
+  TString libraries("libPWG0base");
+  TString packages("PWG0base");
+
+  if (!prepareQuery(libraries, packages, 1))
+    return;
+
+  if (aProof)
+    ProofAddAliRootIncludePath(1, "ITS");
+
+  TChain* chain = CreateESDChain(data, nRuns, offset);
+
+  TList inputList;
+
+  TString selectorName = "AliHighMultiplicitySelector";
+  AliLog::SetClassDebugLevel(selectorName, AliLog::kInfo);
+
+  selectorName += ".cxx+";
+
+  if (aDebug != kFALSE)
+    selectorName += "+g";
+
+  executeQuery(chain, &inputList, selectorName, option);
+}
+