]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG0/highMultiplicity/AliHighMultiplicitySelector.cxx
adding offline triggers
[u/mrichter/AliRoot.git] / PWG0 / highMultiplicity / AliHighMultiplicitySelector.cxx
index 8e5eba028812898c9188aba619d353504ca1dd65..49a1bfa913cb07908fdae8e5e9ef1993e4a97bcd 100644 (file)
@@ -16,6 +16,7 @@
 #include <TGraph.h>
 #include <TLegend.h>
 #include <TLine.h>
+#include <TMath.h>
 
 #include <AliLog.h>
 #include <AliESD.h>
@@ -122,7 +123,7 @@ Bool_t AliHighMultiplicitySelector::Process(Long64_t entry)
     return kFALSE;
   }
 
-  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, AliPWG0Helper::kMB1);
+  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB1);
 
   if (!eventTriggered)
   {
@@ -197,7 +198,7 @@ Bool_t AliHighMultiplicitySelector::Process(Long64_t entry)
 
   TClonesArray* digits = 0;
   treeD->SetBranchAddress("ITSDigitsSPD", &digits);
-  if (digits);
+  if (digits)
     digits->Clear();
 
   // each value for both layers
@@ -863,6 +864,7 @@ void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
 
   /*
 
+  gSystem->Load("libANALYSIS");
   gSystem->Load("libPWG0base");
   .L AliHighMultiplicitySelector.cxx+g
   x = new AliHighMultiplicitySelector();
@@ -886,13 +888,13 @@ void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
   xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
   xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
 
-  // 10^28 lum --> 1.2 kHz
-  // 10^31 lum --> 1200 kHz
-  //Float_t rate = 1200e3;
-  Float_t rate = 250e3;
+  // 10^28 lum --> 1.4 kHz
+  // 10^31 lum --> 1400 kHz
+  //Float_t rate = 1400e3;
+  Float_t rate = 1.4e3;
 
   // time in s
-  Float_t lengthRun = 1e6;
+  Float_t lengthRun = 1e5;
 
   Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
   Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
@@ -916,11 +918,14 @@ void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
     //Float_t ratePerTrigger[] = { 60, 13.3, 13.3, 13.3 };
 
     Int_t nCuts = 3;
-    Int_t cuts[] = { 0, 126, 162 };
-    Float_t ratePerTrigger[] = { 60, 20.0, 20.0 };
+    Int_t cuts[] = { 0, 114, 148 };
+
+    //Int_t nCuts = 3;
+    //Int_t cuts[] = { 0, 126, 162 };
+    //Float_t ratePerTrigger[] = { 60, 20.0, 20.0 };
 
     // desired trigger rate in Hz
-    //Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 };
+    Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 };
 
     xSection->SetStats(kFALSE);
     xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
@@ -974,8 +979,12 @@ void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
       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;
@@ -1000,11 +1009,34 @@ void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
       } else
         proj2->DrawCopy(" P");
 
-      legend2->AddEntry(proj2, Form("%lld evts, FO > %d chips", nTrigger, cut));
+      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);
 
-      TLine* line = new TLine(triggerLimit, proj2->GetMinimum(), triggerLimit, proj2->GetMaximum());
-      line->SetLineColor(colors[currentCut]);
-      line->Draw();
+      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();
@@ -1017,8 +1049,6 @@ void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
 void AliHighMultiplicitySelector::Contamination()
 {
   //
-  // produces a spectrum created with N triggers
-  // number of triggers and thresholds for the moment fixed
   //
 
   /*
@@ -1077,8 +1107,8 @@ void AliHighMultiplicitySelector::Contamination()
       //Double_t rate = rates[3];
 
       // coll. in 100 ns window
-      //Double_t windowSize = 100e-9;
-      Double_t windowSize = 25e-9;
+      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;
@@ -1262,7 +1292,7 @@ void AliHighMultiplicitySelector::Contamination2()
 
       Printf("Raw value for 2 collisions is %e", triggerRate2);
 
-      for (Double_t doubleRate = 0; doubleRate <= 0.2; doubleRate += 0.005)
+      for (Double_t doubleRate = 0; doubleRate <= 0.3; doubleRate += 0.005)
       {
         Float_t totalContamination = (triggerRate2 * doubleRate) / (triggerRate + triggerRate2 * doubleRate);
 
@@ -1277,6 +1307,119 @@ void AliHighMultiplicitySelector::Contamination2()
   }
 }
 
+void AliHighMultiplicitySelector::Contamination3()
+{
+  //
+  //
+
+  /*
+
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libPWG0base");
+  .L AliHighMultiplicitySelector.cxx+g
+  x = new AliHighMultiplicitySelector();
+  x->ReadHistograms("highmult_hijing100k.root");
+  x->Contamination3();
+
+  */
+
+  // 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"));
+
+  // prob for a collision in a bunch crossing
+  Int_t nRates = 3;
+  Float_t rates[] = {0.07, 0.1, 0.2};
+
+  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 currentRate = 0; currentRate<nRates; ++currentRate)
+    {
+      TGraph* graph = new TGraph;
+      graph->SetMarkerColor(colors[currentRate]);
+      graph->SetMarkerStyle(markers[currentRate]);
+
+      Float_t rate = rates[currentRate];
+
+      for (Int_t cut = 100; cut <= 201; cut += 10)
+      {
+        Printf("cut at %d", cut);
+
+      TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
+      Float_t* triggerEff = new Float_t[max];
+      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);
+
+      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", triggerRate3);
+
+      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);
+
+      Float_t totalContamination = (triggerRate2 * doubleRate + triggerRate3 * tripleRate) / (triggerRate * singleRate + triggerRate2 * doubleRate + triggerRate3 * tripleRate);
+
+        //Printf("Total contamination is %.1f%%", totalContamination * 100);
+
+      graph->SetPoint(graph->GetN(), cut, totalContamination);
+      }
+
+      graph->Draw((currentRate == 0) ? "A*" : "* SAME");
+      graph->GetXaxis()->SetTitle("layer 1 threshold");
+      graph->GetYaxis()->SetTitle("contamination");
+      graph->GetYaxis()->SetRangeUser(0, 1);
+    }
+  }
+}
+
 void AliHighMultiplicitySelector::DrawHistograms()
 {
   // draws the histograms
@@ -1302,6 +1445,7 @@ void AliHighMultiplicitySelector::DrawHistograms()
   x->ReadHistograms("highmult_central.root");
   x->DrawHistograms();
 
+  gSystem->Load("libANALYSIS");
   gSystem->Load("libPWG0base");
   .L AliHighMultiplicitySelector.cxx+
   x = new AliHighMultiplicitySelector();
@@ -1336,6 +1480,16 @@ void AliHighMultiplicitySelector::DrawHistograms()
   canvas->SaveAs("L1NoCurve.gif");
   canvas->SaveAs("L1NoCurve.eps");
 
+  TLine* line = new TLine(fMvsL1->GetXaxis()->GetXmin(), 150, fMvsL1->GetXaxis()->GetXmax(), 150);
+  line->SetLineWidth(2);
+  line->SetLineColor(kRed);
+  line->Draw();
+
+  canvas->SaveAs("L1NoCurveCut.gif");
+  canvas->SaveAs("L1NoCurveCut.eps");
+
+  return;
+
   // draw corresponding theoretical curve
   TF1* func = new TF1("func", "[0]*(1-(1-1/[0])**x)", 1, 1000);
   func->SetParameter(0, 400-5*2);