From: ekryshen Date: Thu, 30 Oct 2014 14:53:01 +0000 (+0100) Subject: Physics selection QA added in the QA framework X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=c0d2a42f53dbc056d8ba5f2973a90ad7764b32b1 Physics selection QA added in the QA framework --- diff --git a/PWGPP/EVS/periodLevelQA.C b/PWGPP/EVS/periodLevelQA.C new file mode 100644 index 00000000000..73f21374133 --- /dev/null +++ b/PWGPP/EVS/periodLevelQA.C @@ -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 &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 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;rGetEntry(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;ibitGetEntry(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;ibitWrite(); + 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 &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 index 00000000000..c57291574c4 --- /dev/null +++ b/PWGPP/EVS/runLevelEventStatQA.C @@ -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;itokenGetEntries();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;ibitGetBinContent(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 index 00000000000..d1a8ffb4cf7 --- /dev/null +++ b/PWGPP/EVS/triggerInfo.C @@ -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 index 00000000000..c76d3e4073c --- /dev/null +++ b/PWGPP/QA/detectorQAscripts/EVS.sh @@ -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\")" +}