Physics selection QA added in the QA framework
authorekryshen <evgeny.kryshen@cern.ch>
Thu, 30 Oct 2014 14:53:01 +0000 (15:53 +0100)
committerekryshen <evgeny.kryshen@cern.ch>
Thu, 30 Oct 2014 14:53:01 +0000 (15:53 +0100)
PWGPP/EVS/periodLevelQA.C [new file with mode: 0644]
PWGPP/EVS/runLevelEventStatQA.C [new file with mode: 0644]
PWGPP/EVS/triggerInfo.C [new file with mode: 0644]
PWGPP/QA/detectorQAscripts/EVS.sh [new file with mode: 0755]

diff --git a/PWGPP/EVS/periodLevelQA.C b/PWGPP/EVS/periodLevelQA.C
new file mode 100644 (file)
index 0000000..73f2137
--- /dev/null
@@ -0,0 +1,242 @@
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1D.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "map"
+using namespace std;
+// TODO read number of bits from AliVEvent?
+#define NBITS 29
+TString bitNames[NBITS] = {
+"kMB",
+"kINT7",
+"kMUON",
+"kHighMult",
+"kEMC1",
+"kCINT5",
+"kCMUS5",
+"kMUSH7",
+"kMUL7",
+"kMUU7",
+"kEMC7",
+"kMUS7",
+"kPHI1",
+"kPHI7/kPHI8",
+"kEMCEJE",
+"kEMCEGA",
+"kCentral",
+"kSemiCentral",
+"kDG5",
+"kZED",
+"kSPI7/kSPI8",
+"kINT8",
+"kMuonSingleLowPt8",
+"kMuonSingleHighPt8",
+"kMuonLikeLowPt8",
+"kMuonUnlikeLowPt8",
+"kMuonUnlikeLowPt0",
+"kUserDefined",
+"kTRD"
+};
+
+void SetHisto(TH1D* h);
+void AddFillSeparationLines(TH1D* h, map<Int_t,Int_t> &fills);
+
+void periodLevelQA(TString inputFileName ="trending.root"){
+  TFile* f = new TFile(inputFileName.Data());
+  TTree* t = (TTree*) f->Get("trending");
+  Int_t run;
+  Int_t fill = 0;
+  Int_t nBCsPerOrbit = 0;
+  Double_t duration = 0;
+  Double_t mu = 0;
+  Double_t lumi_seen = 0;
+  UInt_t l0b = 0;
+  Int_t all[NBITS];
+  Int_t accepted[NBITS];
+  
+  t->SetBranchAddress("run",&run);
+  t->SetBranchAddress("fill",&fill);
+  t->SetBranchAddress("bcs",&nBCsPerOrbit);
+  t->SetBranchAddress("duration",&duration);
+  t->SetBranchAddress("mu",&mu);
+  t->SetBranchAddress("l0b",&l0b);
+  t->SetBranchAddress("all",&all);
+  t->SetBranchAddress("accepted",&accepted);
+  t->SetBranchAddress("lumi_seen",&lumi_seen);
+  Int_t nRuns = t->GetEntries();
+  map<Int_t,Int_t> fills;
+  TH1D* hAll[NBITS];
+  TH1D* hAccepted[NBITS] = {0x0};
+  TH1D* hRejected[NBITS] = {0x0};
+  TH1D* hAcceptedFraction[NBITS] = {0x0};
+  TH1D* hRejectedFraction[NBITS] = {0x0};
+  
+  gStyle->SetOptStat(0);
+  gStyle->SetLineScalePS(1.5);
+  gStyle->SetPadGridX(0);
+  gStyle->SetPadGridY(1);
+  TCanvas* dummy = new TCanvas("dummy","dummy",1800,500);
+  gPad->SetMargin(0.05,0.01,0.18,0.06);
+  gPad->Print("global_properties.pdf[");
+  gPad->Print("accepted_event_statistics.pdf[");
+  gPad->Print("accepted_fraction.pdf[");
+  gPad->Print("rejected_fraction.pdf[");
+
+  TH1D* hMu       = new TH1D("hMu","Average number of minimum bias collisions per BC",nRuns,0,nRuns);
+  TH1D* hBCs      = new TH1D("hBCs","Number of colliding bunches",nRuns,0,nRuns);
+  TH1D* hDuration = new TH1D("hDuration","Duration, s",nRuns,0,nRuns);
+  TH1D* hLumiSeen = new TH1D("hLumiSeen","Luminosity seen, nb-1",nRuns,0,nRuns);
+
+  for (Int_t r=0;r<nRuns;r++){
+    t->GetEntry(r);
+    hMu->GetXaxis()->SetBinLabel(r+1,Form("%i",run));
+    hBCs->GetXaxis()->SetBinLabel(r+1,Form("%i",run));
+    hDuration->GetXaxis()->SetBinLabel(r+1,Form("%i",run));
+    hLumiSeen->GetXaxis()->SetBinLabel(r+1,Form("%i",run));
+    hMu->SetBinContent(r+1,mu);
+    hBCs->SetBinContent(r+1,nBCsPerOrbit);
+    hDuration->SetBinContent(r+1,duration);
+    hLumiSeen->SetBinContent(r+1,lumi_seen/1000000.);
+    fills[run]=fill;
+  }
+  SetHisto(hMu);
+  SetHisto(hBCs);
+  SetHisto(hDuration);
+  SetHisto(hLumiSeen);
+  TCanvas* cMu = new TCanvas("mu","mu",1800,500);
+  cMu->SetMargin(0.05,0.01,0.18,0.06);
+  hMu->Draw("h");
+  AddFillSeparationLines(hMu,fills);
+  gPad->Print("global_properties.pdf");
+
+  TCanvas* cBCs = new TCanvas("bcs","bcs",1800,500);
+  cBCs->SetMargin(0.05,0.01,0.18,0.06);
+  hBCs->Draw("h");
+  AddFillSeparationLines(hBCs,fills);
+  gPad->Print("global_properties.pdf");
+  
+  TCanvas* cDuration = new TCanvas("duration","duration",1800,500);
+  cDuration->SetMargin(0.05,0.01,0.18,0.06);
+  hDuration->SetTitle(Form("Duration in seconds: total= %.0f s = %.0f h",hDuration->Integral(),hDuration->Integral()/3600));
+  hDuration->Draw("h");
+  AddFillSeparationLines(hDuration,fills);
+  gPad->Print("global_properties.pdf");
+
+  TCanvas* cLumiSeen = new TCanvas("lumiseen","lumi seen",1800,500);
+  cLumiSeen->SetMargin(0.05,0.01,0.18,0.06);
+  hLumiSeen->SetTitle(Form("Luminosity seen [1/nb]: total= %.0f",hLumiSeen->Integral()));
+  hLumiSeen->Draw("h");
+  AddFillSeparationLines(hLumiSeen,fills);
+  gPad->Print("global_properties.pdf");
+
+
+  dummy->Print("global_properties.pdf]");
+
+  for (Int_t ibit=0;ibit<NBITS;ibit++) {
+    printf("bit=%i\n",ibit);
+    const char* bitName = bitNames[ibit];
+    hAll[ibit]      = new TH1D(Form("hAll%02i"     ,ibit),Form("All: %s"     ,bitName),nRuns,0,nRuns);
+    hAccepted[ibit] = new TH1D(Form("hAccepted%02i",ibit),Form("Accepted: %s",bitName),nRuns,0,nRuns);
+    hRejected[ibit] = new TH1D(Form("hRejected%02i",ibit),Form("Rejected: %s",bitName),nRuns,0,nRuns);
+    for (Int_t r=0;r<nRuns;r++){
+      t->GetEntry(r);
+      hAll[ibit]     ->GetXaxis()->SetBinLabel(r+1,Form("%i",run));
+      hAccepted[ibit]->GetXaxis()->SetBinLabel(r+1,Form("%i",run));
+      hRejected[ibit]->GetXaxis()->SetBinLabel(r+1,Form("%i",run));
+      hAll[ibit]->SetBinContent(r+1,all[ibit]);
+      if (all[ibit]) hAccepted[ibit]->SetBinContent(r+1,accepted[ibit]);
+      if (all[ibit]) hRejected[ibit]->SetBinContent(r+1,all[ibit]-accepted[ibit]);
+    }
+    if (hAll[ibit]->Integral()<1) continue;
+    SetHisto(hAll[ibit]);
+    SetHisto(hAccepted[ibit]);
+    SetHisto(hRejected[ibit]);
+    hAll[ibit]->Sumw2();
+    hAccepted[ibit]->Sumw2();
+    hRejected[ibit]->Sumw2();
+    hAll[ibit]->SetLineColor(kBlue);
+    hAll[ibit]->SetFillColor(kBlue);
+    hAccepted[ibit]->SetLineColor(kGreen+2);
+    hAccepted[ibit]->SetFillColor(kGreen+2);
+    hAcceptedFraction[ibit] = (TH1D*) hAll[ibit]->Clone(Form("hAcceptedFraction%02i",ibit));
+    hAcceptedFraction[ibit]->SetTitle(Form("Accepted fraction: %s",bitName));
+    hAcceptedFraction[ibit]->Divide(hAccepted[ibit],hAll[ibit],1,1,"B");
+    hAcceptedFraction[ibit]->SetFillColor(0);
+    hAcceptedFraction[ibit]->SetLineWidth(2);
+    hRejectedFraction[ibit] = (TH1D*) hAll[ibit]->Clone(Form("hRejectedFraction%02i",ibit));
+    hRejectedFraction[ibit]->SetTitle(Form("Rejected fraction: %s",bitName));
+    hRejectedFraction[ibit]->Divide(hRejected[ibit],hAll[ibit],1,1,"B");
+    hRejectedFraction[ibit]->SetFillColor(0);
+    hRejectedFraction[ibit]->SetLineWidth(2);
+    
+    TCanvas* cAll = new TCanvas(Form("all_%s",bitName),Form("All: %s",bitName),1800,500);
+    cAll->SetMargin(0.05,0.01,0.18,0.06);
+    hAll[ibit]->SetTitle(Form("%s: total=%.0f accepted=%.0f",bitName,hAll[ibit]->Integral(),hAccepted[ibit]->Integral()));
+    hAll[ibit]->Draw("h");
+    hAccepted[ibit]->Draw("h same");
+    AddFillSeparationLines(hAccepted[ibit],fills);
+    gPad->RedrawAxis();
+    gPad->Print("accepted_event_statistics.pdf");
+    
+    TCanvas* cAcceptedFraction = new TCanvas(Form("accepted_fraction_%s",bitName),Form("Accepted Fraction: %s",bitName),1800,500);
+    cAcceptedFraction->SetMargin(0.05,0.01,0.18,0.06);
+    hAcceptedFraction[ibit]->SetTitle(Form("%s: average accepted fraction = %.3f",bitName,hAccepted[ibit]->Integral()/hAll[ibit]->Integral()));
+    hAcceptedFraction[ibit]->Draw();
+    AddFillSeparationLines(hAcceptedFraction[ibit],fills);
+    gPad->Print("accepted_fraction.pdf");
+    
+    TCanvas* cRejectedFraction = new TCanvas(Form("rejected_fraction_%s",bitName),Form("Rejected Fraction: %s",bitName),1800,500);
+    cRejectedFraction->SetMargin(0.05,0.01,0.18,0.06);
+    hRejectedFraction[ibit]->SetTitle(Form("%s: average rejected fraction = %.3f",bitName,hRejected[ibit]->Integral()/hAll[ibit]->Integral()));
+    hRejectedFraction[ibit]->Draw();
+    AddFillSeparationLines(hRejectedFraction[ibit],fills);
+    gPad->Print("rejected_fraction.pdf");
+
+  }
+  dummy->Print("accepted_event_statistics.pdf]");
+  dummy->Print("accepted_fraction.pdf]");
+  dummy->Print("rejected_fraction.pdf]");
+
+  TFile* fout = new TFile("hist.root","recreate");
+  for (Int_t ibit=0;ibit<NBITS;ibit++) {
+    if (hAll[ibit]) hAll[ibit]->Write();
+    if (hAccepted[ibit]) hAccepted[ibit]->Write();
+    if (hAcceptedFraction[ibit]) hAcceptedFraction[ibit]->Write();
+    if (hRejectedFraction[ibit]) hRejectedFraction[ibit]->Write();
+  }
+  fout->Close();
+}
+
+void SetHisto(TH1D* h){
+  h->SetTitleFont(43);
+  h->SetTitleSize(25);
+  h->GetYaxis()->SetTitleFont(43);
+  h->GetXaxis()->SetLabelFont(43);
+  h->GetYaxis()->SetLabelFont(43);
+  h->GetYaxis()->SetTitleSize(25);
+  h->GetXaxis()->SetLabelSize(15);
+  h->GetYaxis()->SetLabelSize(25);
+  h->GetYaxis()->SetTickLength(0.01);
+  h->GetYaxis()->SetTitleOffset(0.6);
+  h->GetYaxis()->SetDecimals(1);
+  h->LabelsOption("av");
+  h->SetLineWidth(2);
+  h->SetLineColor(kBlue);
+  h->SetMinimum(0);
+}
+
+void AddFillSeparationLines(TH1D* h, map<Int_t,Int_t> &fills){
+  gPad->Update();
+  Double_t ymin = gPad->GetUymin();
+  Double_t ymax = gPad->GetUymax();
+  TLine * fillSeparationLine = new TLine();
+  fillSeparationLine->SetLineColor(kRed);
+  fillSeparationLine->SetLineWidth(1);
+  for(Int_t iBin = 1; iBin < h->GetXaxis()->GetNbins(); iBin++) {
+    UInt_t runOld = atoi(h->GetXaxis()->GetBinLabel(iBin));
+    UInt_t runNew = atoi(h->GetXaxis()->GetBinLabel(iBin + 1));
+    if (fills[runOld]==fills[runNew]) continue;
+    fillSeparationLine->DrawLine(iBin,ymin,iBin,ymax);
+  }
+}
diff --git a/PWGPP/EVS/runLevelEventStatQA.C b/PWGPP/EVS/runLevelEventStatQA.C
new file mode 100644 (file)
index 0000000..c572915
--- /dev/null
@@ -0,0 +1,84 @@
+#ifndef __CINT__
+#include "TTree.h"
+#include "TFile.h"
+#include "TMath.h"
+#endif
+#include "triggerInfo.C"
+
+// TODO read number of bits from AliVEvent?
+#define NBITS 29
+
+void runLevelEventStatQA(TString qafilename="/data/alice/2010/LHC10b/pass4/000119061/event_stat.root", Int_t run=168108){
+  printf("runLevelEventStatQA %s %i\n",qafilename.Data(),run);
+  TFile* fin = new TFile(qafilename);
+  
+  TH2D* h = (TH2D*) fin->Get("fHistStatistics");
+  Int_t all[NBITS]      = {0};
+  Int_t accepted[NBITS] = {0};
+  Double_t par[5] = {0};
+
+  TString refClass="";
+  Int_t refSigma;
+  if      (              run<=126437) { refSigma=  62; refClass = "CINT1B-ABCE-NOPF-ALL"; }
+  else if (run>126437 && run<=127718) { refSigma=  62; refClass = "CINT1-B-NOPF-ALLNOTRD";}
+  else if (run>127718 && run<=127730) { refSigma=  62; refClass = "CINT1B-ABCE-NOPF-ALL"; }
+  else if (run>127813 && run<=166475) { refSigma=  62; refClass = "CINT1-B-NOPF-ALLNOTRD";}
+  else if (run>136848 && run<=139517) { refSigma=7640; refClass = "CMBACS2-B-NOPF-ALLNOTRD"; }
+  else if (run>166476 && run<=170593) { refSigma=4100; refClass = "CVLN-B-NOPF-ALLNOTRD"; }
+  else if (run>195144 && run<=197388) { refSigma=1590; refClass = "C0TVX-B-NOPF-ALLNOTRD"; }
+
+  triggerInfo(run,refClass,par);
+  Int_t fill         = TMath::Nint(par[0]);
+  Double_t duration  = par[1];
+  UInt_t l0b         = TMath::Nint(par[2]);
+  Int_t nBCsPerOrbit = TMath::Nint(par[3]);
+  Double_t mu        = par[4];
+  Double_t lumi_seen = l0b/refSigma;
+  return;
+
+  for (Int_t j=1;j<=h->GetNbinsY();j++){
+    TString label = h->GetYaxis()->GetBinLabel(j);
+
+    // skip background triggers
+    // TODO introduce identifier to filter-out background triggers
+    if (!label.Contains("-B-") && !label.Contains("-S-") && !(label.Contains("-ABCE-") && label.Contains("1B-"))) continue;
+
+    // Read mask
+    // TODO think how to propagate mask with TBit aliases
+    UInt_t mask = 0;
+    TObjArray* array = label.Tokenize(" ");
+    for (Int_t itoken=0;itoken<array->GetEntries();itoken++){
+      TString token = array->At(itoken)->GetName();
+      if (token[0]!='&') continue;
+      token.Remove(0,1);
+      mask = token.Atoi();
+      break;
+    }
+    array->Delete();
+    delete array;
+    if (!mask) continue;
+    // Fill all and accepted counters for active bits
+    // TODO can we accidentally double count events?
+    for (Int_t ibit=0;ibit<NBITS;ibit++) {
+      if (!(mask & 1<<ibit)) continue; // to be changed with TBits
+      all[ibit]      += Int_t(h->GetBinContent(1             ,j)); 
+      accepted[ibit] += Int_t(h->GetBinContent(h->GetNbinsX(),j));
+    }
+  }
+  
+  TTree* t = new TTree("trending","tree of trending variables");
+  t->Branch("run",&run);
+  t->Branch("fill",&fill);
+  t->Branch("bcs",&nBCsPerOrbit);
+  t->Branch("duration",&duration);
+  t->Branch("mu",&mu);
+  t->Branch("l0b",&l0b);
+  t->Branch("lumi_seen",&lumi_seen);
+  t->Branch("all",&all,Form("all[%i]/I",NBITS));
+  t->Branch("accepted",&accepted,Form("accepted[%i]/I",NBITS));
+  t->Fill();
+  
+  TFile* fout = new TFile(Form("trending_%i.root",run),"recreate");
+  t->Write();
+  fout->Close();
+}
diff --git a/PWGPP/EVS/triggerInfo.C b/PWGPP/EVS/triggerInfo.C
new file mode 100644 (file)
index 0000000..d1a8ffb
--- /dev/null
@@ -0,0 +1,89 @@
+#ifndef __CINT__
+#include "TMath.h"
+#include "AliCDBManager.h"
+#include "AliTriggerScalers.h"
+#include "AliTriggerRunScalers.h"
+#include "AliTimeStamp.h"
+#include "AliTriggerScalersRecord.h"
+#include "AliTriggerConfiguration.h"
+#include "AliLHCData.h"
+#include "AliTriggerClass.h"
+#include "AliTriggerBCMask.h"
+#include "AliCDBPath.h"
+#include "AliCDBEntry.h"
+#endif
+
+UInt_t dif(UInt_t stop, UInt_t start){
+  UInt_t d;
+  if(stop >= start) d=stop-start;
+  else d = stop+(0xffffffff-start)+1;
+  return d;
+};
+
+Int_t triggerInfo(Int_t run, TString refClassName, Double_t* par){
+  AliCDBManager* man = AliCDBManager::Instance();
+  //  man->SetDefaultStorage("local:///data/alice/OCDB");
+  man->SetRun(run);
+  AliLHCData* lhc = (AliLHCData*) man->Get("GRP/GRP/LHCData")->GetObject();
+  par[0] = lhc->GetFillNumber();
+  
+  AliTriggerConfiguration* cfg = (AliTriggerConfiguration*) man->Get("GRP/CTP/Config")->GetObject();
+
+  // Get scalers 
+  AliTriggerRunScalers* scalers = (AliTriggerRunScalers*) man->Get("GRP/CTP/Scalers")->GetObject();
+  Int_t nEntries = scalers->GetScalersRecords()->GetEntriesFast();
+  
+  // Get SOR and EOR scaler records
+  AliTriggerScalersRecord* record1 = scalers->GetScalersRecord(0);
+  AliTriggerScalersRecord* record2 = scalers->GetScalersRecord(nEntries-1);
+  if (!record1) { printf("Null pointer to start scalers record\n"); return 1; }
+  if (!record2) { printf("Null pointer to  stop scalers record\n"); return 1; }
+
+  // Extract SOR and EOR times
+  const AliTimeStamp* stemp1 = record1->GetTimeStamp();
+  const AliTimeStamp* stemp2 = record2->GetTimeStamp();
+  if (!stemp1) { printf("Null pointer to start timestemp\n"); return 1; }
+  if (!stemp2) { printf("Null pointer to stop timestemp\n");  return 1; }
+  Double_t duration = stemp2->GetSeconds()-stemp1->GetSeconds();
+  par[1] = duration;
+  if (TMath::Abs(duration)<1) return 2;
+
+  // Extract SOR and EOR trigger counts
+  Int_t classid = cfg->GetClassIndexFromName(refClassName);
+  const AliTriggerScalers* scaler1 = record1->GetTriggerScalersForClass(classid);
+  const AliTriggerScalers* scaler2 = record2->GetTriggerScalersForClass(classid);
+  if (!scaler1) { printf("Null pointer to start scalers for reference class\n"); return 1; }
+  if (!scaler2) { printf("Null pointer to stop scalers for reference class\n");  return 1; }
+  UInt_t l0b = dif(scaler2->GetLOCB(),scaler1->GetLOCB());
+  par[2] = l0b;
+
+  // Get number of colliding bunches per orbit
+  Double_t orbitRate = 11245.; // Hz
+  Double_t nBCsPerOrbit = -1;
+  if (refClassName.Contains("1B-ABCE-")){
+    nBCsPerOrbit = lhc->GetNInteractingBunchesMeasured();
+    Printf("Number of BCs from LHC data=%i",nBCsPerOrbit);
+    if (nBCsPerOrbit<0) {
+      Int_t emptyclassid = cfg->GetClassIndexFromName("CBEAMB-ABCE-NOPF-ALL");
+      if (emptyclassid<0) return 3;
+      const AliTriggerScalers* emptyScaler1 = record1->GetTriggerScalersForClass(emptyclassid);
+      const AliTriggerScalers* emptyScaler2 = record2->GetTriggerScalersForClass(emptyclassid);
+      if (!scaler1) { printf("Null pointer to start scalers for reference class\n"); return 1; }
+      if (!scaler2) { printf("Null pointer to stop scalers for reference class\n");  return 1; }
+      UInt_t l0bempty = dif(emptyScaler2->GetLOCB(),emptyScaler1->GetLOCB());
+      if (l0bempty==0) return 4;
+      nBCsPerOrbit = l0bempty/orbitRate/duration;
+    }
+  }
+  else {
+    // Extract number of bunches per orbit
+    AliTriggerClass* cl = cfg->GetTriggerClass(classid);
+    AliTriggerBCMask* mask = cl->GetBCMask();
+    nBCsPerOrbit = mask->GetNUnmaskedBCs();
+  }
+  par[3] = nBCsPerOrbit;
+  
+  Double_t totalBCs = duration*orbitRate*nBCsPerOrbit;
+  par[4] = -TMath::Log(1-l0b/totalBCs); // mu
+  return 0;
+}
diff --git a/PWGPP/QA/detectorQAscripts/EVS.sh b/PWGPP/QA/detectorQAscripts/EVS.sh
new file mode 100755 (executable)
index 0000000..c76d3e4
--- /dev/null
@@ -0,0 +1,15 @@
+runLevelEventStatQA()
+{
+  eventStatFile=$1
+
+  cp $ALICE_ROOT/PWGPP/EVS/runLevelEventStatQA.C .
+  aliroot -b -q -l "runLevelEventStatQA.C(\"$eventStatFile\",${runNumber})" 
+}
+
+periodLevelQA()
+{
+  trendingFile=$1
+
+  cp $ALICE_ROOT/PWGPP/EVS/periodLevelQA.C .
+  aliroot -b -q -l "periodLevelQA.C(\"trending.root\")"
+}