diff --git a/QC/timeseries/DoMerge.C b/QC/timeseries/DoMerge.C new file mode 100644 index 0000000..8506aa3 --- /dev/null +++ b/QC/timeseries/DoMerge.C @@ -0,0 +1,23 @@ +#include "TFileMerger.h" + +// macro to merge all the files in a list (they should belong to the same school) +// The header of the first file will be kept +void DoMerge(char* nomelista="listaM",char* output="filtered.root"){ +// system("ls calib/*.root >listaM"); + FILE *f = fopen(nomelista,"r"); + +// TGrid::Connect("alien://"); + + TFileMerger m(kFALSE); + m.OutputFile(output); + + Int_t i=0; + char nome[100]; + while (fscanf(f,"%s",nome)==1) { + m.AddFile(nome); + i++; + } + if (i) + m.Merge(); +} + diff --git a/QC/timeseries/README b/QC/timeseries/README new file mode 100644 index 0000000..98734fe --- /dev/null +++ b/QC/timeseries/README @@ -0,0 +1,18 @@ +# creating our tree from TPC timeseries accessing inputs from remote (alien), note that in the script apass5 is assumed (in case change it) +./dorun.sh PERIOD RUN +#example ./dorun.sh LHC23zzh 544091 +# output -> 544091.root + + +# process the output (files list taken from listaT) +# requirements +# 1)dedx peramaterization (provided a default here) -> expDeDx.root +# 2) mismatch plots from qcdb -> etaMism.root, mism.root +# qc_async/TOF/MO/PID/mHistMismatchVsEta +# qc_async/TOF/MO/PID/LoverCvsEta + +ls YOUR_RUNS > listaT +root -b -q -l process.C + +# check paramerization of track contribution to expected time resolution (only antiparticles used) +root getRes.C diff --git a/QC/timeseries/checkrun.sh b/QC/timeseries/checkrun.sh new file mode 100755 index 0000000..8f7cb86 --- /dev/null +++ b/QC/timeseries/checkrun.sh @@ -0,0 +1,6 @@ +#!/bin/bash +export period=$2 +export year=$(echo $2|awk -F"LHC" '{print "20"substr($2,1,2)}') + +echo $year $period $1 +alien.py find /alice/data/$year/${period}/$1/apass5/ time_series_tracks_0.root |grep Stage_5 diff --git a/QC/timeseries/dorun.sh b/QC/timeseries/dorun.sh new file mode 100755 index 0000000..7c7d0e3 --- /dev/null +++ b/QC/timeseries/dorun.sh @@ -0,0 +1,12 @@ +#!/bin/bash +export period=$2 +export year=$(echo $2|awk -F"LHC" '{print "20"substr($2,1,2)}') + +echo alien.py find /alice/data/$year/${period}/$1/apass5/ time_series_tracks_0.root + +alien.py find /alice/data/$year/${period}/$1/apass5/ time_series_tracks_0.root |grep Stage_5 >lista +rm -rf $1 +mkdir $1 +cat -b lista|awk '{print "./filter.sh","alien://"$2,"'$1'/"$1".root"}'|bash + +ls $1/*.root >listaM ; root -b -q -l DoMerge.C\(\"listaM\",\"$1.root\"\) diff --git a/QC/timeseries/expDeDx.root b/QC/timeseries/expDeDx.root new file mode 100644 index 0000000..0304e45 Binary files /dev/null and b/QC/timeseries/expDeDx.root differ diff --git a/QC/timeseries/filter b/QC/timeseries/filter new file mode 100644 index 0000000..4da85c7 --- /dev/null +++ b/QC/timeseries/filter @@ -0,0 +1,14 @@ +double 1 mTOFSignal tof +float 9 mTOFLength.mT[9] texp +float 1 pt pt +float 1 mDXatTOF dx +float 1 mDZatTOF dz +float 1 vertex_time vertex_time +int 1 ncl ncl +UShort_t 1 vertex_nContributors vertex_nContributors +UInt_t 1 TOFmask TOFmask +float 1 phi phi +float 1 tgl tgl +float 1 dedxTPC.dEdxTotTPC dedx +UInt_t 1 firstTFOrbit firstTFOrbit +float 1 qpt qpt diff --git a/QC/timeseries/filter.C b/QC/timeseries/filter.C new file mode 100644 index 0000000..65fd9a1 --- /dev/null +++ b/QC/timeseries/filter.C @@ -0,0 +1,118 @@ +void filter(const char* fname="time_series_tracks_0.root",const char* fnameOut="filtered.root"){ + int typoV[1000]; + int intV[1000]; // 0 + float floatV[1000]; // 1 + double doubleV[1000]; // 2 + UShort_t ushortV[1000]; // 3 + UInt_t uintV[1000]; // 4 + + float *pt; + int *ncl; + + // float downscaleextra = 1; + + TGrid::Connect("alien://"); + + TFile *fts = TFile::Open(fname); + TTree *t = (TTree *) fts->Get("treeTimeSeries"); + + TFile *fout = new TFile(fnameOut,"RECREATE"); + TTree *tout = new TTree("tree","tree"); + + char hasITS; + t->SetBranchAddress("hasITSTPC",&hasITS); + + FILE *f = fopen("filter","r"); + char typeV[30],labelIn[100],labelOut[100]; + int ndim; + int np=0; + while(fscanf(f,"%s %d %s %s",typeV,&ndim,labelIn,labelOut) == 4){ + printf("new\n"); + TString a(typeV); + TString b(labelIn); + if(a.Contains("int")){ + printf("%s -> int\n",typeV); + typoV[np] = 0; + t->SetBranchAddress(labelIn,&(intV[np])); + if(ndim==1){ + tout->Branch(labelOut,&(intV[np])); + } else { + tout->Branch(labelOut,&(intV[np]),Form("%s[%d]",labelOut,ndim)); + } + + if(b.Contains("ncl")){ + ncl = &(intV[np]); + } + + } else if(a.Contains("float")){ + printf("%s -> float\n",typeV); + typoV[np] = 1; + t->SetBranchAddress(labelIn,&(floatV[np])); + if(ndim==1){ + tout->Branch(labelOut,&(floatV[np])); + } else { + tout->Branch(labelOut,&(floatV[np]),Form("%s[%d]",labelOut,ndim)); + } + + if(b.Contains("pt") && !b.Contains("qpt")){ + pt = &(floatV[np]); + } + + } else if(a.Contains("double")){ + printf("%s -> double\n",typeV); + typoV[np] = 2; + t->SetBranchAddress(labelIn,&(doubleV[np])); + if(ndim==1){ + tout->Branch(labelOut,&(doubleV[np])); + } else { + tout->Branch(labelOut,&(doubleV[np]),Form("%s[%d]",labelOut,ndim)); + } + } else if(a.Contains("UShort_t")){ + printf("%s -> ushort\n",typeV); + typoV[np] = 2; + t->SetBranchAddress(labelIn,&(ushortV[np])); + if(ndim==1){ + tout->Branch(labelOut,&(ushortV[np])); + } else { + tout->Branch(labelOut,&(ushortV[np]),Form("%s[%d]",labelOut,ndim)); + } + } else if(a.Contains("UInt_t")){ + printf("%s -> uint\n",typeV); + typoV[np] = 2; + t->SetBranchAddress(labelIn,&(uintV[np])); + if(ndim==1){ + tout->Branch(labelOut,&(uintV[np])); + } else { + tout->Branch(labelOut,&(uintV[np]),Form("%s[%d]",labelOut,ndim)); + } + } + np+=ndim; + } + + // factor 10 downscale for pT < 1 + for(int i=0; i < t->GetEntries();i++){ + t->GetEvent(i); + + /* + if(gRandom->Rndm() < downscaleextra){ + continue; + } + */ + + if(! hasITS){ + continue; + } + + if(*pt > 40 || *ncl < 100){ + continue; + } + if(1 || *pt > 1 || gRandom->Rndm() < 0.1){ + tout->Fill(); + } + } + + fout->cd(); + tout->Write(); + fout->Close(); + +} diff --git a/QC/timeseries/filter.sh b/QC/timeseries/filter.sh new file mode 100755 index 0000000..8fda5c7 --- /dev/null +++ b/QC/timeseries/filter.sh @@ -0,0 +1,4 @@ +#!/bin/bash + +root -b -q -l filter.C\(\"$1\",\"$2\"\) + diff --git a/QC/timeseries/getRes.C b/QC/timeseries/getRes.C new file mode 100644 index 0000000..acc5495 --- /dev/null +++ b/QC/timeseries/getRes.C @@ -0,0 +1,232 @@ +TF1 *func,*fp0; +float sigma=65; +TFile *f; + +Double_t fun2(Double_t *x, Double_t *par) { + double trk = par[0]*TMath::Power((TMath::Max(x[0]-0.417,0.1))*(1-par[1]*x[1]*x[1]),-par[2]); + return sqrt(trk*trk + sigma*sigma); +} + +void tuneBack(TH1D *hB, TH1D *hS){ + TH1D hr(*hB); + hr.SetName("htemp"); + hr.Divide(hS); + hr.Fit(fp0,"WW"); + hB->Scale(1./fp0->GetParameter(0)); +} + +float getResEta(float par[6],const char* species="Pi",int binEta=0, int binPt=10){ + TH2F *hS = (TH2F *) f->Get(Form("h%sEta_%d",species,binEta)); + TH2F *hB = (TH2F *) f->Get(Form("h%sEtaM_%d",species,binEta)); + + TH1D *hSig = hS->ProjectionY("hSig",binPt,binPt); + TH1D *hBack = hB->ProjectionY("hBack",binPt,binPt); + tuneBack(hBack, hSig); + + hSig->Draw(); + + hSig->Add(hBack,-1); + + hSig->Fit(func,"","",-200,150); + hSig->Fit(func,"WW","",-200,150); + + hBack->Draw("SAME"); + hBack->SetLineColor(2); + + par[0] = func->GetParameter(0); + par[1] = func->GetParameter(1); + par[2] = func->GetParameter(2); + par[3] = hS->GetXaxis()->GetBinCenter(binPt); + par[4] = func->GetParError(1); + par[5] = func->GetParError(2); + + return par[2]; +} + +void getRes(){ + /* from O2Physics/Common/Core/PID/PIDTOF.h + static constexpr std::array mDefaultResoParams{"14.3*TMath::Power((TMath::Max(x-0.319,0.1))*(1-0.4235*y*y),-0.8467)", + "14.3*TMath::Power((TMath::Max(x-0.319,0.1))*(1-0.4235*y*y),-0.8467)", + "14.3*TMath::Power((TMath::Max(x-0.319,0.1))*(1-0.4235*y*y),-0.8467)", + "42.66*TMath::Power((TMath::Max(x-0.417,0.1))*(1-0.4235*y*y),-0.7145)", + "99.46*TMath::Power((TMath::Max(x-0.447,0.1))*(1-0.4235*y*y),-0.8094)", + "216*TMath::Power((TMath::Max(x-0.647,0.1))*(1-0.4235*y*y),-0.76)", + "315*TMath::Power((TMath::Max(x-0.811,0.1))*(1-0.4235*y*y),-0.783)", + "157*TMath::Power((TMath::Max(x-0.556,0.1))*(1-0.4235*y*y),-0.783)", + "216*TMath::Power((TMath::Max(x-0.647,0.1))*(1-0.4235*y*y),-0.76)"}; + */ + + static constexpr std::array mDefaultResoParams{"14.3*TMath::Power((TMath::Max(x-0.319,0.1))*(1-1.0235*y*y),-0.8467)", + "14.3*TMath::Power((TMath::Max(x-0.319,0.1))*(1-1.0235*y*y),-0.8467)", + "23*TMath::Power((TMath::Max(x-0.319,0.2))*(1-1.3*y*y),-0.373601)", + "42.3505*TMath::Power((TMath::Max(x-0.417,0.1))*(1-1.13022*y*y),-0.588796)", + "108.89*TMath::Power((TMath::Max(x-0.447,0.1))*(1-0.573154*y*y),-0.580692)", + "236*TMath::Power((TMath::Max(x-0.647,0.1))*(1-0.4235*y*y),-0.76)", + "335*TMath::Power((TMath::Max(x-0.811,0.1))*(1-0.4235*y*y),-0.783)", + "170*TMath::Power((TMath::Max(x-0.556,0.1))*(1-0.4235*y*y),-0.783)", + "236*TMath::Power((TMath::Max(x-0.647,0.1))*(1-0.4235*y*y),-0.76)"}; + + TF2 *ffit = new TF2("ffit",fun2,0,3,-0.8,0.8,3); + + TF2 *fpi = new TF2("fpi",mDefaultResoParams[2],0,3,-0.8,0.8); + TF2 *fka = new TF2("fka",mDefaultResoParams[3],0,3,-0.8,0.8); + TF2 *fpr = new TF2("fpr",mDefaultResoParams[4],0,3,-0.8,0.8); + TF2 *fde = new TF2("fde",mDefaultResoParams[5],0,3,-0.8,0.8); + + f = new TFile("summary.root"); + + fp0 = new TF1("fp0","pol0"); + + // func = new TF1("func","[0]*TMath::Exp(-(x-[1])*(x-[1])/(2*([2]*[2]+[3]*[3])))",-150,150); + func = new TF1("func","gaus"); + func->SetParameter(0,1); + func->SetParameter(1,0); + func->SetParLimits(1,-50,50); + func->SetParameter(2,100); + func->SetParLimits(2,sigma,2000); + // func->FixParameter(3,sigma); + func->SetNpx(2000); + + float par[6]; + + int imin = 25; + int imax = 48; + + double xPt[2000]; + double xPtE[2000]; + double xEta[2000]; + double xEtaE[2000]; + double xMean[4][2000]; + double xSigma[4][2000]; + double xMeanE[4][2000]; + double xSigmaE[4][2000]; + int np = 0; + + for(int i=imax; i>= imin; i--){ + for(int j=0; j< 16; j++){ + xEta[np] = -0.75 + 0.1*j; + xEtaE[np] = 0; + getResEta(par, "Pi", j, i); + xPt[np] = -par[3]; + xPtE[np] = 0; + float piRes = sqrt(fpi->Eval(xPt[np],xEta[np]) * fpi->Eval(xPt[np],xEta[np]) + sigma*sigma); + xMean[0][np] = par[1]; + xSigma[0][np] = par[2] / piRes; + xMeanE[0][np] = par[4]; + xSigmaE[0][np] = par[5] / piRes; + + getResEta(par, "Ka", j, i); + float kaRes = sqrt(fka->Eval(xPt[np],xEta[np]) * fka->Eval(xPt[np],xEta[np]) + sigma*sigma); + xMean[1][np] = par[1]; + xSigma[1][np] = par[2] / kaRes; + xMeanE[1][np] = par[4]; + xSigmaE[1][np] = par[5] / kaRes; + + getResEta(par, "Pr", j, i); + float prRes = sqrt(fpr->Eval(xPt[np],xEta[np]) * fpr->Eval(xPt[np],xEta[np]) + sigma*sigma); + xMean[2][np] = par[1]; + xSigma[2][np] = par[2] / prRes; + xMeanE[2][np] = par[4]; + xSigmaE[2][np] = par[5] / prRes; + + /* + getResEta(par, "De", j, i); + float deRes = sqrt(fde->Eval(xPt[np],xEta[np]) * fde->Eval(xPt[np],xEta[np]) + sigma*sigma); + xMean[3][np] = par[1]; + xSigma[3][np] = par[2] / deRes; + xMeanE[3][np] = par[4]; + xSigmaE[3][np] = par[5] / deRes; + */ + np++; + } + } + + TCanvas *c = new TCanvas("cpions","cpions"); + c->Divide(2,1); + c->cd(1); + TGraph2DErrors *gPiM = new TGraph2DErrors(np,xPt,xEta,xMean[0],xPtE,xEtaE,xMeanE[0]); + gPiM->SetName("gPiM"); + gPiM->SetMarkerStyle(20); + gPiM->Draw("SURF2"); + gPiM->SetLineColor(1); + gPiM->SetMarkerColor(1); + c->cd(2); + TGraph2DErrors *gPiS = new TGraph2DErrors(np,xPt,xEta,xSigma[0],xPtE,xEtaE,xSigmaE[0]); + gPiS->SetName("gPiS"); + gPiS->SetMarkerStyle(20); + gPiS->Draw("SURF2"); + gPiS->SetLineColor(1); + gPiS->SetMarkerColor(1); + + gPiM->SetMinimum(-50); + gPiM->SetMaximum(50); + gPiS->SetMinimum(0); + gPiS->SetMaximum(2); + + c = new TCanvas("ckaons","ckaons"); + c->Divide(2,1); + c->cd(1); + TGraph2DErrors *gKaM = new TGraph2DErrors(np-16*7,xPt,xEta,xMean[1],xPtE,xEtaE,xMeanE[1]); + gKaM->SetName("gKaM"); + gKaM->SetMarkerStyle(20); + gKaM->Draw("SURF2"); + gKaM->SetLineColor(2); + gKaM->SetMarkerColor(2); + c->cd(2); + TGraph2DErrors *gKaS = new TGraph2DErrors(np- 16*7,xPt,xEta,xSigma[1],xPtE,xEtaE,xSigmaE[1]); + gKaS->SetName("gKaS"); + gKaS->SetMarkerStyle(20); + gKaS->Draw("SURF2"); + gKaS->SetLineColor(2); + gKaS->SetMarkerColor(2); + + gKaM->SetMinimum(-50); + gKaM->SetMaximum(50); + gKaS->SetMinimum(0); + gKaS->SetMaximum(2); + + c = new TCanvas("cprotons","cprotons"); + c->Divide(2,1); + c->cd(1); + TGraph2DErrors *gPrM = new TGraph2DErrors(np,xPt,xEta,xMean[2],xPtE,xEtaE,xMeanE[2]); + gPrM->SetName("gPrM"); + gPrM->SetMarkerStyle(20); + gPrM->Draw("SURF2"); + gPrM->SetLineColor(4); + gPrM->SetMarkerColor(4); + c->cd(2); + TGraph2DErrors *gPrS = new TGraph2DErrors(np,xPt,xEta,xSigma[2],xPtE,xEtaE,xSigmaE[2]); + gPrS->SetName("gPrS"); + gPrS->SetMarkerStyle(20); + gPrS->Draw("SURF2"); + gPrS->SetLineColor(4); + gPrS->SetMarkerColor(4); + + gPrM->SetMinimum(-50); + gPrM->SetMaximum(50); + gPrS->SetMinimum(0); + gPrS->SetMaximum(2); + + return; + c = new TCanvas("cdeuterons","cdeuterons"); + c->Divide(2,1); + c->cd(1); + TGraph2DErrors *gDeM = new TGraph2DErrors(np,xPt,xEta,xMean[3],xPtE,xEtaE,xMeanE[3]); + gDeM->SetName("gDeM"); + gDeM->SetMarkerStyle(20); + gDeM->Draw("SURF2"); + gDeM->SetLineColor(6); + gDeM->SetMarkerColor(6); + c->cd(2); + TGraph2DErrors *gDeS = new TGraph2DErrors(np,xPt,xEta,xSigma[3],xPtE,xEtaE,xSigmaE[3]); + gDeS->SetName("gDeS"); + gDeS->SetMarkerStyle(20); + gDeS->Draw("SURF2"); + gDeS->SetLineColor(6); + gDeS->SetMarkerColor(6); + + gDeM->SetMinimum(-50); + gDeM->SetMaximum(50); + gDeS->SetMinimum(0); + gDeS->SetMaximum(2); +} diff --git a/QC/timeseries/process.C b/QC/timeseries/process.C new file mode 100644 index 0000000..69c5ac2 --- /dev/null +++ b/QC/timeseries/process.C @@ -0,0 +1,341 @@ +void process(){ + float tpcthr = 100000000; + float minBetaGamma = 0.4; + + TChain *t = new TChain("tree"); + FILE *f = fopen("listaT","r"); + char fname[100]; + while(fscanf(f,"%s",fname) == 1){ + t->AddFile(fname); + } + TFile *feloss = TFile::Open("expDeDx.root"); + TH1D *eloss = nullptr; + if(feloss){ + eloss = (TH1D *) feloss->Get("hDeDxExp_px"); + eloss->SetName("eloss"); + int iminBetaGamma = eloss->FindBin(minBetaGamma); + for(int i=1; i <= eloss->GetNbinsX(); i++){ + if(i < iminBetaGamma) { + eloss->SetBinContent(i,eloss->GetBinContent(iminBetaGamma)); + } + } + } + + TFile *fEtamism = new TFile("etaMism.root"); + TProfile *hEtaMism = (TProfile *) fEtamism->Get("ccdb_object"); + TFile *fMism = new TFile("mism.root"); + TH2F *hMism = (TH2F *) fMism->Get("ccdb_object"); + + const int netaBins = 21; + TH1D *hMismPro[netaBins]; + for(int i=0; i < netaBins; i++){ + hMismPro[i] = hMism->ProjectionY(Form("mism_pro%d",i),i,i); + } + + double tof; t->SetBranchAddress("tof",&tof); + float texp[9]; t->SetBranchAddress("texp",texp); + float pt; t->SetBranchAddress("pt",&pt); + float dx; t->SetBranchAddress("dx",&dx); + float dz; t->SetBranchAddress("dz",&dz); + float vtime; t->SetBranchAddress("vertex_time",&vtime); + int ncl; t->SetBranchAddress("ncl",&ncl); + UShort_t mult; t->SetBranchAddress("vertex_nContributors",&mult); + UInt_t mask; t->SetBranchAddress("TOFmask",&mask); + float tgl; t->SetBranchAddress("tgl",&tgl); + float phi; t->SetBranchAddress("phi",&phi); + float dedx; t->SetBranchAddress("dedx",&dedx); + float qpt; t->SetBranchAddress("qpt",&qpt); + + double bctime = o2::tof::Geo::BC_TIME_INPS; + double bctimeInv = o2::tof::Geo::BC_TIME_INPS_INV; + + long nev = t->GetEntries(); + + const int nbinEta = 16; + float etaBin = 0.1; + float minEta = -etaBin*nbinEta/2; + float maxEta = etaBin*nbinEta/2; + + TH2F *hPiEta[nbinEta]; + TH2F *hKaEta[nbinEta]; + TH2F *hPrEta[nbinEta]; + TH2F *hDeEta[nbinEta]; + TH2F *hPiEtaM[nbinEta]; + TH2F *hKaEtaM[nbinEta]; + TH2F *hPrEtaM[nbinEta]; + TH2F *hDeEtaM[nbinEta]; + for(int i=0; i < nbinEta; i++){ + hPiEta[i] = new TH2F(Form("hPiEta_%d",i),Form("%.1f < #eta < %.1f;p (GeV/c); t_{TOF} - t_{exp}^{#pi} - t_{0}^{FT0-AC} (ps)",minEta,maxEta),100,-5,5,400,-5000,5000); + hKaEta[i] = new TH2F(Form("hKaEta_%d",i),Form("%.1f < #eta < %.1f;p (GeV/c); t_{TOF} - t_{exp}^{#pi} - t_{0}^{FT0-AC} (ps)",minEta,maxEta),100,-5,5,400,-5000,5000); + hPrEta[i] = new TH2F(Form("hPrEta_%d",i),Form("%.1f < #eta < %.1f;p (GeV/c); t_{TOF} - t_{exp}^{#pi} - t_{0}^{FT0-AC} (ps)",minEta,maxEta),100,-5,5,400,-5000,5000); + hDeEta[i] = new TH2F(Form("hDeEta_%d",i),Form("%.1f < #eta < %.1f;p (GeV/c); t_{TOF} - t_{exp}^{#pi} - t_{0}^{FT0-AC} (ps)",minEta,maxEta),100,-5,5,400,-5000,5000); + hPiEtaM[i] = new TH2F(Form("hPiEtaM_%d",i),Form("%.1f < #eta < %.1f;p (GeV/c); t_{TOF} - t_{exp}^{#pi} - t_{0}^{FT0-AC} (ps)",minEta,maxEta),100,-5,5,400,-5000,5000); + hKaEtaM[i] = new TH2F(Form("hKaEtaM_%d",i),Form("%.1f < #eta < %.1f;p (GeV/c); t_{TOF} - t_{exp}^{#pi} - t_{0}^{FT0-AC} (ps)",minEta,maxEta),100,-5,5,400,-5000,5000); + hPrEtaM[i] = new TH2F(Form("hPrEtaM_%d",i),Form("%.1f < #eta < %.1f;p (GeV/c); t_{TOF} - t_{exp}^{#pi} - t_{0}^{FT0-AC} (ps)",minEta,maxEta),100,-5,5,400,-5000,5000); + hDeEtaM[i] = new TH2F(Form("hDeEtaM_%d",i),Form("%.1f < #eta < %.1f;p (GeV/c); t_{TOF} - t_{exp}^{#pi} - t_{0}^{FT0-AC} (ps)",minEta,maxEta),100,-5,5,400,-5000,5000); + } + + TProfile *hEff = new TProfile("hEff","|charge|=1;p_{T} (GeV/c);#varepsilon",1000,-5,5); + TProfile *hEff2 = new TProfile("hEff2","|charge|=2;p_{T} (GeV/c);#varepsilon",1000,-5,5); + TH1F *hQPt = new TH1F("hQPt",";Q/p_{T}",100,-10,10); + TH2F *hPi = new TH2F("hPi",";p_{T} (GeV/c); t_{TOF} - t_{exp}^{#pi} - t_{0}^{FT0-AC} (ps)",100,-5,5,300,-2000,10000); + TH2F *hKa = new TH2F("hKa",";p_{T} (GeV/c); t_{TOF} - t_{exp}^{K} - t_{0}^{FT0-AC} (ps)",100,-5,5,300,-2000,10000); + TH2F *hPr = new TH2F("hPr",";p_{T} (GeV/c); t_{TOF} - t_{exp}^{p} - t_{0}^{FT0-AC} (ps)",100,-5,5,300,-2000,10000); + TH2F *hDe = new TH2F("hDe",";p_{T} (GeV/c); t_{TOF} - t_{exp}^{D} - t_{0}^{FT0-AC} (ps)",100,-5,5,300,-2000,10000); + TH2F *hTr = new TH2F("hTr",";p_{T} (GeV/c); t_{TOF} - t_{exp}^{T} - t_{0}^{FT0-AC} (ps)",100,-5,5,300,-10000,10000); + TH2F *hHe = new TH2F("hHe",";p_{T} (GeV/c); t_{TOF} - t_{exp}^{^{3}He} - t_{0}^{FT0-AC} (ps)",100,-5,5,300,-2000,10000); + TH2F *hAl = new TH2F("hAl",";p_{T} (GeV/c); t_{TOF} - t_{exp}^{^{4}He} - t_{0}^{FT0-AC} (ps)",100,-5,5,300,-2000,10000); + TH2F *hPiTPC = new TH2F("hPiTPC",";p_{T} (GeV/c); dE/dx - dE/dx^{exp}",100,-5,5,300,-100,100); + TH2F *hKaTPC = new TH2F("hKaTPC",";p_{T} (GeV/c); dE/dx - dE/dx^{exp}",100,-5,5,300,-100,100); + TH2F *hPrTPC = new TH2F("hPrTPC",";p_{T} (GeV/c); dE/dx - dE/dx^{exp}",100,-5,5,300,-100,100); + TH2F *hDeTPC = new TH2F("hDeTPC",";p_{T} (GeV/c); dE/dx - dE/dx^{exp}",100,-5,5,300,-100,100); + TH2F *hTrTPC = new TH2F("hTrTPC",";p_{T} (GeV/c); dE/dx - dE/dx^{exp}",100,-5,5,300,-100,100); + TH2F *hHeTPC = new TH2F("hHeTPC",";p_{T} (GeV/c); dE/dx - dE/dx^{exp}",100,-10,10,300,-100,100); + TH2F *hAlTPC = new TH2F("hAlTPC",";p_{T} (GeV/c); dE/dx - dE/dx^{exp}",100,-10,10,300,-100,100); + TH2F *hMass = new TH2F("hMass","|charge|=1;Mass (GeV/c^{2});p_{T}",2000,-4,4,50,0,5); + TH2F *hMass2 = new TH2F("hMass2","|charge|=2;Mass (GeV/c^{2});p_{T}",3000,-6,6,50,0,10); + TH2F *hPiM = new TH2F("hPiM",";p_{T} (GeV/c); t_{mism} - t_{exp}^{#pi} - t_{0}^{FT0-AC} (ps)",100,-5,5,300,-2000,10000); + TH2F *hKaM = new TH2F("hKaM",";p_{T} (GeV/c); t_{mism} - t_{exp}^{K} - t_{0}^{FT0-AC} (ps)",100,-5,5,300,-2000,10000); + TH2F *hPrM = new TH2F("hPrM",";p_{T} (GeV/c); t_{mism} - t_{exp}^{p} - t_{0}^{FT0-AC} (ps)",100,-5,5,300,-2000,10000); + TH2F *hDeM = new TH2F("hDeM",";p_{T} (GeV/c); t_{mism} - t_{exp}^{D} - t_{0}^{FT0-AC} (ps)",100,-5,5,300,-2000,10000); + TH2F *hTrM = new TH2F("hTrM",";p_{T} (GeV/c); t_{mism} - t_{exp}^{T} - t_{0}^{FT0-AC} (ps)",100,-5,5,300,-10000,10000); + TH2F *hMassM = new TH2F("hMassM","|charge|=1;Mass (GeV/c^{2});p_{T}",2000,-4,4,50,0,5); + TH2F *hDeDx = new TH2F("hDeDx","|charge|=1;#beta#gamma = p/(M_{TOF} c);dE/dx_{TPC}",400,-10,10,1600,0,800); + TProfile *hDeDxExp = new TProfile("hDeDxExp","|charge|=1;#beta#gamma = p/(M_{TOF} c);dE/dx_{TPC}",400,0,20); + TH2F *hDeDxTh = new TH2F("hDeDxTh","|charge|=1;#beta#gamma = p/(M_{TOF} c);dE/dx_{TPC} - dE/dx_{expected} ",400,-10,10,1000,-200,200); + TH2F *hDeDxThM = new TH2F("hDeDxThM","|charge|=1, p < 2 GeV/c;M_{TOF} (GeV/c^{2});dE/dx_{TPC} - dE/dx_{expected}",1000,-4,4,500,-200,200); + TH2F *hDeDxThMr = new TH2F("hDeDxThMr",";M_{TOF} (GeV/c^{2});dE/dx_{TPC} / dE/dx_{expected}",1000,-5,5,500,0,10); + + for(long i=0; i < nev; i++){ + t->GetEvent(i); + hQPt->Fill(qpt); + float charge = qpt/std::abs(qpt); + float pt_z = 1./qpt; + + float tglHalf = (sqrt(1 + tgl*tgl) - 1)/tgl; + float eta = - TMath::Log((1 + tglHalf)/(1 - tglHalf)); + float etaAbs = TMath::Abs(eta); + if(etaAbs > 0.8){ + continue; + } + + float p_z = std::abs(pt_z)*TMath::CosH(eta); + if(dedx > 200 && dedx > 600./p_z/p_z){ // this is He + charge *= 2; + } + pt = charge*pt_z; + + if(!(i % (nev/20))){ + printf("%ld/%ld\n",i,nev); + } + + if(ncl < 130){ + continue; + } + + if(std::abs(charge) < 1.5){ + hEff->Fill(pt_z,texp[2] > 1); + } else { + hEff2->Fill(pt_z,texp[2] > 1); + } + + if(texp[2] < 1){ + continue; + } + + float dvtime = vtime*1E6 - tof; + if(std::abs(dvtime) > 20E3){ +// continue; + } + + bool multiHitX = mask&1; + bool multiHitZ = mask&2; + bool badDy = mask&4; + bool multistrip = mask&8; + bool notinpad = mask&16; + bool chiGT3 = mask&32; + bool chiGT5 = mask&64; + bool hasT0sameBC = mask&128; + bool hasT0_1BCbefore = mask&256; + bool hasT0_2BCbefore = mask&512; + + if(badDy || notinpad){ + continue; + } + + if(!hasT0sameBC){ + continue; + } + + if(hasT0_1BCbefore || hasT0_2BCbefore){ + continue; + } + + if(multiHitX || multiHitZ){ + continue; + } + + if(multiHitZ){ + tof -= 0; + } + if(multiHitX){ + tof -= 0; + } + + double t0 = int((tof-10000) * bctimeInv)*bctime; + double beta = texp[0]/(tof-t0); + + int ibinY = hMism->GetXaxis()->FindBin(eta); + float loverc = hEtaMism->GetBinContent(ibinY); + float tmism = hMismPro[ibinY-1]->GetRandom() + loverc; + while(tmism > texp[6] + 2000){ + tmism = hMismPro[ibinY-1]->GetRandom() + loverc; + } + double betaM = texp[0]/(tmism); + if(std::abs(charge) < 1.5){ + float m_z = p_z/beta*sqrt(std::abs(1-beta*beta))*charge; + float nsigmComp = (dedx - eloss->Interpolate(p_z/std::abs(m_z)))/10; + if((nsigmComp > -3 && nsigmComp < 5) || 1){//std::abs(nsigmComp) < 5){ + hPi->Fill(pt_z, tof - t0 - texp[2]); + hKa->Fill(pt_z, tof - t0 - texp[3]); + hPr->Fill(pt_z, tof - t0 - texp[4]); + if(eta > minEta && eta < maxEta){ + int ibin = int((eta - minEta)/etaBin); + int ibin2 = 15 - int((eta - minEta)/etaBin); + + hPiEta[ibin]->Fill(p_z*charge, tof - t0 - texp[2]); + hPiEta[ibin2]->Fill(p_z*charge, tof - t0 - texp[2]); + hPiEtaM[ibin]->Fill(p_z*charge, tmism - texp[2]); + hPiEtaM[ibin2]->Fill(p_z*charge, tmism - texp[2]); + + hKaEta[ibin]->Fill(p_z*charge, tof - t0 - texp[3]); + hKaEta[ibin2]->Fill(p_z*charge, tof - t0 - texp[3]); + hKaEtaM[ibin]->Fill(p_z*charge, tmism - texp[3]); + hKaEtaM[ibin2]->Fill(p_z*charge, tmism - texp[3]); + + hPrEta[ibin]->Fill(p_z*charge, tof - t0 - texp[4]); + hPrEta[ibin2]->Fill(p_z*charge, tof - t0 - texp[4]); + hPrEtaM[ibin]->Fill(p_z*charge, tmism - texp[4]); + hPrEtaM[ibin2]->Fill(p_z*charge, tmism - texp[4]); + + hDeEta[ibin]->Fill(p_z*charge, tof - t0 - texp[5]); + hDeEta[ibin2]->Fill(p_z*charge, tof - t0 - texp[5]); + hDeEtaM[ibin]->Fill(p_z*charge, tmism - texp[5]); + hDeEtaM[ibin2]->Fill(p_z*charge, tmism - texp[5]); + } + + if(std::abs(dedx - eloss->Interpolate(p_z/1.875)) < 15){ + hDe->Fill(pt_z, tof - t0 - texp[5]); + } + if(std::abs(dedx - eloss->Interpolate(p_z/2.809)) < 15){ + hTr->Fill(pt_z, tof - t0 - texp[6]); + } + if(std::abs(tof - t0 - texp[2]) < 200 && std::abs(dedx- eloss->Interpolate(TMath::Min(p_z/0.139,6.))) < eloss->Interpolate(TMath::Min(p_z/0.139,6.))*tpcthr){ + hPiTPC->Fill(p_z*charge, dedx- eloss->Interpolate(p_z/0.139)); + hDeDx->Fill(p_z/0.139*charge,dedx); + hDeDxExp->Fill(p_z/0.139,dedx); + } + if(std::abs(tof - t0 - texp[3]) < 200 && std::abs(dedx- eloss->Interpolate(p_z/0.493)) < eloss->Interpolate(p_z/0.493)*tpcthr && p_z < 2E6){ + hKaTPC->Fill(p_z*charge, dedx- eloss->Interpolate(p_z/0.493)); + hDeDx->Fill(p_z/0.493*charge,dedx); + hDeDxExp->Fill(p_z/0.493,dedx); + } + if(std::abs(tof - t0 - texp[4]) < 200 && std::abs(dedx- eloss->Interpolate(p_z/0.938)) < eloss->Interpolate(p_z/0.938)*tpcthr && p_z < 3.5E6){ + hPrTPC->Fill(p_z*charge, dedx- eloss->Interpolate(p_z/0.938)); + hDeDx->Fill(p_z/0.938*charge,dedx); + hDeDxExp->Fill(p_z/0.938,dedx); + } + if(std::abs(tof - t0 - texp[5]) < 1000 && std::abs(dedx- eloss->Interpolate(p_z/1.875)) < eloss->Interpolate(p_z/1.875)*tpcthr && p_z < 1.5E6){ + hDeTPC->Fill(p_z*charge, dedx- eloss->Interpolate(p_z/1.875)); + hDeDx->Fill(p_z/1.875*charge,dedx); + hDeDxExp->Fill(p_z/1.875,dedx); + } + if(std::abs(tof - t0 - texp[6]) < 1000 && std::abs(dedx- eloss->Interpolate(p_z/2.808)) < eloss->Interpolate(p_z/2.808)*tpcthr && p_z < 1.5E6){ + hTrTPC->Fill(p_z*charge, dedx- eloss->Interpolate(p_z/2.808)); + hDeDx->Fill(p_z/2.808*charge,dedx); + hDeDxExp->Fill(p_z/2.808,dedx); + } + hMass->Fill(m_z, pt); + // if(std::abs(p_z/m_z) > 1 || dedx > 70) + hDeDxTh->Fill(p_z/m_z,dedx - eloss->Interpolate(p_z/std::abs(m_z))); + if(std::abs(p_z)<2){ + hDeDxThM->Fill(m_z,dedx - eloss->Interpolate(p_z/std::abs(m_z))); + } + hDeDxThMr->Fill(m_z,dedx / eloss->Interpolate(p_z/std::abs(m_z))); + } + m_z = p_z/betaM*sqrt(std::abs(1-betaM*betaM))*charge; + nsigmComp = (dedx - eloss->Interpolate(p_z/std::abs(m_z)))/5.86; + if((nsigmComp > -3 && nsigmComp < 5) || 1){//std::abs(nsigmComp) < 5){ + hPiM->Fill(pt_z, tmism- texp[2]); + hKaM->Fill(pt_z, tmism - texp[3]); + hPrM->Fill(pt_z, tmism - texp[4]); + hDeM->Fill(pt_z, tmism - texp[5]); + hTrM->Fill(pt_z, tmism - texp[6]); + hMassM->Fill(m_z, pt); + } + } else { + float m_z = p_z/beta*sqrt(std::abs(1-beta*beta))*charge; + hHe->Fill(pt_z*2, tof - t0 - texp[7]); + hAl->Fill(pt_z*2, tof - t0 - texp[8]); + if(std::abs(tof - t0 - texp[7]) < 1000 || 1){ + hHeTPC->Fill(p_z*charge, dedx- eloss->Interpolate(p_z*2/2.81)*5); + } + if(std::abs(tof - t0 - texp[8]) < 1000 || 1){ + hAlTPC->Fill(p_z*charge, dedx- eloss->Interpolate(p_z*2/3.73)*5); + } + + hMass2->Fill(m_z, pt); + hDeDxThMr->Fill(m_z,dedx / eloss->Interpolate(p_z/std::abs(m_z)*2)); + } + + + } + + TFile *fout = new TFile("summary.root","RECREATE"); + hEff->Write(); + hEff2->Write(); + hQPt->Write(); + hPi->Write(); + hKa->Write(); + hPr->Write(); + hDe->Write(); + hTr->Write(); + hHe->Write(); + hAl->Write(); + hPiM->Write(); + hKaM->Write(); + hPrM->Write(); + hDeM->Write(); + hTrM->Write(); + hMass->Write(); + hMass2->Write(); + hMassM->Write(); + hDeDx->Write(); + hDeDxTh->Write(); + hDeDxThM->Write(); + hDeDxThMr->Write(); + hPiTPC->Write(); + hKaTPC->Write(); + hPrTPC->Write(); + hDeTPC->Write(); + hTrTPC->Write(); + hHeTPC->Write(); + hAlTPC->Write(); + for(int i=0; i < nbinEta; i++){ + hPiEta[i]->Write(); + hKaEta[i]->Write(); + hPrEta[i]->Write(); + hDeEta[i]->Write(); + hPiEtaM[i]->Write(); + hKaEtaM[i]->Write(); + hPrEtaM[i]->Write(); + hDeEtaM[i]->Write(); + } + fout->Close(); + + /* + TFile *expdedx = new TFile("expDeDx.root","RECREATE"); + hDeDxExp->Write(); + expdedx->Close(); + */ +} diff --git a/calibration/checkmap/README b/calibration/checkmap/README new file mode 100644 index 0000000..4e8d98c --- /dev/null +++ b/calibration/checkmap/README @@ -0,0 +1,16 @@ +# download FEELIGHT and channelCalib via timestamp (e.g. run 564400) +o2-ccdb-downloadccdbfile -p TOF/Calib/FEELIGHT -t 1751772182000 +o2-ccdb-downloadccdbfile -p TOF/Calib/ChannelCalib -t 1751772182000 + +# download QC with hit map +alien.py find /alice/data/2025/LHC25ae/564400/cpass0/ QC_fullrun.root|awk '{print "alien.py cp",$1,"file:QC.root"}'|bash + +# run check +root checkmap.C +# or (if you have QC) +root checkmapqc.C + +# if you want to fix the map by masking additional modules you can edit the macro (masking method) and switch off what you want following examples +root checkmap.C +# this will produce a new snapshot2.root in TOF/Calib/FEELIGHT dir + diff --git a/calibration/checkmap/checkmap.C b/calibration/checkmap/checkmap.C new file mode 100644 index 0000000..e48f2cb --- /dev/null +++ b/calibration/checkmap/checkmap.C @@ -0,0 +1,103 @@ +bool checkNew = false; + +bool masking(int ch){ + if(checkNew){ + return false; + } + + int ech = o2::tof::Geo::getECHFromCH(ch); + int crate = o2::tof::Geo::getCrateFromECH(ech); + int indexes[5]; + o2::tof::Geo::getVolumeIndices(ch, indexes); + int& module = indexes[1]; + + // for masking one crate (e.g. crate = 53) +// if(crate == 53){ +// return true; +// } + + // for masking one module (e.g. crate=53, module=1 +// if(crate == 52 && indexes[1] == 1){// && (module == 0 && module == 1)){ +// return true; +// } + + return false; +} + +void checkmap(){ + TFile *f = new TFile("TOF/Calib/ChannelCalib/snapshot.root"); + o2::dataformats::CalibTimeSlewingParamTOF *o = (o2::dataformats::CalibTimeSlewingParamTOF *) f->Get("ccdb_object"); + int nch = 91*96*18; + int ninsec = 91*96; + int nact=0; + int nact2=0; + int nact3=0; + TH2F *h = new TH2F("h","calibrated;sector;",48*18,0,18,91*2,0,91); + TH2F *h2 = new TH2F("h2","FEE map;sector;",48*18,0,18,91*2,0,91); + TH2F *h3 = new TH2F("h2","FEE map && calibrated;sector;",48*18,0,18,91*2,0,91); + + TFile *f2; + if(! checkNew){ + f2 = new TFile("TOF/Calib/FEELIGHT/snapshot.root"); + } else { + f2 = new TFile("TOF/Calib/FEELIGHT/snapshot2.root"); + } + + o2::tof::TOFFEElightInfo *o2 = (o2::tof::TOFFEElightInfo *) f2->Get("ccdb_object"); + + + for(int i=0; i < nch; i++){ + float sector = i/ninsec; + float istrip = (i%ninsec) / 96; + int chLoc = i%96; + int padz = chLoc/48; + int padx = chLoc%48; + sector += (padx+0.5)/48.; + istrip += (padz+0.5)/2.; + + if(masking(i)){ + o2->mChannelEnabled[i] = false; + } + + bool isGoodFEE = o2->getChannelEnabled(i); + bool isGoodCal = ! o->isProblematic(i); + + if(isGoodFEE){ + h2->Fill(sector,istrip); + nact2++; + } + if(isGoodCal){ + h->Fill(sector,istrip); + nact++; + } + + if(isGoodCal && isGoodFEE){ + h3->Fill(sector,istrip); + nact3++; + } + + } + TCanvas *c = new TCanvas; + c->Divide(3,1); + c->cd(1); + h2->Draw("colz"); + h2->SetStats(0); + c->cd(2); + h->Draw("colz"); + h->SetStats(0); + c->cd(3); + h3->Draw("colz"); + h3->SetStats(0); + printf("nActiveFEE = %d\n",nact2); + printf("nActiveCal = %d\n",nact); + printf("nActiveFEE&Cal = %d\n",nact); + + new TCanvas; + h3->Draw("colz"); + + if(! checkNew){ + TFile *fout = new TFile("TOF/Calib/FEELIGHT/snapshot2.root","RECREATE"); + fout->WriteObjectAny(o2, o2->Class(), "ccdb_object"); + fout->Close(); + } +} diff --git a/calibration/checkmap/checkmap.sh b/calibration/checkmap/checkmap.sh new file mode 100755 index 0000000..ce50d0f --- /dev/null +++ b/calibration/checkmap/checkmap.sh @@ -0,0 +1,13 @@ +#!/bin/bash +export run=$2 +export period=$1 + +export ts=$(cat timestamps|grep $run|awk '{print $2}') +o2-ccdb-downloadccdbfile -p TOF/Calib/FEELIGHT -t $ts +o2-ccdb-downloadccdbfile -p TOF/Calib/ChannelCalib -t $ts + +# download QC with hit map +rm QC.root +alien.py find /alice/data/2025/$period/$run/cpass0/ QC_fullrun.root|awk '{print "alien.py cp",$1,"file:QC.root"}'|bash + +root -b -q -l checkmapqc.C\($run,\"$period\"\) diff --git a/calibration/checkmap/checkmapLocal.sh b/calibration/checkmap/checkmapLocal.sh new file mode 100755 index 0000000..aef8f1f --- /dev/null +++ b/calibration/checkmap/checkmapLocal.sh @@ -0,0 +1,14 @@ +#!/bin/bash +export run=$2 +export period=$1 + +export ts=$(cat timestamps|grep $run|awk '{print $2}') +o2-ccdb-downloadccdbfile -p TOF/Calib/FEELIGHT -t $ts +#o2-ccdb-downloadccdbfile -p TOF/Calib/ChannelCalib -t $ts +cp TOF/Calib/ChannelCalib/${run}_ts.root TOF/Calib/ChannelCalib/snapshot.root + +# download QC with hit map +rm QC.root +alien.py find /alice/data/2025/$period/$run/cpass0/ QC_fullrun.root|awk '{print "alien.py cp",$1,"file:QC.root"}'|bash + +root -b -q -l checkmapqc.C\($run,\"$period\"\) diff --git a/calibration/checkmap/checkmapqc.C b/calibration/checkmap/checkmapqc.C new file mode 100644 index 0000000..5375b46 --- /dev/null +++ b/calibration/checkmap/checkmapqc.C @@ -0,0 +1,190 @@ +bool checkNew = false; +bool hitmap = true; + +bool masking(int ch){ + if(checkNew){ + return false; + } + + int ech = o2::tof::Geo::getECHFromCH(ch); + int crate = o2::tof::Geo::getCrateFromECH(ech); + int indexes[5]; + o2::tof::Geo::getVolumeIndices(ch, indexes); + int& module = indexes[1]; + + // for masking one crate (e.g. crate = 53) +// if(crate == 53){ +// return true; +// } + + // for masking one module (e.g. crate=53, module=1 +// if(crate == 52 && indexes[1] == 1){// && (module == 0 && module == 1)){ +// return true; +// } + + return false; +} + +void checkmapqc(int run = 0, const char *period="no"){ + TFile *f = new TFile("TOF/Calib/ChannelCalib/snapshot.root"); + o2::dataformats::CalibTimeSlewingParamTOF *o = (o2::dataformats::CalibTimeSlewingParamTOF *) f->Get("ccdb_object"); + int nch = 91*96*18; + int ninsec = 91*96; + int nact=0; + int nact2=0; + int nact3=0; + TH2F *h = new TH2F("h","calibrated;sector;",48*18,0,18,91*2,0,91); + TH2F *h2 = new TH2F("h2","FEE map;sector;",48*18,0,18,91*2,0,91); + TH2F *h3 = new TH2F("h3","FEE map && calibrated;sector;",48*18,0,18,91*2,0,91); + + TFile *f2; + if(! checkNew){ + f2 = new TFile("TOF/Calib/FEELIGHT/snapshot.root"); + } else { + f2 = new TFile("TOF/Calib/FEELIGHT/snapshot2.root"); + } + + o2::tof::TOFFEElightInfo *oo2 = (o2::tof::TOFFEElightInfo *) f2->Get("ccdb_object"); + + TH2F *hhit = nullptr; + + if(hitmap){ + TFile *fqc = TFile::Open("QC.root"); + if(fqc){ + o2::quality_control::core::MonitorObjectCollection *mon = (o2::quality_control::core::MonitorObjectCollection *) fqc->Get("int/TOF/Digits"); + o2::quality_control::core::MonitorObject *m1 = (o2::quality_control::core::MonitorObject *) mon->FindObject("HitMapNoiseFiltered"); + hhit = (TH2F *) m1->getObject(); + } else { + hitmap = false; + } + } + + for(int i=0; i < nch; i++){ + float sector = i/ninsec; + float istrip = (i%ninsec) / 96; + int chLoc = i%96; + int padz = chLoc/48; + int padx = chLoc%48; + sector += (padx+0.5)/48.; + istrip += (padz+0.5)/2.; + + if(masking(i)){ + oo2->mChannelEnabled[i] = false; + } + + bool isGoodFEE = oo2->getChannelEnabled(i); + bool isGoodCal = ! o->isProblematic(i); + + if(isGoodFEE){ + h2->Fill(sector,istrip); + nact2++; + } + if(isGoodCal){ + h->Fill(sector,istrip); + nact++; + } + + if(isGoodCal && isGoodFEE){ + h3->Fill(sector,istrip); + nact3++; + } + + } + TCanvas *c = new TCanvas; + c->Divide(3,1); + c->cd(1); + h2->Draw("colz"); + h2->SetStats(0); + c->cd(2); + h->Draw("colz"); + h->SetStats(0); + c->cd(3); + h3->Draw("colz"); + h3->SetStats(0); + printf("nActiveFEE = %d\n",nact2); + printf("nActiveCal = %d\n",nact); + printf("nActiveFEE&Cal = %d\n",nact3); + + TCanvas *c2 = new TCanvas; + c2->Divide(1 + hitmap,1); + c2->cd(1); + h3->Draw("colz"); + + if(! hhit){ + return; + } + + c2->cd(2); + int nbinsx = hhit->GetNbinsX(); + int nbinsy = hhit->GetNbinsY(); + int rebinx = 48*18 / nbinsx; + int rebiny = 91*2 / nbinsy; + h3->RebinX(rebinx); + h3->RebinY(rebiny); + h3->Divide(h3); + h->RebinX(rebinx); + h->RebinY(rebiny); + h->Divide(h); + h2->RebinX(rebinx); + h2->RebinY(rebiny); + h2->Divide(h2); + + hhit->Scale(nbinsx*nbinsy*1./hhit->Integral()); + for(int i=0; i < nbinsx; i++){ + for(int j=0; j < nbinsy; j++){ + if(hhit->GetBinContent(i,j) < 0.1){ + hhit->SetBinContent(i,j,0); + } + } + } + + + hhit->Divide(hhit); + // combine with problematics + hhit->Multiply(h); + + hhit->Draw("colz"); + + printf("nActiveHitMap (extrapolated) = %d - ratio wrt expected = %.3f\n",int(hhit->Integral() * nact3 / h3->Integral()), hhit->Integral() / h3->Integral()); + + hhit->SetTitle(Form("%s %d: ratio wrt expected = %.4f",period,run,hhit->Integral() / h3->Integral())); + system(Form("echo %s %d %.4f>>checkmap_summary",period,run,hhit->Integral() / h3->Integral())); + + hhit->Add(h3, -1); + + TFile *fres = new TFile(Form("checkmap/%s_%d.root",period,run),"RECREATE"); + c2->Write(); + fres->Close(); + + if(! checkNew){ + // re-masking + for(int i=0; i < nch; i++){ + float sector = i/ninsec; + float istrip = (i%ninsec) / 96; + int chLoc = i%96; + int padz = chLoc/48; + int padx = chLoc%48; + sector += (padx+0.5)/48.; + istrip += (padz+0.5)/2.; + + if(hhit->GetBinContent(hhit->GetXaxis()->FindBin(sector),hhit->GetYaxis()->FindBin(istrip)) < 0){ +// printf("switch off %.2f %.2f\n",sector,istrip); + oo2->mChannelEnabled[i] = false; + } + else if(hhit->GetBinContent(hhit->GetXaxis()->FindBin(sector),hhit->GetYaxis()->FindBin(istrip)) > 0.1){ +// printf("switch off %.2f %.2f\n",sector,istrip); + oo2->mChannelEnabled[i] = true; + } + } + + TFile *fout; + if(run > 0){ + fout = new TFile(Form("TOF/Calib/FEELIGHT/%d.root",run),"RECREATE"); + } else { + fout = new TFile("TOF/Calib/FEELIGHT/snapshot2.root","RECREATE"); + } + fout->WriteObjectAny(oo2, oo2->Class(), "ccdb_object"); + fout->Close(); + } + +}