#include "AliLog.h"
#include "AliTRDCalibChamberStatus.h"
-#include "AliTRDdEdxUtils.h"
+#include "AliTRDdEdxBaseUtils.h"
+#include "AliTRDdEdxCalibHistArray.h"
+#include "AliTRDdEdxCalibUtils.h"
#include "AliTRDCalibTask.h"
if(fPH2dTest) delete fPH2dTest;
if(fLinearVdriftTest) delete fLinearVdriftTest;
if(IsPHQon()){
- AliTRDdEdxUtils::DeleteCalibHist();
+ AliTRDdEdxCalibUtils::DeleteHistArray();
}
if(fCalDetGain) delete fCalDetGain;
if(IsPHQon()){
printf("\n AliTRDCalibTask PHQ is on!! \n\n");
- AliTRDdEdxUtils::PrintControl();
- AliTRDdEdxUtils::IniCalibHist(fListHist, kTRUE);
+ AliTRDdEdxBaseUtils::PrintControl();
+ AliTRDdEdxCalibUtils::IniHistArray(fListHist, kTRUE);
}
else{
printf("\n AliTRDCalibTask PHQ is off!! \n\n");
}
if(IsPHQon()){
- const Double_t mag = AliTRDdEdxUtils::IsExBOn() ? fESD->GetMagneticField() : -1;
- const Int_t charge = AliTRDdEdxUtils::IsExBOn() ? fkEsdTrack->Charge() : -1;
- const Double_t toTPCscale = AliTRDdEdxUtils::GetCalibTPCscale(fkEsdTrack->GetTPCncls(), fkEsdTrack->GetTPCsignal());
+ const Double_t mag = AliTRDdEdxBaseUtils::IsExBOn() ? fESD->GetMagneticField() : -1;
+ const Int_t charge = AliTRDdEdxBaseUtils::IsExBOn() ? fkEsdTrack->Charge() : -1;
+ const Double_t toTPCscale = AliTRDdEdxCalibUtils::GetCalibTPCscale(fkEsdTrack->GetTPCncls(), fkEsdTrack->GetTPCsignal());
if(toTPCscale>0){
- AliTRDdEdxUtils::FillCalibHist(fTrdTrack, 0, mag, charge, toTPCscale);
+ AliTRDdEdxCalibUtils::FillHist(fTrdTrack, 0, mag, charge, toTPCscale);
}
}
#include "AliTRDCommonParam.h"
#include "AliCDBManager.h"
#include "AliCDBEntry.h"
-#include "AliTRDdEdxUtils.h"
+#include "AliTRDdEdxBaseUtils.h"
+#include "AliTRDdEdxCalibHistArray.h"
+#include "AliTRDdEdxCalibUtils.h"
ClassImp(AliTRDPreprocessorOffline)
if(fSparse) delete fSparse;
if(IsPHQon()){
- AliTRDdEdxUtils::DeleteCalibHist();
- AliTRDdEdxUtils::DeleteCalibObj();
+ AliTRDdEdxCalibUtils::DeleteHistArray();
+ AliTRDdEdxCalibUtils::DeleteObjArray();
}
if(fAliTRDCalibraVdriftLinearFit) delete fAliTRDCalibraVdriftLinearFit;
if(IsPHQon()){
printf("\n AliTRDPreprocessorOffline PHQ on!!\n\n");
- AliTRDdEdxUtils::PrintControl();
+ AliTRDdEdxBaseUtils::PrintControl();
CalibPHQ(file, startRunNumber, endRunNumber, ocdbStorage);
}
else{
// read calibration entries from file
//
- return AliTRDdEdxUtils::ReadCalibHist(fileName, fNameList);
+ return AliTRDdEdxCalibUtils::ReadHistArray(fileName, fNameList);
}
//___________________________________________________________________________________________________________________
//
//Produce PHQ calibration results
//
- for(Int_t iter=0; iter<8; iter++){
- THnF *hi = (THnF*) AliTRDdEdxUtils::GetHistPHQ()->At(iter);
- TObjArray *obji = AliTRDdEdxUtils::GetCalibObj(hi, startRunNumber);
+ for(Int_t iter=0; iter<AliTRDdEdxCalibUtils::GetHistArray()->GetSize(); iter++){
+ THnBase *hi = (THnBase*) AliTRDdEdxCalibUtils::GetHistAt(iter);
+ TObjArray *obji = AliTRDdEdxCalibUtils::HistToObj(hi, startRunNumber);
//printf("test analyze %s\n", obji->GetName());
- AliTRDdEdxUtils::GetObjPHQ()->AddAt(obji, iter);
+ AliTRDdEdxCalibUtils::GetObjArray()->AddAt(obji, iter);
}
- fCalibObjects->AddAt(AliTRDdEdxUtils::GetObjPHQ(), kPHQ);
+ fCalibObjects->AddAt(AliTRDdEdxCalibUtils::GetObjArray(), kPHQ);
return kTRUE;
}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//
+// TRD dEdx base utils
+// xx
+// xx
+// xx
+// xx
+//
+// Xianguo Lu
+// lu@physi.uni-heidelberg.de
+// Xianguo.Lu@cern.ch
+//
+//
+
+#include "TF1.h"
+#include "TFile.h"
+#include "TH1D.h"
+#include "TH2D.h"
+#include "THnSparse.h"
+#include "TMath.h"
+#include "TMatrixD.h"
+#include "TMinuit.h"
+#include "TObjArray.h"
+#include "TRandom3.h"
+#include "TStopwatch.h"
+#include "TVectorD.h"
+
+#include "TTreeStream.h"
+
+#include "AliCDBId.h"
+#include "AliCDBMetaData.h"
+#include "AliCDBStorage.h"
+#include "AliESDEvent.h"
+#include "AliESDfriendTrack.h"
+#include "AliESDtrack.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDCalROC.h"
+#include "AliTRDtrackV1.h"
+
+#include "AliTRDdEdxBaseUtils.h"
+
+#define EPSILON 1e-12
+
+Double_t AliTRDdEdxBaseUtils::fgQ0Frac = 0.3;
+Double_t AliTRDdEdxBaseUtils::fgQ1Frac = 0.5;
+Double_t AliTRDdEdxBaseUtils::fgTimeBinCountCut = 0.0;
+Int_t AliTRDdEdxBaseUtils::fgCalibTPCnclsCut = 70;
+Bool_t AliTRDdEdxBaseUtils::fgExBOn = kTRUE;
+Bool_t AliTRDdEdxBaseUtils::fgPadGainOn = kTRUE;
+Double_t AliTRDdEdxBaseUtils::fgQScale = 45;
+
+//===================================================================================
+// Math and Histogram
+//===================================================================================
+void AliTRDdEdxBaseUtils::GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres)
+{
+ //
+ //counts of hh is sorted
+ //
+
+ for(Int_t ii=0; ii<ncut; ii++){
+ cuts[ii] = -999;
+ }
+
+ Int_t nsel = 0;
+ const Int_t nbin = hh->GetNbinsX();
+ Double_t datas[nbin];
+ for(Int_t ii=1; ii<=nbin; ii++){
+ const Double_t res = hh->GetBinContent(ii);
+ if(res<thres){
+ continue;
+ }
+
+ datas[nsel] = res;
+ nsel++;
+ }
+ if(!nsel)
+ return;
+
+ Int_t id[nsel];
+ TMath::Sort(nsel, datas, id, kFALSE);
+
+ for(Int_t ii=0; ii<ncut; ii++){
+ const Double_t icdf = cdfs[ii];
+ if(icdf<0 || icdf>1){
+ printf("AliTRDdEdxBaseUtils::GetCDFCuts error cdfs[%d] %15f out of range!\n", ii, icdf); exit(1);
+ }
+ cuts[ii] = datas[id[Int_t(icdf*nsel)]];
+ }
+}
+
+Double_t AliTRDdEdxBaseUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr)
+{
+ //
+ //calculate mean (with error) and rms from sum, w2s, nn
+ //if nn=0, mean, error, and rms all = 0
+ //
+
+ Double_t tmean = 0, trms = 0, terr = 0;
+
+ if(nn>EPSILON){
+ tmean = sum/nn;
+
+ const Double_t arg = w2s/nn-tmean*tmean;
+ if(TMath::Abs(arg)<EPSILON){
+ trms = 0;
+ }
+ else{
+ if( arg <0 ){
+ printf("AliTRDdEdxBaseUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1);
+ }
+
+ trms = TMath::Sqrt(arg);
+ }
+
+ terr = trms/TMath::Sqrt(nn);
+ }
+
+ if(grms){
+ (*grms) = trms;
+ }
+
+ if(gerr){
+ (*gerr) = terr;
+ }
+
+ return tmean;
+}
+
+Double_t AliTRDdEdxBaseUtils::TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr, Double_t *wws)
+{
+ //
+ //calculate truncated mean
+ //return <x*w>_{low-high according to x}
+ //
+
+ /*
+ //test->
+ for(Int_t ii=0; ii<nx; ii++){
+ printf("test %d/%d %f\n", ii, nx, xdata[ii]);
+ }
+ //<--test
+ */
+
+ Int_t index[nx];
+ TMath::Sort(nx, xdata, index, kFALSE);
+
+ Int_t nused = 0;
+ Double_t sum = 0;
+ Double_t w2s = 0;
+ const Int_t istart = Int_t (nx*lowfrac);
+ const Int_t istop = Int_t (nx*highfrac);
+
+ //=,< correct, because when low=0, high=1 it is correct
+ for(Int_t ii=istart; ii<istop; ii++){
+ Double_t weight = 1;
+ if(wws){
+ weight = wws[index[ii]];
+ }
+ const Double_t sx = xdata[index[ii]]*weight;
+
+ sum += sx;
+ w2s += sx*sx;
+
+ nused++;
+ //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s);
+
+ }
+
+ return GetMeanRMS(nused, sum, w2s, grms, gerr);
+}
+
+Double_t AliTRDdEdxBaseUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr)
+{
+ //
+ //do truncation on histogram
+ //
+ //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct
+
+ //with under-/over-flow
+ Double_t npreTrunc = 0;
+ for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
+ const Double_t be = hh->GetBinError(itmp);
+ const Double_t bc = hh->GetBinContent(itmp);
+ if(be<EPSILON){
+ if(bc>EPSILON){
+ printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
+ }
+ continue;
+ }
+ npreTrunc += bc*bc/be/be;
+ }
+
+ const Double_t nstart = npreTrunc*lowfrac;
+ const Double_t nstop = npreTrunc*highfrac;
+
+ //with Double_t this should also handle normalized hist
+ Double_t ntot = 0;
+ Double_t nused = 0;
+ Double_t sum = 0;
+ Double_t w2s = 0;
+ for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
+ const Double_t be = hh->GetBinError(itmp);
+ const Double_t bc = hh->GetBinContent(itmp);
+ if(be<EPSILON){
+ if(bc>EPSILON){
+ printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
+ }
+ continue;
+ }
+ const Double_t weight = bc*bc/be/be;
+ ntot+=weight;
+ //<= correct, because when high=1, nstop has to be included
+ if(ntot>nstart && ntot<=nstop){
+
+ const Double_t val = hh->GetBinCenter(itmp);
+ sum += weight*val;
+ w2s += weight*val*val;
+
+ nused += weight;
+
+ //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample);
+ }
+ else if(ntot>nstop){
+ if(itmp>=hh->GetNbinsX()){
+ printf("AliTRDdEdxBaseUtils::TruncatedMean warning hist range too small %s %f %f %d %d, %15f %15f %15f; nused w2s sum set to 0\n", hh->GetName(), hh->GetBinLowEdge(1), hh->GetBinLowEdge(itmp), itmp, hh->GetNbinsX(), hh->GetBinContent(hh->GetNbinsX())/hh->Integral(0,hh->GetNbinsX()+1), hh->GetBinContent(hh->GetNbinsX()), hh->Integral(0,hh->GetNbinsX()+1)); //exit(1);
+ nused = 0;
+ w2s = sum = 0;
+ }
+ break;
+ }
+ }
+
+ return GetMeanRMS(nused, sum, w2s, grms, gerr);
+}
+
+void AliTRDdEdxBaseUtils::FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, const Double_t thres, const Double_t lowfrac, const Double_t highfrac)
+{
+ //
+ //fit slices of hh using truncation
+ //
+
+ const Int_t x0 = hh->GetXaxis()->GetFirst();
+ const Int_t x1 = hh->GetXaxis()->GetLast();
+ const Int_t y0 = hh->GetYaxis()->GetFirst();
+ const Int_t y1 = hh->GetYaxis()->GetLast();
+
+ const Int_t nx = hh->GetNbinsX();
+ const Int_t ny = hh->GetNbinsY();
+ const Double_t xmin = hh->GetXaxis()->GetXmin();
+ const Double_t xmax = hh->GetXaxis()->GetXmax();
+ const Double_t ymin = hh->GetYaxis()->GetXmin();
+ const Double_t ymax = hh->GetYaxis()->GetXmax();
+
+ hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax);
+ hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax);
+ hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax);
+ hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax);
+
+ for(Int_t ix=x0; ix<=x1; ix++){
+ //to speed up
+ const Double_t rawcount = hh->Integral(ix,ix,0, ny+1);
+ if(rawcount<EPSILON){
+ continue;
+ }
+
+ TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
+ Double_t ntot = 0;
+ for(Int_t iy=y0; iy<=y1; iy++){
+ const Double_t be = hh->GetBinError(ix,iy);
+ const Double_t bc = hh->GetBinContent(ix, iy);
+
+ if(be<EPSILON){
+ if(bc>EPSILON){
+ printf("AliTRDdEdxBaseUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
+ }
+ continue;
+ }
+
+ htmp->SetBinContent(iy, bc);
+ htmp->SetBinError(iy, be);
+
+ ntot += (bc/be)*(bc/be);
+
+ //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2));
+ }
+
+ hnor->SetBinContent(ix, ntot);
+ hnor->SetBinError( ix, 0);
+
+ if(ntot<thres || htmp->GetRMS()<EPSILON){
+ delete htmp;
+ continue;
+ }
+
+ //test htmp->Draw();
+ Double_t trms = -999, terr = -999;
+ const Double_t tmean = TruncatedMean(htmp, lowfrac, highfrac, &trms, &terr);
+
+ hmpv->SetBinContent(ix, tmean);
+ hmpv->SetBinError( ix, terr);
+
+ hwid->SetBinContent(ix, trms);
+ hwid->SetBinError( ix, 0);
+
+ hres->SetBinContent(ix, tmean>EPSILON ? trms/tmean:0);
+ hres->SetBinError( ix, 0);
+
+ delete htmp;
+ }
+
+ TH1 *hhs[]={hnor, hmpv, hwid, hres};
+ const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV"};
+ const Int_t nh = sizeof(hhs)/sizeof(TH1*);
+ for(Int_t ii=0; ii<nh; ii++){
+ hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle()));
+ hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle());
+ hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset());
+ hhs[ii]->SetTitle(hh->GetTitle());
+ }
+}
+
+//===================================================================================
+// TRD Analysis Fast Tool
+//===================================================================================
+
+Int_t AliTRDdEdxBaseUtils::GetNtracklet(const AliESDEvent *esd)
+{
+ //
+ //number of trd tracklet in one esd event
+ //
+ const Int_t ntrk0 = esd->GetNumberOfTracks();
+ Int_t ntrdv1=0;
+ for(Int_t ii=0; ii<ntrk0; ii++){
+ ntrdv1 += esd->GetTrack(ii)->GetTRDntracklets();
+ }
+ return ntrdv1;
+}
+
+AliTRDtrackV1 * AliTRDdEdxBaseUtils::GetTRDtrackV1(const AliESDtrack * esdtrack)
+{
+ //
+ //Get TRD friend track
+ //
+
+ AliESDfriendTrack * friendtrk = (AliESDfriendTrack *)esdtrack->GetFriendTrack();
+ if(!friendtrk){
+ //printf("xlulog AliAnalysisTaskCosmicTRD::GetTRDtrack no friend!!\n"); exit(1);
+ return 0x0;
+ }
+
+ TObject *calibObject=0x0;
+ AliTRDtrackV1 * trdtrack=0x0;
+ for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) {
+ if( (trdtrack=dynamic_cast<AliTRDtrackV1*>(calibObject)) )
+ break;
+ }
+
+ return trdtrack;
+}
+
+Bool_t AliTRDdEdxBaseUtils::IsInSameStack(const AliTRDtrackV1 *trdtrack)
+{
+ //
+ // to check if all tracklets are in the same stack, useful in cosmic
+ //
+
+ TVectorD secs(18), stks(5);
+
+ for(Int_t ilayer = 0; ilayer < 6; ilayer++){
+ AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
+ if(!tracklet)
+ continue;
+
+ const Int_t det = tracklet->GetDetector();
+ const Int_t isector = AliTRDgeometry::GetSector(det);
+ const Int_t istack = AliTRDgeometry::GetStack(det);
+
+ secs[isector] = 1;
+ stks[istack] = 1;
+ }
+
+ if(secs.Sum()!=1 || stks.Sum()!=1){
+ return kFALSE;
+ }
+ else
+ return kTRUE;
+}
+
+Bool_t AliTRDdEdxBaseUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
+{
+ //
+ //as function name
+ //
+ isec = istk = -999;
+ mom = -999;
+
+ for(Int_t ilayer = 0; ilayer < 6; ilayer++){
+ AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
+ if(!tracklet)
+ continue;
+
+ const Int_t det = tracklet->GetDetector();
+ isec = AliTRDgeometry::GetSector(det);
+ istk = AliTRDgeometry::GetStack(det);
+
+ mom = tracklet->GetMomentum();
+
+ break;
+ }
+
+ if(isec<0)
+ return kFALSE;
+ else
+ return kTRUE;
+}
+
+//===================================================================================
+// Detector and Data Constant
+//===================================================================================
+Int_t AliTRDdEdxBaseUtils::ToDetector(const Int_t gtb)
+{
+ //
+ //gtb = det*Ntb+itb
+ //
+ return gtb/AliTRDseedV1::kNtb;
+}
+
+Int_t AliTRDdEdxBaseUtils::ToTimeBin(const Int_t gtb)
+{
+ //
+ //gtb = det*Ntb+itb
+ //
+ return gtb%AliTRDseedV1::kNtb;
+}
+
+Int_t AliTRDdEdxBaseUtils::ToSector(const Int_t gtb)
+{
+ //
+ //return sector
+ //
+ return AliTRDgeometry::GetSector(ToDetector(gtb));
+}
+
+Int_t AliTRDdEdxBaseUtils::ToStack(const Int_t gtb)
+{
+ //
+ //return stack
+ //
+ return AliTRDgeometry::GetStack(ToDetector(gtb));
+}
+
+Int_t AliTRDdEdxBaseUtils::ToLayer(const Int_t gtb)
+{
+ //
+ //return layer
+ //
+ return AliTRDgeometry::GetLayer(ToDetector(gtb));
+}
+
+TString AliTRDdEdxBaseUtils::GetRunType(const Int_t run)
+{
+ //
+ //return run type
+ //
+
+ TString type;
+ if(run>=121527 && run<= 126460)//LHC10d
+ type="pp2010LHC10d";
+ else if(run>=126461 && run<= 130930)//LHC10e
+ type="pp2010LHC10e";
+ else if(run>=136782 && run <= 139846)//LHC10h
+ type="PbPb2010LHC10h";
+ else if(run>= 143224 && run<= 143237)//2011Feb
+ type="cosmic2011Feb";
+ else if(run>= 150587 && run<= 154930){
+ type="cosmic2011MayJun";
+
+ TString runstr(Form("%d",run));
+ const TString listrun1kg("154601 154602 154629 154634 154636 154639 154643");
+ if(listrun1kg.Contains(runstr)){
+ type+="1kG";
+ }
+ else{
+ type+="5kG";
+ }
+ }
+ else{
+ type="unknown";
+ }
+
+ type.ToUpper();
+ return type;
+}
+
+void AliTRDdEdxBaseUtils::PrintControl()
+{
+ //
+ //print out control variable
+ //
+ printf("\nAliTRDdEdxBaseUtils::PrintControl Q0Frac %.1f, Q1Frac %.1f, TimeBinCountCut %.2f, CalibTPCnclsCut %d, IsExBOn %d, IsPadGainOn %d, QScale %.2f\n\n", Q0Frac(), Q1Frac(), TimeBinCountCut(), CalibTPCnclsCut(), IsExBOn(), IsPadGainOn(), QScale());
+}
+
+//===================================================================================
+// dEdx Parameterization
+//===================================================================================
+
+Double_t AliTRDdEdxBaseUtils::Q0MeanTRDpp(const Double_t bg)
+{
+ //
+ //truncated Mean Q_{xx} in TRD
+ //
+
+ Double_t par[8];
+ //03132012161150
+ //opt: ppQ0
+par[0]= 2.397001e-01;
+par[1]= 1.334697e+00;
+par[2]= 6.967470e+00;
+par[3]= 9.055289e-02;
+par[4]= 9.388760e+00;
+par[5]= 9.452742e-04;
+par[6]= -1.866091e+00;
+par[7]= 1.403545e+00;
+
+ ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see2.log:hhtype2Q0b2c2 scale 0.428543 at ltbg 0.650000
+ return 0.428543*MeandEdxTR(&bg, par);
+}
+
+Double_t AliTRDdEdxBaseUtils::Q0MeanTRDPbPb(const Double_t bg)
+{
+ //
+ //truncated Mean Q_{xx} in TRD
+ //
+
+ Double_t par[8];
+
+ //03132012161259
+ //opt: PbPbQ0
+par[0]= 1.844912e-01;
+par[1]= 2.509702e+00;
+par[2]= 6.744031e+00;
+par[3]= 7.355123e-02;
+par[4]= 1.166023e+01;
+par[5]= 1.736186e-04;
+par[6]= -1.716063e+00;
+par[7]= 1.611366e+00;
+
+ ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale 0.460994 at ltbg 0.650000
+ return 0.460994*MeandEdxTR(&bg, par);
+}
+
+Double_t AliTRDdEdxBaseUtils::Q1MeanTRDpp(const Double_t bg)
+{
+ //
+ //truncated Mean 1/(1/Q)_{xx} in TRD
+ //
+
+ Double_t par[8];
+
+ //So 4. Mär 13:30:51 CET 2012
+ //opt: trdppQ1
+ par[0]= 2.434646e-01;
+ par[1]= 1.400211e+00;
+ par[2]= 6.937471e+00;
+ par[3]= 7.758118e-02;
+ par[4]= 1.097372e+01;
+ par[5]= 4.297518e-04;
+ par[6]= -1.806266e+00;
+ par[7]= 1.543811e+00;
+
+ //hhtype2Q1b2c2 scale 0.418629 at ltbg 0.650000
+
+ return 0.418629*MeandEdxTR(&bg, par);
+}
+
+Double_t AliTRDdEdxBaseUtils::Q1MeanTRDPbPb(const Double_t bg)
+{
+ //
+ //truncated Mean 1/(1/Q)_{xx} in TRD
+ //
+
+ Double_t par[8];
+
+ //So 4. Mär 13:30:52 CET 2012
+ //opt: trdPbPbQ1
+ par[0]= 2.193660e-01;
+ par[1]= 2.051864e+00;
+ par[2]= 6.825112e+00;
+ par[3]= 6.151693e-02;
+ par[4]= 1.390343e+01;
+ par[5]= 6.010032e-05;
+ par[6]= -1.676324e+00;
+ par[7]= 1.838873e+00;
+
+ //hhtype4Q1b2c2 scale 0.457988 at ltbg 0.650000
+
+ return 0.457988*MeandEdxTR(&bg, par);
+}
+
+Double_t AliTRDdEdxBaseUtils::QMeanTPC(const Double_t bg)
+{
+ //
+ //bethe bloch in TPC
+ //
+
+ Double_t par[5];
+ //Mi 15. Feb 14:48:05 CET 2012
+ //train_2012-02-13_1214.12001, tpcsig
+ par[0]= 4.401269;
+ par[1]= 9.725370;
+ par[2]= 0.000178;
+ par[3]= 1.904962;
+ par[4]= 1.426576;
+
+ return MeandEdx(&bg, par);
+}
+
+Double_t AliTRDdEdxBaseUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
+{
+ //
+ //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
+ //npar = 8 = 3+5
+ //
+ Double_t ptr[4]={0};
+ for(int ii=0; ii<3; ii++){
+ ptr[ii+1]=pin[ii];
+ }
+ return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
+}
+
+Double_t AliTRDdEdxBaseUtils::MeanTR(const Double_t * xx, const Double_t * par)
+{
+ //
+ //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
+ //npar = 4
+ //
+
+ const Double_t bg = xx[0];
+ const Double_t gamma = sqrt(1+bg*bg);
+
+ const Double_t p0 = TMath::Abs(par[1]);
+ const Double_t p1 = TMath::Abs(par[2]);
+ const Double_t p2 = TMath::Abs(par[3]);
+
+ const Double_t zz = TMath::Log(gamma);
+ const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
+
+ return par[0]+tryield;
+}
+
+Double_t AliTRDdEdxBaseUtils::MeandEdx(const Double_t * xx, const Double_t * par)
+{
+ //
+ //ALEPH parametrization for dEdx
+ //npar = 5
+ //
+
+ const Double_t bg = xx[0];
+ const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
+
+ const Double_t p0 = TMath::Abs(par[0]);
+ const Double_t p1 = TMath::Abs(par[1]);
+ const Double_t p2 = TMath::Abs(par[2]);
+ const Double_t p3 = TMath::Abs(par[3]);
+ const Double_t p4 = TMath::Abs(par[4]);
+
+ const Double_t aa = TMath::Power(beta, p3);
+ const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
+
+ //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
+
+ return (p1-aa-bb)*p0/aa;
+}
+
+Double_t AliTRDdEdxBaseUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
+{
+ //
+ //f(x)-> f(y) with y=log10(x)
+ //
+ const Double_t x2[]={TMath::Power(10, xx[0])};
+ return func(x2, par);
+}
+
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+//
+//
+//
+// Xianguo Lu
+// lu@physi.uni-heidelberg.de
+// Xianguo.Lu@cern.ch
+//
+/*
+grep " AliTRDdEdxBaseUtils::" AliTRDdEdxBaseUtils.cxx | grep "=" -v | grep -v "[6]" | grep -v printf |wc
+grep "(" AliTRDdEdxBaseUtils.h | grep ";" | grep -v grep | grep -v ClassDef | grep -v "{" | grep -v typedef | wc
+*/
+
+
+#ifndef ALITRDDEDXBASEUTILS_H
+#define ALITRDDEDXBASEUTILS_H
+
+#ifndef TVECTORD_H
+#include "TVectorD.h"
+#endif
+
+#ifndef THNSPARSE_H
+#include "THnBase.h"
+#endif
+
+#ifndef TTREESTREAM_H
+#include "TTreeStream.h"
+#endif
+
+class TH1D;
+class TH2D;
+class TObjArray;
+
+class AliESDEvent;
+class AliESDtrack;
+class AliTRDcluster;
+class AliTRDtrackV1;
+class AliTRDseedV1;
+
+class AliTRDdEdxBaseUtils
+{
+ public:
+ //===================================================================================
+ // Math and Histogram
+ //===================================================================================
+ static void GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres);
+ static Double_t GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms=0x0, Double_t * gerr=0x0);
+ static Double_t TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms=0x0, Double_t * gerr=0x0, Double_t *wws=0x0);
+ static Double_t TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms=0x0, Double_t * gerr=0x0);
+ static void FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, const Double_t thres, const Double_t lowfrac, const Double_t highfrac);
+
+ //===================================================================================
+ // TRD Analysis Fast Tool
+ //===================================================================================
+ static Int_t GetNtracklet(const AliESDEvent *esd);
+ static AliTRDtrackV1 * GetTRDtrackV1(const AliESDtrack * esdtrack);
+ static Bool_t IsInSameStack(const AliTRDtrackV1 *trdtrack);
+ static Bool_t GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom);
+
+ //===================================================================================
+ // Detector, Data and Control Constant
+ //===================================================================================
+
+ static Int_t NTRDchamber(){return 18*5*6;} //540
+ static Int_t NTRDtimebin(){return NTRDchamber()*31;} //16740
+ static Int_t ToDetector(const Int_t gtb);
+ static Int_t ToTimeBin(const Int_t gtb);
+ static Int_t ToSector(const Int_t gtb);
+ static Int_t ToStack(const Int_t gtb);
+ static Int_t ToLayer(const Int_t gtb);
+
+ static TString GetRunType(const Int_t run);
+
+ static void SetQ0Frac(const Double_t q0){ fgQ0Frac = q0; }
+ static void SetQ1Frac(const Double_t q1){ fgQ1Frac = q1; }
+ static void SetTimeBinCountCut(const Double_t tbc){ fgTimeBinCountCut = tbc; }
+ static void SetCalibTPCnclsCut(const Int_t tpc){ fgCalibTPCnclsCut = tpc; }
+ static void SetExBOn(const Bool_t kon){ fgExBOn = kon; }
+ static void SetPadGainOn(const Bool_t kon){ fgPadGainOn = kon; }
+ static void SetQScale(const Double_t scale){ fgQScale = scale; }
+
+ static Double_t Q0Frac(){return fgQ0Frac;}
+ static Double_t Q1Frac(){return fgQ1Frac;}
+ static Double_t TimeBinCountCut(){return fgTimeBinCountCut;}
+ static Int_t CalibTPCnclsCut(){return fgCalibTPCnclsCut;}
+ static Bool_t IsExBOn(){return fgExBOn;}
+ static Bool_t IsPadGainOn(){return fgPadGainOn;}
+ static Double_t QScale(){return fgQScale;}
+
+ static void PrintControl();
+
+ //===================================================================================
+ // dEdx Parameterization
+ //===================================================================================
+
+ static Double_t MeandEdx(const Double_t * xx, const Double_t * par);
+ static Double_t MeanTR(const Double_t * xx, const Double_t * par);
+ static Double_t MeandEdxTR(const Double_t * xx, const Double_t * par);
+
+ static Double_t QMeanTPC(const Double_t bg);
+ static Double_t Q0MeanTRDpp(const Double_t bg);
+ static Double_t Q1MeanTRDpp(const Double_t bg);
+ static Double_t Q0MeanTRDPbPb(const Double_t bg);
+ static Double_t Q1MeanTRDPbPb(const Double_t bg);
+
+ typedef Double_t (*FFunc)(const Double_t *xx, const Double_t *par);
+
+ static Double_t MeandEdxLogx(const Double_t * xx, const Double_t * par){return ToLogx(MeandEdx, xx, par);}
+ static Double_t MeanTRLogx(const Double_t * xx, const Double_t * par){return ToLogx(MeanTR, xx, par);}
+ static Double_t MeandEdxTRLogx(const Double_t * xx, const Double_t * par){return ToLogx(MeandEdxTR, xx, par);}
+
+ private:
+ //dEdx Parameterization
+ static Double_t ToLogx(FFunc func, const Double_t * xx, const Double_t * par);
+
+ //Control Constant
+ static Double_t fgQ0Frac; //q0frac
+ static Double_t fgQ1Frac; //q1frac
+ static Double_t fgTimeBinCountCut; //tbcut
+ static Int_t fgCalibTPCnclsCut; //tpccut
+ static Bool_t fgExBOn; //exbon
+ static Bool_t fgPadGainOn; //pad gain
+ static Double_t fgQScale; //Qscale
+
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//
+// xx
+// xx
+// xx
+// xx
+// xx
+//
+// Xianguo Lu
+// lu@physi.uni-heidelberg.de
+// Xianguo.Lu@cern.ch
+//
+//
+#include "THnBase.h"
+#include "THnSparse.h"
+#include "TCollection.h"
+
+#include "AliTRDdEdxBaseUtils.h"
+#include "AliTRDdEdxCalibHistArray.h"
+
+ClassImp(AliTRDdEdxCalibHistArray);
+
+AliTRDdEdxCalibHistArray::AliTRDdEdxCalibHistArray(const Bool_t kNoInv):
+ TObjArray(kNoInv ? 4: 8)
+{
+ //
+ //constructor
+ //
+ SetName(GetArrayName());
+ SetOwner(kTRUE);
+
+ const Int_t nbin[2]={AliTRDdEdxBaseUtils::NTRDtimebin(), 11250};
+ const Double_t xmin[2]={0, 0};
+ const Double_t xmax[2]={nbin[0], 20};
+
+ for(Int_t iter=0; iter<GetSize(); iter++){
+ THnBase *hi = new THnSparseS(GetNameAt(iter), "", 2, nbin, xmin, xmax);
+ AddAt(hi, iter);
+ }
+}
+
+AliTRDdEdxCalibHistArray::AliTRDdEdxCalibHistArray(const AliTRDdEdxCalibHistArray &obj):
+ TObjArray(obj)
+{
+ //
+ //copy constructor
+ //
+}
+
+AliTRDdEdxCalibHistArray & AliTRDdEdxCalibHistArray::operator=(const AliTRDdEdxCalibHistArray &obj)
+{
+ //
+ //assignment operator
+ //
+
+ if(&obj == this) return *this;
+
+ TObjArray::operator=(obj);
+
+ return *this;
+}
+
+Long64_t AliTRDdEdxCalibHistArray::Merge(const TCollection* list)
+{
+ //
+ // Merge list of objects (needed by PROOF)
+ //
+
+ if(!list)
+ return 0;
+
+ if(list->IsEmpty())
+ return 1;
+
+ TIterator* iter = list->MakeIterator();
+ TObject* obj = 0;
+
+ Int_t count=0;
+ while((obj = iter->Next()) != 0)
+ {
+ AliTRDdEdxCalibHistArray * entry = dynamic_cast<AliTRDdEdxCalibHistArray*>(obj);
+ if (entry == 0) continue;
+
+ if(GetSize()!= entry->GetSize()){
+ printf("AliTRDdEdxCalibHistArray::Merge GetSize()!= entry->GetSize() %d %d\n", GetSize(), entry->GetSize()); exit(1);
+ }
+
+ for(Int_t ii=0; ii<GetSize(); ii++){
+ THnBase *h0 = (THnBase*) At(ii);
+ THnBase *h1 = (THnBase*) entry->At(ii);
+ h0->Add(h1);
+ }
+
+ count++;
+ }
+
+ return count;
+
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+//
+//
+//
+// Xianguo Lu
+// lu@physi.uni-heidelberg.de
+// Xianguo.Lu@cern.ch
+//
+
+#ifndef ALITRDDEDXCALIBHISTARRAY_H
+#define ALITRDDEDXCALIBHISTARRAY_H
+
+
+#ifndef THNSPARSE_H
+#include "THnBase.h"
+#endif
+
+class TObjArray;
+class TCollection;
+
+class AliTRDdEdxCalibHistArray: public TObjArray
+{
+ public:
+ AliTRDdEdxCalibHistArray(const Bool_t kNoInv=kTRUE);
+ AliTRDdEdxCalibHistArray(const AliTRDdEdxCalibHistArray &obj);
+ AliTRDdEdxCalibHistArray & operator=(const AliTRDdEdxCalibHistArray &obj);
+ virtual ~AliTRDdEdxCalibHistArray(){} //definition {} important for virtual
+ virtual Long64_t Merge(const TCollection* list);
+
+ static TString GetArrayName(){ return "TRDdEdxCalibHistArray"; }
+ static TString GetNameAt(const Int_t iter){ return Form("TRDdEdxCalibHist%d", iter); }
+ static Int_t GetIterator(const Bool_t kinvq, const Double_t mag, const Int_t charge){ return kinvq*4 + (mag>0)*2 + (charge>0); }
+
+
+ THnBase * GetHist(const Bool_t kinvq, const Double_t mag, const Int_t charge){ return (THnBase*) At(GetIterator(kinvq, mag, charge)); }
+
+ private:
+
+ ClassDef(AliTRDdEdxCalibHistArray,1);
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//
+// class to calculate TRD dEdx
+// xx
+// xx
+// xx
+// xx
+//
+// Xianguo Lu
+// lu@physi.uni-heidelberg.de
+// Xianguo.Lu@cern.ch
+//
+//
+
+
+#include "TF1.h"
+#include "TFile.h"
+#include "TH1D.h"
+#include "TH2D.h"
+#include "THnSparse.h"
+#include "TMath.h"
+#include "TMatrixD.h"
+#include "TMinuit.h"
+#include "TObjArray.h"
+#include "TRandom3.h"
+#include "TStopwatch.h"
+#include "TVectorD.h"
+
+#include "TTreeStream.h"
+
+#include "AliCDBId.h"
+#include "AliCDBMetaData.h"
+#include "AliCDBStorage.h"
+#include "AliESDEvent.h"
+#include "AliESDfriendTrack.h"
+#include "AliESDtrack.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDCalROC.h"
+#include "AliTRDtrackV1.h"
+
+#include "AliTRDdEdxBaseUtils.h"
+#include "AliTRDdEdxReconUtils.h"
+#include "AliTRDdEdxCalibHistArray.h"
+#include "AliTRDdEdxCalibUtils.h"
+
+#define EPSILON 1e-12
+
+AliTRDdEdxCalibHistArray * AliTRDdEdxCalibUtils::fgHistArray = 0x0;
+TObjArray * AliTRDdEdxCalibUtils::fgObjArray = 0x0;
+
+//============================================================
+// CalibObj
+//============================================================
+TObjArray * AliTRDdEdxCalibUtils::GetObjArray()
+{
+ //
+ //return fgObjArray, initialized if null
+ //
+
+ if(!fgObjArray){
+ fgObjArray = new TObjArray(8);
+ }
+
+ return fgObjArray;
+}
+
+TObjArray * AliTRDdEdxCalibUtils::GetObj(const Bool_t kinvq, const Double_t mag, const Int_t charge)
+{
+ //
+ //return calib obj
+ //
+ if(!fgObjArray){
+ printf("AliTRDdEdxCalibUtils::GetObjArray(kinvq, mag, charge) error fgObjArray null!!\n"); exit(1);
+ }
+
+ return (TObjArray*) fgObjArray->At(AliTRDdEdxCalibHistArray::GetIterator(kinvq, mag, charge));
+}
+
+void AliTRDdEdxCalibUtils::DeleteObjArray()
+{
+ //
+ //delete calib obj
+ //
+ if(fgObjArray){
+ fgObjArray->SetOwner();
+ delete fgObjArray;
+ fgObjArray = 0x0;
+ }
+}
+
+Bool_t AliTRDdEdxCalibUtils::GenerateDefaultOCDB(const TString path)
+{
+ //
+ //generate default OCDB object PHQ, do like
+ //AliTRDdEdxCalibUtils::GenerateDefaultPHQOCDB("local://./")
+ //
+
+ TObjArray * arr8 = new TObjArray(8);
+ arr8->SetOwner();
+
+ for(Int_t ii=0; ii<8; ii++){
+ TObjArray * arr1 = new TObjArray(1);
+ arr1->SetOwner();
+ TString objn(AliTRDdEdxCalibHistArray::GetNameAt(ii));
+ objn.ReplaceAll("Hist","Obj");
+ arr1->SetName(objn);
+
+ const Int_t nbins = AliTRDdEdxBaseUtils::NTRDtimebin();
+ TVectorD * vec = new TVectorD(nbins);
+ for(Int_t jj=0; jj<nbins; jj++){
+ (*vec)[jj] = 1;
+ }
+ arr1->Add(vec);
+ arr8->Add(arr1);
+ }
+
+ AliCDBMetaData *metaData= new AliCDBMetaData();
+ metaData->SetObjectClassName("TObjArray");
+ metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu");
+
+ AliCDBId id1("TRD/Calib/PHQ", 0, 999999999);
+ AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage(path);
+ gStorage->Put(arr8, id1, metaData);
+
+ delete metaData;
+ delete arr8;
+
+ return kTRUE;
+}
+
+//============================================================
+// CalibHist
+//============================================================
+void AliTRDdEdxCalibUtils::DeleteHistArray()
+{
+ //
+ //delete calib hist
+ //
+ delete fgHistArray;
+ fgHistArray = 0x0;
+}
+
+THnBase * AliTRDdEdxCalibUtils::GetHistAt(const Int_t iter)
+{
+ //
+ //
+ //
+ if(iter<fgHistArray->GetSize())
+ return (THnBase*)fgHistArray->At(iter);
+ else
+ return 0x0;
+}
+
+void AliTRDdEdxCalibUtils::IniHistArray(TList *list, const Bool_t kNoInv)
+{
+ //
+ //initialize calib hist, list should not own the hist, or list->Clear/delete hist should not be called
+ //
+
+ delete fgHistArray;
+ fgHistArray = new AliTRDdEdxCalibHistArray(kNoInv);
+ list->Add(fgHistArray);
+}
+
+Bool_t AliTRDdEdxCalibUtils::ReadHistArray(const TString filename, const TString listname)
+{
+ //
+ //used in AliTRDPreprocessorOffline
+ //read in calib hist from file, only for PHQ
+ //
+
+ //maybe already open by others... don't close
+ TFile fcalib(filename);
+ TObjArray * array = (TObjArray*)fcalib.Get(listname);
+ const TString histname = AliTRDdEdxCalibHistArray::GetArrayName();
+
+ AliTRDdEdxCalibHistArray * tmph=0x0;
+ if(array){
+ tmph = (AliTRDdEdxCalibHistArray *) array->FindObject(histname);
+ }
+ else{
+ tmph = (AliTRDdEdxCalibHistArray *) fcalib.Get(histname);
+ }
+ if(!tmph){
+ printf("AliTRDdEdxCalibUtils::ReadCalibHist warning calib hist not found! %s %s %s\n", filename.Data(), listname.Data(), histname.Data());
+ fcalib.ls();
+ if(array){
+ array->ls();
+ }
+ return kFALSE;
+ }
+
+ delete fgHistArray;
+ fgHistArray = new AliTRDdEdxCalibHistArray(*tmph);
+
+ return kTRUE;
+}
+
+void AliTRDdEdxCalibUtils::FillHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnBase * hcalib, const Double_t scale)
+{
+ //
+ //fill calibration hist
+ //
+ if(!hcalib){printf("AliTRDdEdxCalibUtils::FillCalibHist errro hcalib null!!\n"); exit(1);}
+
+ for(Int_t ii=0; ii<ncls; ii++){
+ const Double_t dq = (*arrayQ)[ii];
+ const Double_t xx = (*arrayX)[ii];
+
+ const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
+ const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
+ const Double_t xmax = hcalib->GetAxis(0)->GetXmax();
+
+ if(xx>=xmax || xx<xmin){
+ printf("AliTRDdEdxCalibUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(), xx, xmin, xmax); exit(1);
+ }
+
+ const Double_t var[]={xx, TMath::Min(dq, qmax)/scale};
+ hcalib->Fill(var);
+ }
+}
+
+void AliTRDdEdxCalibUtils::FillHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale)
+{
+ //
+ //get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled
+ //
+ if(!fgHistArray){
+ printf("AliTRDdEdxCalibUtils::FillHist fgHistArray not initialized!!\n"); exit(1);
+ }
+
+ THnBase * hcalib = fgHistArray->GetHist(kinvq, mag, charge);
+
+ TVectorD arrayQ(200), arrayX(200);
+ const Int_t ncls = AliTRDdEdxReconUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
+ FillHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);
+
+ static Int_t kprint = 100;
+ if(kprint<0){
+ printf("\nAliTRDdEdxCalibUtils::FillHist summary: \n");
+ printf("\nkinvq= %d;\n", kinvq);
+ for(Int_t iq=0; iq<ncls; iq++){
+ printf("arrayX[%3d] = %15.0f; arrayQ[%3d] = %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
+ }
+ printf("\n");
+ }
+ kprint++;
+}
+
+//============================================================
+//
+//============================================================
+
+Double_t AliTRDdEdxCalibUtils::GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig)
+{
+ //
+ //the scale used in calibration
+ //
+
+ if(tpcncls < AliTRDdEdxBaseUtils::CalibTPCnclsCut())
+ return -999;
+
+ if(tpcsig<EPSILON)
+ return -999;
+
+ return tpcsig/120;
+
+}
+
+void AliTRDdEdxCalibUtils::GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean)
+{
+ //
+ //calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output
+ //
+ const Int_t ndet = 540;
+ TObjArray *obj=new TObjArray(ndet);
+ obj->SetOwner();
+ for(Int_t ii=0; ii<ndet; ii++){
+ obj->Add(new TVectorD(AliTRDseedV1::kNtb));
+ }
+
+ //ibin = binlowedge of bin(ibin+1) = the number fills this bin
+ for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
+ const Double_t stat = hnor->GetBinContent(ibin+1);
+ if(stat<EPSILON){
+ continue;
+ }
+
+ const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
+ const Int_t itb = AliTRDdEdxBaseUtils::ToTimeBin(ibin);
+ TVectorD *vec=(TVectorD *)obj->At(idet);
+ (*vec)[itb] = stat;
+ }
+
+ hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
+ for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
+ const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
+ const TVectorD *vec=(TVectorD *)obj->At(idet);
+
+ Int_t nonzero = 0;
+ for(Int_t ii=0; ii<vec->GetNrows(); ii++){
+ if((*vec)[ii]>EPSILON){
+ nonzero++;
+ }
+ }
+
+ Double_t mean = 0;
+ const Double_t lowfrac = 0.6;
+ //if there are too many 0's, reject this chamber by settig mean=rms=0
+ if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){
+ //only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large!
+ mean = AliTRDdEdxBaseUtils::TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
+ }
+
+ hmean->SetBinContent(ibin+1, mean);
+ }
+
+ delete obj;
+}
+
+void AliTRDdEdxCalibUtils::Output(const TList *lin, Int_t run)
+{
+ //
+ //produce calibration objects
+ //
+
+ TString objnames;
+ for(Int_t iter=0; iter<8; iter++){
+ objnames+= AliTRDdEdxCalibHistArray::GetNameAt(iter)+" ";
+ }
+
+ TList * lout = new TList;
+ lout->SetOwner();
+
+ TTreeSRedirector *calibStream = new TTreeSRedirector(Form("TRDCalibStream_%010d.root", run));
+
+ const Int_t nh=lin->GetEntries();
+ for(Int_t ii=0; ii<nh; ii++){
+ const THnBase *hh=(THnBase*)lin->At(ii);
+ const TString hname = hh->GetName();
+ if(!objnames.Contains(hname))
+ continue;
+
+ TObjArray * cobj0 = HistToObj(hh, run, lout, calibStream);
+ lout->Add(cobj0);
+ }
+
+ //lout->ls();
+
+ //=============================================================
+ //=============================================================
+
+ TFile *fout=new TFile(Form("TRDCalibObj_%010d.root", run),"recreate");
+ fout->cd();
+ const Int_t nout=lout->GetEntries();
+ for(Int_t ii=0; ii<nout; ii++){
+ const TString oname = lout->At(ii)->GetName();
+ if(oname.Contains("Obj")){
+ TObjArray * cobj = (TObjArray*) lout->At(ii);
+ cobj->Write(oname, TObjArray::kSingleKey);
+ }
+ }
+ fout->Save();
+ fout->Close();
+ delete fout;
+
+ fout=new TFile(Form("TRDCalibList_%010d.root", run),"recreate");
+ fout->cd();
+ lin->Write();
+ lout->Write();
+ fout->Save();
+ fout->Close();
+ delete fout;
+
+ delete calibStream;
+
+ /*
+ http://root.cern.ch/root/html/TH1.html
+ When an histogram is created, a reference to it is automatically added to the list of in-memory objects for the current file or directory. This default behaviour can be changed by:
+
+ h->SetDirectory(0); for the current histogram h
+ TH1::AddDirectory(kFALSE); sets a global switch disabling the reference
+
+ When the histogram is deleted, the reference to it is removed from the list of objects in memory. When a file is closed, all histograms in memory associated with this file are automatically deleted.
+ */
+ delete lout;
+}
+
+TObjArray* AliTRDdEdxCalibUtils::HistToObj(const THnBase *hh, Int_t run, TList *lout, TTreeSRedirector *calibStream)
+{
+ //
+ //produce calibration objects
+ //
+
+ const TString hname = hh->GetName();
+ const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;
+
+ //----------------------------------------
+ // Define nbin, tag, cobj0
+ //----------------------------------------
+ const Int_t nbin = AliTRDdEdxBaseUtils::NTRDtimebin();
+
+
+ TString tag(hname);
+ tag.ReplaceAll("Hist","Obj");
+
+ TObjArray * cobj0 = new TObjArray(1);
+ cobj0->SetOwner();
+ cobj0->SetName(tag);
+ cobj0->Add(new TVectorD(nbin));
+
+ //----------------------------------------
+ // Define lowFrac, highFrac
+ //----------------------------------------
+ Double_t lowFrac = -999, highFrac = -999;
+ if(!kinvq) {
+ lowFrac = 0.01; highFrac = AliTRDdEdxBaseUtils::Q0Frac();
+ }
+ else{
+ lowFrac = AliTRDdEdxBaseUtils::Q1Frac(); highFrac = 0.99;
+ }
+
+ //----------------------------------------
+ // Get analysis result
+ //----------------------------------------
+ TH1::AddDirectory(kFALSE);//important!
+ TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
+ TH2D *hpj = hh->Projection(1,0);
+ AliTRDdEdxBaseUtils::FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
+ GetPHCountMeanRMS(hnor, htrdphmean);
+ if(lout){
+ lout->Add(htrdphmean);
+ }
+ delete hpj;
+
+ if(lout){
+ lout->Add(hnor);
+ lout->Add(hmpv);
+ lout->Add(hwid);
+ lout->Add(hres);
+ }
+
+ //----------------------------------------
+ // Define Counter
+ //----------------------------------------
+ TVectorD *countDet=0x0;
+ TObjArray *countSSL=0x0;
+
+ if(!kinvq){
+ countDet=new TVectorD(540);
+ countSSL=new TObjArray(90);//SectorStackLayer
+ countSSL->SetOwner();
+ for(Int_t ii=0; ii<90; ii++){
+ countSSL->Add(new TVectorD(6));
+ }
+ }
+
+ //----------------------------------------
+ // Fill cobj0
+ //----------------------------------------
+
+ //ibin = binlowedge of bin(ibin+1) = the number fills this bin
+ for(Int_t ibin=0; ibin<nbin; ibin++){
+ Double_t gnor = hnor->GetBinContent(ibin+1);
+ Double_t gmpv = hmpv->GetBinContent(ibin+1);
+ Double_t gwid = hwid->GetBinContent(ibin+1);
+ Double_t gres = hres->GetBinContent(ibin+1);
+
+ //--- set additional cut by kpass
+ Bool_t kpass = kTRUE;
+ Double_t gtrdphmean = -999;
+ if(htrdphmean){
+ gtrdphmean = htrdphmean->GetBinContent(ibin+1);
+ //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
+ if(gtrdphmean<EPSILON){
+ kpass = kFALSE;
+ }
+ if(gnor<AliTRDdEdxBaseUtils::TimeBinCountCut()*gtrdphmean){
+ kpass = kFALSE;
+ }
+ }
+
+ //--- set calibration constant p0
+ Double_t p0= 0;
+
+ //reason for gmpv=0:
+ //1)gnor<=3; truncation in hist: (0, 0.6*ntot=1.8 with ntot=3]={1}, in hist entries can pile up so that ntot=2, or 3, and (ntot>nstart && ntot<=nstop) is skipped;
+ //2)TruncatedMean(hist) out of range (only for Q0, not Q1).
+
+ if(gmpv>EPSILON && kpass){
+ if(tag.Contains("T0")){
+ p0 = gmpv;
+ }
+ else{
+ p0 = 1/gmpv;
+ }
+ //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
+ }
+
+ (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0;
+
+ //--- save optional record
+ if(p0>EPSILON && countDet && countSSL){
+ const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
+ (*countDet)[idet]=1;
+
+ const Int_t isector = AliTRDdEdxBaseUtils::ToSector(ibin);
+ const Int_t istack = AliTRDdEdxBaseUtils::ToStack(ibin);
+ const Int_t ilayer = AliTRDdEdxBaseUtils::ToLayer(ibin);
+ TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector);
+ (*vecsectorstack)[ilayer]=1;
+ }
+
+ if(calibStream){
+ (*calibStream)<<tag<<
+ "run="<<run<<
+ "p0="<<p0<<
+
+ "nor="<<gnor<<
+ "mpv="<<gmpv<<
+ "wid="<<gwid<<
+ "res="<<gres<<
+ "gtrdphmean="<<gtrdphmean<<
+
+ "ibin="<<ibin<<
+ "\n";
+ }
+ }
+
+ //----------------------------------------
+ // Status Report
+ //----------------------------------------
+ if(countDet && countSSL){
+ TVectorD count2Dstack(90);
+ for(Int_t ii=0; ii<90; ii++){
+ TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
+ const Int_t nlayer = (Int_t)vecsectorstack->Sum();
+ if(nlayer==6){
+ count2Dstack[ii]=1;
+ }
+ }
+
+ printf("\nAliTRDdEdxCalibUtils::GetCalibObj Summary run: %d name: %s entries: %.0f ndetector: %03.0f n2dstack %02.0f\n\n", run, hname.Data(), hh->GetEntries(), countDet->Sum(), count2Dstack.Sum());
+ }
+
+ //----------------------------------------
+ // Clean Up
+ //----------------------------------------
+
+ TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
+ const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
+ for(Int_t ihh=0; ihh<nhh; ihh++){
+ if(!lout){
+ delete (*hhs[ihh]);
+ }
+ }
+
+ delete countDet;
+ delete countSSL;
+
+ //----------------------------------------
+
+ return cobj0;
+}
+
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+//
+//
+//
+// Xianguo Lu
+// lu@physi.uni-heidelberg.de
+// Xianguo.Lu@cern.ch
+//
+/*
+grep " AliTRDdEdxCalibUtils::" AliTRDdEdxCalibUtils.cxx | grep "=" -v | grep -v "[6]" | grep -v printf |wc
+grep "(" AliTRDdEdxCalibUtils.h | grep ";" | grep -v grep | grep -v ClassDef | grep -v "{" | grep -v typedef | wc
+*/
+
+
+#ifndef ALITRDDEDXCALIBUTILS_H
+#define ALITRDDEDXCALIBUTILS_H
+
+#ifndef TVECTORD_H
+#include "TVectorD.h"
+#endif
+
+#ifndef THNSPARSE_H
+#include "THnBase.h"
+#endif
+
+#ifndef TTREESTREAM_H
+#include "TTreeStream.h"
+#endif
+
+#ifndef ALITRDDEDXCALIBHISTARRAY_H
+#include "AliTRDdEdxCalibHistArray.h"
+#endif
+
+class TH1D;
+class TH2D;
+class TObjArray;
+
+class AliESDEvent;
+class AliESDtrack;
+class AliTRDcluster;
+class AliTRDtrackV1;
+class AliTRDseedV1;
+
+class AliTRDdEdxCalibUtils
+{
+ public:
+
+ static void SetObjArray(TObjArray * obj){fgObjArray = obj;}
+ static TObjArray * GetObjArray();
+ static TObjArray * GetObj(const Bool_t kinvq, const Double_t mag, const Int_t charge);
+ static TObjArray* HistToObj(const THnBase *hh, Int_t run=-999, TList *lout=0x0, TTreeSRedirector *calibStream=0x0);
+ static void DeleteObjArray();
+ static Bool_t GenerateDefaultOCDB(const TString path="local://./");
+
+ static AliTRDdEdxCalibHistArray * GetHistArray(){return fgHistArray;}
+ static THnBase * GetHistAt(const Int_t iter);
+ static void IniHistArray(TList *list, const Bool_t kNoInv);
+ static Bool_t ReadHistArray(const TString filename, const TString listname);
+ static void FillHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale) ;
+ static void DeleteHistArray();
+
+ static Double_t GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig);
+ static void Output(const TList *lin, Int_t run);
+
+ private:
+ static void FillHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnBase * hcalib, const Double_t scale);
+ static void GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean);
+
+
+ static AliTRDdEdxCalibHistArray * fgHistArray; //array containing 8 THnBase!
+ static TObjArray * fgObjArray; //array containing 8 TObjArray!
+
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//
+// TRD dEdx recon utils
+// xx
+// xx
+// xx
+// xx
+//
+// Xianguo Lu
+// lu@physi.uni-heidelberg.de
+// Xianguo.Lu@cern.ch
+//
+//
+
+#include "TF1.h"
+#include "TFile.h"
+#include "TH1D.h"
+#include "TH2D.h"
+#include "THnSparse.h"
+#include "TMath.h"
+#include "TMatrixD.h"
+#include "TMinuit.h"
+#include "TObjArray.h"
+#include "TRandom3.h"
+#include "TStopwatch.h"
+#include "TVectorD.h"
+
+#include "TTreeStream.h"
+
+#include "AliCDBId.h"
+#include "AliCDBMetaData.h"
+#include "AliCDBStorage.h"
+#include "AliESDEvent.h"
+#include "AliESDfriendTrack.h"
+#include "AliESDtrack.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDCalROC.h"
+#include "AliTRDtrackV1.h"
+
+#include "AliTRDdEdxBaseUtils.h"
+#include "AliTRDdEdxReconUtils.h"
+
+#define EPSILON 1e-12
+
+Int_t AliTRDdEdxReconUtils::ApplyCalib(const Int_t nc0, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
+{
+ //
+ //apply calibration on arrayQ
+ //
+ if(!cobj){ printf("AliTRDdEdxReconUtils::ApplyCalib error gain array null!!\n"); exit(1);}
+
+ TVectorD tmpq(arrayQ->GetNrows());
+ TVectorD tmpx(arrayX->GetNrows());
+ Int_t ncls = 0;
+
+ const TVectorD * gain = (TVectorD*) cobj->At(0);
+ for(Int_t ii=0; ii<nc0; ii++){
+ const Double_t dq = (*arrayQ)[ii];
+ const Int_t xx = (Int_t)(*arrayX)[ii];
+ const Double_t gg = (*gain)[xx];
+
+ if(gg<EPSILON){
+ continue;
+ }
+
+ tmpq[ncls] = dq*gg;
+ tmpx[ncls] = xx;
+ ncls++;
+ }
+
+ (*arrayQ)=tmpq;
+ (*arrayX)=tmpx;
+
+ return ncls;
+}
+
+Double_t AliTRDdEdxReconUtils::ToyCook(const Bool_t kinvq, Int_t &ncluster, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
+{
+ //
+ //template for cookdedx
+ //
+ if(cobj){
+ if(arrayQ && arrayX){
+ ncluster = ApplyCalib(ncluster, arrayQ, arrayX, cobj);
+ }
+ else{
+ printf("AliTRDdEdxReconUtils::ToyCook arrayQ arrayX null, applycalib can not be applied!\n"); exit(1);
+ }
+ }
+
+ Double_t lowFrac =-999, highFrac = -999;
+ if(kinvq){
+ lowFrac = AliTRDdEdxBaseUtils::Q1Frac(); highFrac = 0.99;
+ }
+ else{
+ lowFrac = 0.01; highFrac = AliTRDdEdxBaseUtils::Q0Frac();
+ }
+
+ Double_t meanQ = AliTRDdEdxBaseUtils::TruncatedMean(ncluster, arrayQ->GetMatrixArray(), lowFrac, highFrac);
+ if(kinvq){
+ if(meanQ>EPSILON){
+ meanQ = 1/meanQ;
+ }
+ }
+
+ return meanQ;
+}
+
+Double_t AliTRDdEdxReconUtils::CombineddEdx(const Bool_t kinvq, Int_t &concls, TVectorD *coarrayQ, TVectorD *coarrayX, const Int_t tpcncls, const TVectorD *tpcarrayQ, const TVectorD *tpcarrayX, const Int_t trdncls, const TVectorD *trdarrayQ, const TVectorD *trdarrayX)
+{
+ //
+ //combine tpc and trd dedx
+ //
+
+ for(Int_t iq=0; iq<tpcncls; iq++){
+ (*coarrayQ)[iq]=(*tpcarrayQ)[iq];
+ if(tpcarrayX && trdarrayX && coarrayX){
+ (*coarrayX)[iq]=(*tpcarrayX)[iq];
+ }
+ }
+ for(Int_t iq=0; iq<trdncls; iq++){
+ (*coarrayQ)[tpcncls+iq]=(*trdarrayQ)[iq];
+ if(tpcarrayX && trdarrayX && coarrayX){
+ (*coarrayX)[tpcncls+iq]=159+(*trdarrayX)[iq];
+ }
+ }
+
+ concls=trdncls+tpcncls;
+
+ const Double_t coQ = ToyCook(kinvq, concls, coarrayQ, coarrayX);
+
+ return coQ;
+}
+
+Double_t AliTRDdEdxReconUtils::GetAngularCorrection(const AliTRDseedV1 *seed)
+{
+ //
+ //return angular normalization factor
+ //
+
+ return TMath::Sqrt(1+seed->GetYref(1)*seed->GetYref(1)+seed->GetZref(1)*seed->GetZref(1));
+}
+
+Double_t AliTRDdEdxReconUtils::GetPadGain(const Int_t det, const Int_t icol, const Int_t irow)
+{
+ //
+ //get pad calibration
+ //
+ AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
+ if (!calibration) {
+ printf("AliTRDdEdxReconUtils::GetPadCalib No AliTRDcalibDB instance available\n"); exit(1);
+ }
+ AliTRDCalROC * calGainFactorROC = calibration->GetGainFactorROC(det);
+ if(!calGainFactorROC){
+ printf("AliTRDdEdxReconUtils::GetPadCalib no calGainFactorROC!\n"); exit(1);
+ }
+
+ Double_t padgain = -999;
+ if( icol >= 0 &&
+ icol < calGainFactorROC->GetNcols() &&
+ irow >=0 &&
+ irow < calGainFactorROC->GetNrows()){
+ padgain = calGainFactorROC->GetValue(icol, irow);
+ if(padgain<EPSILON){
+ printf("AliTRDdEdxReconUtils::GetPadGain padgain 0! %f %f -- %d %d %d -- %d %d\n", padgain, EPSILON, det, icol, irow, calGainFactorROC->GetNcols(), calGainFactorROC->GetNrows()); exit(1);
+ }
+ }
+ else{
+ //printf("\nAliTRDdEdxReconUtils::GetPadGain warning!! indices out of range %d %d %d -- %d %d\n\n", det, icol, irow, calGainFactorROC->GetNcols(), calGainFactorROC->GetNrows() );
+ }
+
+ return padgain;
+}
+
+Double_t AliTRDdEdxReconUtils::GetRNDClusterQ(AliTRDcluster *cl)
+{
+ //
+ //get cluter q from GetRawQ, apply baseline and Kr pad-calibration
+ //
+
+ const Int_t det = cl->GetDetector();
+ const Int_t pad3col = cl->GetPadCol();
+ const Int_t padrow = cl->GetPadRow();
+
+ const Double_t baseline = 10;
+
+ Double_t rndqsum = 0;
+ for(Int_t ii=0; ii<7; ii++){
+ if(cl->GetSignals()[ii] < EPSILON){//bad pad marked by electronics
+ continue;
+ }
+
+ const Int_t icol = pad3col+(ii-3);
+ const Double_t padgain = GetPadGain(det, icol, padrow);
+ if(padgain<0){//indices out of range, pad3col near boundary case
+ continue;
+ }
+
+ const Double_t rndsignal = (cl->GetSignals()[ii] - baseline)/(AliTRDdEdxBaseUtils::IsPadGainOn()? padgain : 1);
+
+ //sum it anyway even if signal below baseline, as long as the total is positive
+ rndqsum += rndsignal;
+ }
+
+ return rndqsum;
+}
+
+Double_t AliTRDdEdxReconUtils::GetClusterQ(const Bool_t kinvq, const AliTRDseedV1 * seed, const Int_t itb)
+{
+ //
+ //get cluster charge
+ //
+ Double_t dq = 0;
+ AliTRDcluster *cl = 0x0;
+
+ //GetRNDClusterQ(cl)>0 ensures that the total sum of q is above baseline*NsignalPhysical.
+ cl = seed->GetClusters(itb); if(cl && GetRNDClusterQ(cl)>0 ) dq+= GetRNDClusterQ(cl);//cl->GetRawQ();
+ cl = seed->GetClusters(itb+AliTRDseedV1::kNtb); if(cl && GetRNDClusterQ(cl)>0 ) dq+= GetRNDClusterQ(cl);//cl->GetRawQ();
+
+ dq /= GetAngularCorrection(seed);
+
+ dq /= AliTRDdEdxBaseUtils::QScale();
+
+ if(kinvq){
+ if(dq>EPSILON){
+ dq = 1/dq;
+ }
+ }
+
+ return dq;
+}
+
+Int_t AliTRDdEdxReconUtils::GetArrayClusterQ(const Bool_t kinvq, TVectorD *arrayQ, TVectorD *arrayX, const AliTRDtrackV1 *trdtrack, Int_t timeBin0, Int_t timeBin1, Int_t tstep)
+{
+ //
+ //return nclustter
+ //(if kinvq, return 1/q array), size of array must be larger than 31*6
+ //
+ if(!arrayQ || arrayQ->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
+ printf("AliTRDdEdxReconUtils::GetArrayClusterQ error arrayQ null or size too small! %d\n", arrayQ? arrayQ->GetNrows() : -999); exit(1);
+ }
+ if(!arrayX || arrayX->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
+ printf("AliTRDdEdxReconUtils::GetArrayClusterQ error arrayX null or size too small! %d\n", arrayX? arrayX->GetNrows() : -999); exit(1);
+ }
+
+ const Int_t mintb = 0;
+ const Int_t maxtb = AliTRDseedV1::kNtb-1;
+ if(timeBin0<mintb) timeBin0=mintb;
+ if(timeBin1>maxtb) timeBin1=maxtb;
+ if(tstep<=0) tstep=1;
+
+ //============
+ Int_t tbN=0;
+ Double_t tbQ[200];
+ Int_t tbBin[200];
+
+ for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
+ const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
+ if(!seed)
+ continue;
+
+ const Int_t det = seed->GetDetector();
+
+ for(Int_t itb=timeBin0; itb<=timeBin1; itb+=tstep){
+ const Double_t dq = GetClusterQ(kinvq, seed, itb);
+ if(dq<EPSILON)
+ continue;
+
+ const Int_t gtb = det * AliTRDseedV1::kNtb + itb;
+
+ tbQ[tbN]=dq;
+ tbBin[tbN]=gtb;
+ tbN++;
+ }
+ }
+
+ Int_t ncls = 0;
+ for(Int_t iq=0; iq<tbN; iq++){
+ if(tbQ[iq]<EPSILON)
+ continue;
+
+ (*arrayQ)[ncls] = tbQ[iq];
+ (*arrayX)[ncls] = tbBin[iq];
+
+ ncls++;
+ }
+
+ static Int_t kprint = 100;
+ if(kprint<0){
+ printf("\nAliTRDdEdxReconUtils::GetArrayClusterQ raw cluster-Q\n");
+ for(Int_t iq=0; iq<ncls; iq++){
+ const Int_t ichamber = AliTRDdEdxBaseUtils::ToLayer((*arrayX)[iq]);
+ const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
+ if(!seed){
+ printf("error seed null!!\n"); exit(1);
+ }
+ const Double_t rawq = (*arrayQ)[iq] * 45. * GetAngularCorrection(seed);
+ printf("esdid=%d; chamber=%d; timebin=%d; rawq= %.3f; myq[%d]= %e;\n", trdtrack->GetESDid(), ichamber, AliTRDdEdxBaseUtils::ToTimeBin((*arrayX)[iq]), rawq, iq, (*arrayQ)[iq]);
+ }
+ printf("\n");
+ }
+ kprint++;
+
+ return ncls;
+}
+
+Int_t AliTRDdEdxReconUtils::UpdateArrayX(const Int_t ncls, TVectorD* arrayX)
+{
+ //
+ //arrayX det*Ntb+itb -> itb
+ //
+
+ TVectorD countChamber(6);
+ for(Int_t ii = 0; ii<ncls; ii++){
+ const Int_t xx = (Int_t)(*arrayX)[ii];
+ const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(xx);
+
+ const Double_t ich = AliTRDgeometry::GetLayer(idet);
+ const Double_t itb = AliTRDdEdxBaseUtils::ToTimeBin(xx);
+ (*arrayX)[ii] = ich*AliTRDseedV1::kNtb+itb;
+
+ countChamber[ich] = 1;
+ }
+
+ const Double_t nch = countChamber.Sum();
+ return (Int_t) nch;
+}
+
+
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+//
+//
+//
+// Xianguo Lu
+// lu@physi.uni-heidelberg.de
+// Xianguo.Lu@cern.ch
+//
+/*
+grep " AliTRDdEdxReconUtils::" AliTRDdEdxReconUtils.cxx | grep "=" -v | grep -v "[6]" | grep -v printf |wc
+grep "(" AliTRDdEdxReconUtils.h | grep ";" | grep -v grep | grep -v ClassDef | grep -v "{" | grep -v typedef | wc
+*/
+
+
+#ifndef ALITRDDEDXRECONUTILS_H
+#define ALITRDDEDXRECONUTILS_H
+
+#ifndef TVECTORD_H
+#include "TVectorD.h"
+#endif
+
+#ifndef THNSPARSE_H
+#include "THnBase.h"
+#endif
+
+#ifndef TTREESTREAM_H
+#include "TTreeStream.h"
+#endif
+
+class TH1D;
+class TH2D;
+class TObjArray;
+
+class AliESDEvent;
+class AliESDtrack;
+class AliTRDcluster;
+class AliTRDtrackV1;
+class AliTRDseedV1;
+
+class AliTRDdEdxReconUtils
+{
+ public:
+
+
+ static Int_t ApplyCalib(const Int_t nc0, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj);
+ static Double_t ToyCook(const Bool_t kinvq, Int_t &ncluster, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj=0x0);
+ static Double_t CombineddEdx(const Bool_t kinvq, Int_t &concls, TVectorD *coarrayQ, TVectorD *coarrayX, const Int_t tpcncls, const TVectorD *tpcarrayQ, const TVectorD *tpcarrayX, const Int_t trdncls, const TVectorD *trdarrayQ, const TVectorD *trdarrayX);
+ static Int_t GetArrayClusterQ(const Bool_t kinvq, TVectorD *arrayQ, TVectorD *arrayX, const AliTRDtrackV1 *trdtrack, Int_t timeBin0=-1, Int_t timeBin1=1000, Int_t tstep=1);
+ static Int_t UpdateArrayX(const Int_t ncls, TVectorD* arrayX);
+
+ private:
+ static Double_t GetAngularCorrection(const AliTRDseedV1 *seed);
+ static Double_t GetPadGain(const Int_t det, const Int_t icol, const Int_t irow);
+ static Double_t GetRNDClusterQ(AliTRDcluster *cl);
+ static Double_t GetClusterQ(const Bool_t kinvq, const AliTRDseedV1 * seed, const Int_t itb);
+
+};
+
+#endif
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-
-//
-// class to calculate TRD dEdx
-// xx
-// xx
-// xx
-// xx
-//
-// Xianguo Lu
-// lu@physi.uni-heidelberg.de
-// Xianguo.Lu@cern.ch
-//
-//
-
-
-#include "TF1.h"
-#include "TFile.h"
-#include "TH1D.h"
-#include "TH2D.h"
-#include "THn.h"
-#include "THnSparse.h"
-#include "TMath.h"
-#include "TMatrixD.h"
-#include "TMinuit.h"
-#include "TObjArray.h"
-#include "TRandom3.h"
-#include "TStopwatch.h"
-#include "TVectorD.h"
-
-#include "TTreeStream.h"
-
-#include "AliCDBId.h"
-#include "AliCDBMetaData.h"
-#include "AliCDBStorage.h"
-#include "AliESDEvent.h"
-#include "AliESDfriendTrack.h"
-#include "AliESDtrack.h"
-#include "AliTRDcalibDB.h"
-#include "AliTRDCalROC.h"
-#include "AliTRDtrackV1.h"
-
-
-#include "AliTRDdEdxUtils.h"
-
-#define EPSILON 1e-12
-
-THnF * AliTRDdEdxUtils::fgHistGain=0x0;
-THnF * AliTRDdEdxUtils::fgHistT0=0x0;
-THnF * AliTRDdEdxUtils::fgHistVd=0x0;
-TObjArray * AliTRDdEdxUtils::fgHistPHQ=new TObjArray(8);
-
-TString AliTRDdEdxUtils::fgCalibFile;
-TObjArray * AliTRDdEdxUtils::fgObjGain = 0x0;
-TObjArray * AliTRDdEdxUtils::fgObjT0 = 0x0;
-TObjArray * AliTRDdEdxUtils::fgObjVd = 0x0;
-TObjArray * AliTRDdEdxUtils::fgObjPHQ = 0x0;
-
-Int_t AliTRDdEdxUtils::fgNchamber = -999;
-Double_t AliTRDdEdxUtils::fgChamberQ[6];
-Double_t AliTRDdEdxUtils::fgChamberTmean[6];
-
-Double_t AliTRDdEdxUtils::fgTrackTmean = -999;
-
-Bool_t AliTRDdEdxUtils::fgPadGainOn = kTRUE;
-Bool_t AliTRDdEdxUtils::fgExBOn = kTRUE;
-Double_t AliTRDdEdxUtils::fgQScale = 45;
-Double_t AliTRDdEdxUtils::fgQ0Frac = 0.3;
-Double_t AliTRDdEdxUtils::fgQ1Frac = 0.5;
-Double_t AliTRDdEdxUtils::fgTimeBinCountCut = 0.0;
-Int_t AliTRDdEdxUtils::fgCalibTPCnclsCut = 70;
-
-//===================================================================================
-// Math and Histogram
-//===================================================================================
-void AliTRDdEdxUtils::GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres)
-{
- //
- //counts of hh is sorted
- //
-
- for(Int_t ii=0; ii<ncut; ii++){
- cuts[ii] = -999;
- }
-
- Int_t nsel = 0;
- const Int_t nbin = hh->GetNbinsX();
- Double_t datas[nbin];
- for(Int_t ii=1; ii<=nbin; ii++){
- const Double_t res = hh->GetBinContent(ii);
- if(res<thres){
- continue;
- }
-
- datas[nsel] = res;
- nsel++;
- }
- if(!nsel)
- return;
-
- Int_t id[nsel];
- TMath::Sort(nsel, datas, id, kFALSE);
-
- for(Int_t ii=0; ii<ncut; ii++){
- const Double_t icdf = cdfs[ii];
- if(icdf<0 || icdf>1){
- printf("AliTRDdEdxUtils::GetCDFCuts error cdfs[%d] %15f out of range!\n", ii, icdf); exit(1);
- }
- cuts[ii] = datas[id[Int_t(icdf*nsel)]];
- }
-}
-
-Double_t AliTRDdEdxUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr)
-{
- //
- //calculate mean (with error) and rms from sum, w2s, nn
- //if nn=0, mean, error, and rms all = 0
- //
-
- Double_t tmean = 0, trms = 0, terr = 0;
-
- if(nn>EPSILON){
- tmean = sum/nn;
-
- const Double_t arg = w2s/nn-tmean*tmean;
- if(TMath::Abs(arg)<EPSILON){
- trms = 0;
- }
- else{
- if( arg <0 ){
- printf("AliTRDdEdxUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1);
- }
-
- trms = TMath::Sqrt(arg);
- }
-
- terr = trms/TMath::Sqrt(nn);
- }
-
- if(grms){
- (*grms) = trms;
- }
-
- if(gerr){
- (*gerr) = terr;
- }
-
- return tmean;
-}
-
-Double_t AliTRDdEdxUtils::TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr, Double_t *wws)
-{
- //
- //calculate truncated mean
- //return <x*w>_{low-high according to x}
- //
-
- /*
- //test->
- for(Int_t ii=0; ii<nx; ii++){
- printf("test %d/%d %f\n", ii, nx, xdata[ii]);
- }
- //<--test
- */
-
- Int_t index[nx];
- TMath::Sort(nx, xdata, index, kFALSE);
-
- Int_t nused = 0;
- Double_t sum = 0;
- Double_t w2s = 0;
- const Int_t istart = Int_t (nx*lowfrac);
- const Int_t istop = Int_t (nx*highfrac);
-
- //=,< correct, because when low=0, high=1 it is correct
- for(Int_t ii=istart; ii<istop; ii++){
- Double_t weight = 1;
- if(wws){
- weight = wws[index[ii]];
- }
- const Double_t sx = xdata[index[ii]]*weight;
-
- sum += sx;
- w2s += sx*sx;
-
- nused++;
- //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s);
-
- }
-
- return GetMeanRMS(nused, sum, w2s, grms, gerr);
-}
-
-Double_t AliTRDdEdxUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr)
-{
- //
- //do truncation on histogram
- //
- //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct
-
- //with under-/over-flow
- Double_t npreTrunc = 0;
- for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
- const Double_t be = hh->GetBinError(itmp);
- const Double_t bc = hh->GetBinContent(itmp);
- if(be<EPSILON){
- if(bc>EPSILON){
- printf("AliTRDdEdxUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
- }
- continue;
- }
- npreTrunc += bc*bc/be/be;
- }
-
- const Double_t nstart = npreTrunc*lowfrac;
- const Double_t nstop = npreTrunc*highfrac;
-
- //with Double_t this should also handle normalized hist
- Double_t ntot = 0;
- Double_t nused = 0;
- Double_t sum = 0;
- Double_t w2s = 0;
- for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
- const Double_t be = hh->GetBinError(itmp);
- const Double_t bc = hh->GetBinContent(itmp);
- if(be<EPSILON){
- if(bc>EPSILON){
- printf("AliTRDdEdxUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
- }
- continue;
- }
- const Double_t weight = bc*bc/be/be;
- ntot+=weight;
- //<= correct, because when high=1, nstop has to be included
- if(ntot>nstart && ntot<=nstop){
-
- const Double_t val = hh->GetBinCenter(itmp);
- sum += weight*val;
- w2s += weight*val*val;
-
- nused += weight;
-
- //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample);
- }
- else if(ntot>nstop){
- if(itmp>=hh->GetNbinsX()){
- printf("AliTRDdEdxUtils::TruncatedMean warning hist range too small %s %f %f %d %d, %15f %15f %15f; nused w2s sum set to 0\n", hh->GetName(), hh->GetBinLowEdge(1), hh->GetBinLowEdge(itmp), itmp, hh->GetNbinsX(), hh->GetBinContent(hh->GetNbinsX())/hh->Integral(0,hh->GetNbinsX()+1), hh->GetBinContent(hh->GetNbinsX()), hh->Integral(0,hh->GetNbinsX()+1)); //exit(1);
- nused = 0;
- w2s = sum = 0;
- }
- break;
- }
- }
-
- return GetMeanRMS(nused, sum, w2s, grms, gerr);
-}
-
-void AliTRDdEdxUtils::FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, const Double_t thres, const Double_t lowfrac, const Double_t highfrac)
-{
- //
- //fit slices of hh using truncation
- //
-
- const Int_t x0 = hh->GetXaxis()->GetFirst();
- const Int_t x1 = hh->GetXaxis()->GetLast();
- const Int_t y0 = hh->GetYaxis()->GetFirst();
- const Int_t y1 = hh->GetYaxis()->GetLast();
-
- const Int_t nx = hh->GetNbinsX();
- const Int_t ny = hh->GetNbinsY();
- const Double_t xmin = hh->GetXaxis()->GetXmin();
- const Double_t xmax = hh->GetXaxis()->GetXmax();
- const Double_t ymin = hh->GetYaxis()->GetXmin();
- const Double_t ymax = hh->GetYaxis()->GetXmax();
-
- hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax);
- hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax);
- hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax);
- hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax);
-
- for(Int_t ix=x0; ix<=x1; ix++){
- //to speed up
- const Double_t rawcount = hh->Integral(ix,ix,0, ny+1);
- if(rawcount<EPSILON){
- continue;
- }
-
- TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
- Double_t ntot = 0;
- for(Int_t iy=y0; iy<=y1; iy++){
- const Double_t be = hh->GetBinError(ix,iy);
- const Double_t bc = hh->GetBinContent(ix, iy);
-
- if(be<EPSILON){
- if(bc>EPSILON){
- printf("AliTRDdEdxUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
- }
- continue;
- }
-
- htmp->SetBinContent(iy, bc);
- htmp->SetBinError(iy, be);
-
- ntot += (bc/be)*(bc/be);
-
- //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2));
- }
-
- hnor->SetBinContent(ix, ntot);
- hnor->SetBinError( ix, 0);
-
- if(ntot<thres || htmp->GetRMS()<EPSILON){
- delete htmp;
- continue;
- }
-
- //test htmp->Draw();
- Double_t trms = -999, terr = -999;
- const Double_t tmean = TruncatedMean(htmp, lowfrac, highfrac, &trms, &terr);
-
- hmpv->SetBinContent(ix, tmean);
- hmpv->SetBinError( ix, terr);
-
- hwid->SetBinContent(ix, trms);
- hwid->SetBinError( ix, 0);
-
- hres->SetBinContent(ix, tmean>EPSILON ? trms/tmean:0);
- hres->SetBinError( ix, 0);
-
- delete htmp;
- }
-
- TH1 *hhs[]={hnor, hmpv, hwid, hres};
- const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV"};
- const Int_t nh = sizeof(hhs)/sizeof(TH1*);
- for(Int_t ii=0; ii<nh; ii++){
- hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle()));
- hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle());
- hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset());
- hhs[ii]->SetTitle(hh->GetTitle());
- }
-}
-
-//===================================================================================
-// TRD Analysis Fast Tool
-//===================================================================================
-
-Int_t AliTRDdEdxUtils::GetNtracklet(const AliESDEvent *esd)
-{
- //
- //number of trd tracklet in one esd event
- //
- const Int_t ntrk0 = esd->GetNumberOfTracks();
- Int_t ntrdv1=0;
- for(Int_t ii=0; ii<ntrk0; ii++){
- ntrdv1 += esd->GetTrack(ii)->GetTRDntracklets();
- }
- return ntrdv1;
-}
-
-AliTRDtrackV1 * AliTRDdEdxUtils::GetTRDtrackV1(const AliESDtrack * esdtrack)
-{
- //
- //Get TRD friend track
- //
-
- AliESDfriendTrack * friendtrk = (AliESDfriendTrack *)esdtrack->GetFriendTrack();
- if(!friendtrk){
- //printf("xlulog AliAnalysisTaskCosmicTRD::GetTRDtrack no friend!!\n"); exit(1);
- return 0x0;
- }
-
- TObject *calibObject=0x0;
- AliTRDtrackV1 * trdtrack=0x0;
- for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) {
- if( (trdtrack=dynamic_cast<AliTRDtrackV1*>(calibObject)) )
- break;
- }
-
- return trdtrack;
-}
-
-Bool_t AliTRDdEdxUtils::IsInSameStack(const AliTRDtrackV1 *trdtrack)
-{
- //
- // to check if all tracklets are in the same stack, useful in cosmic
- //
-
- TVectorD secs(18), stks(5);
-
- for(Int_t ilayer = 0; ilayer < 6; ilayer++){
- AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
- if(!tracklet)
- continue;
-
- const Int_t det = tracklet->GetDetector();
- const Int_t isector = AliTRDgeometry::GetSector(det);
- const Int_t istack = AliTRDgeometry::GetStack(det);
-
- secs[isector] = 1;
- stks[istack] = 1;
- }
-
- if(secs.Sum()!=1 || stks.Sum()!=1){
- return kFALSE;
- }
- else
- return kTRUE;
-}
-
-Bool_t AliTRDdEdxUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
-{
- //
- //as function name
- //
- isec = istk = -999;
- mom = -999;
-
- for(Int_t ilayer = 0; ilayer < 6; ilayer++){
- AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
- if(!tracklet)
- continue;
-
- const Int_t det = tracklet->GetDetector();
- isec = AliTRDgeometry::GetSector(det);
- istk = AliTRDgeometry::GetStack(det);
-
- mom = tracklet->GetMomentum();
-
- break;
- }
-
- if(isec<0)
- return kFALSE;
- else
- return kTRUE;
-}
-
-//===================================================================================
-// Calibration
-//===================================================================================
-Double_t AliTRDdEdxUtils::GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig)
-{
- //
- //the scale used in calibration
- //
-
- if(tpcncls < CalibTPCnclsCut())
- return -999;
-
- if(tpcsig<EPSILON)
- return -999;
-
- return tpcsig/120;
-
-}
-
-Int_t AliTRDdEdxUtils::GetPHQIterator(const Bool_t kinvq, const Double_t mag, const Int_t charge)
-{
- //
- //iterator for calib obj and hist
- //
- return kinvq*4 + (mag>0)*2 + (charge>0);
-}
-
-TObjArray * AliTRDdEdxUtils::GetObjPHQ()
-{
- //
- //return fgObjPHQ, initialized if null
- //
-
- if(!fgObjPHQ){
- fgObjPHQ = new TObjArray(8);
- }
-
- return fgObjPHQ;
-}
-
-TObjArray * AliTRDdEdxUtils::GetObjPHQ(const Bool_t kinvq, const Double_t mag, const Int_t charge)
-{
- //
- //return calib obj
- //
- if(!fgObjPHQ){
- printf("AliTRDdEdxUtils::GetObjPHQ(kinvq, mag, charge) error fgObjPHQ null!!\n"); exit(1);
- }
-
- return (TObjArray*) fgObjPHQ->At(GetPHQIterator(kinvq, mag, charge));
-}
-
-THnF * AliTRDdEdxUtils::GetHistPHQ(const Bool_t kinvq, const Double_t mag, const Int_t charge)
-{
- //
- //return calib hist
- //
- return (THnF*) fgHistPHQ->At(GetPHQIterator(kinvq, mag, charge));
-}
-
-TString AliTRDdEdxUtils::GetPHQName(const Bool_t kobj, const Int_t iter)
-{
- //
- //get name of calib obj/hist of PHQ
- //
- return Form("TRDCalib%sPHQ%d", kobj?"Obj":"Hist", iter);
-}
-
-void AliTRDdEdxUtils::DeleteCalibObj()
-{
- //
- //delete calib obj
- //
- delete fgObjGain;
- delete fgObjT0;
- delete fgObjVd;
-
- fgObjGain = 0x0;
- fgObjT0 = 0x0;
- fgObjVd = 0x0;
-
- if(fgObjPHQ){
- fgObjPHQ->SetOwner();
- delete fgObjPHQ;
- fgObjPHQ = 0x0;
- }
-}
-
-Bool_t AliTRDdEdxUtils::GenerateDefaultPHQOCDB(const TString path)
-{
- //
- //generate default OCDB object PHQ, do like
- //AliTRDdEdxUtils::GenerateDefaultPHQOCDB("local://./")
- //
-
- TObjArray * arr8 = new TObjArray(8);
- arr8->SetOwner();
-
- for(Int_t ii=0; ii<8; ii++){
- TObjArray * arr1 = new TObjArray(1);
- arr1->SetOwner();
- arr1->SetName(GetPHQName(1, ii));
-
- const Int_t nbins = NTRDtimebin();
- TVectorD * vec = new TVectorD(nbins);
- for(Int_t jj=0; jj<nbins; jj++){
- (*vec)[jj] = 1;
- }
- arr1->Add(vec);
- arr8->Add(arr1);
- }
-
- AliCDBMetaData *metaData= new AliCDBMetaData();
- metaData->SetObjectClassName("TObjArray");
- metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu");
-
- AliCDBId id1("TRD/Calib/PHQ", 0, 999999999);
- AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage(path);
- gStorage->Put(arr8, id1, metaData);
-
- delete metaData;
- delete arr8;
-
- return kTRUE;
-}
-
-void AliTRDdEdxUtils::IniCalibObj()
-{
- //
- //set CalibObj from file, clone to static calib obj
- //
-
- DeleteCalibObj();
-
- TFile *cfile=new TFile(fgCalibFile);
- if(!cfile->IsOpen()){
- printf("AliTRDdEdxUtils::IniCalibObj error fgCalibFile not open! %s\n", fgCalibFile.Data());exit(1);
- }
-
- printf("\nAliTRDdEdxUtils::IniCalibObj file: %s\n", fgCalibFile.Data());
-
- //---
- const TString objnames[] ={"TRDCalibObjGain", "TRDCalibObjT0", "TRDCalibObjVd"};
- TObjArray ** gobjs[]={ &fgObjGain, &fgObjT0, &fgObjVd};
-
- const Int_t nobj = sizeof(objnames)/sizeof(TString);
- for(Int_t iobj=0; iobj<nobj; iobj++){
- TObjArray *tmpo=0x0;
- cfile->GetObject(objnames[iobj], tmpo);
- if(!tmpo){
- printf("AliTRDdEdxUtils::IniCalibObj error obj %s not found!\n", objnames[iobj].Data()); exit(1);
- }
-
- (*gobjs[iobj])=(TObjArray*)tmpo->Clone();
- (*gobjs[iobj])->SetOwner();
- }
-
- fgObjPHQ = new TObjArray(8);
- for(Int_t iter=0; iter<8; iter++){
- const TString objn = GetPHQName(1, iter);
- TObjArray *tmpo=0x0;
- cfile->GetObject(objn, tmpo);
- if(!tmpo){
- printf("AliTRDdEdxUtils::IniCalibObj error obj %s not found!\n", objn.Data()); exit(1);
- }
-
- TObjArray *obji=(TObjArray*) tmpo->Clone();
- obji->SetOwner();
- fgObjPHQ->AddAt(obji, iter);
- }
-
- //---
-
- cfile->Close();
- delete cfile;
-}
-
-void AliTRDdEdxUtils::DeleteCalibHist()
-{
- //
- //delete calib hist
- //
- delete fgHistGain;
- delete fgHistT0;
- delete fgHistVd;
-
- fgHistGain = 0x0;
- fgHistT0 = 0x0;
- fgHistVd = 0x0;
-
- //fgHistPHQ owns the hists
- fgHistPHQ->SetOwner();
- fgHistPHQ->Clear();
-}
-
-void AliTRDdEdxUtils::IniCalibHist(TList *list, const Bool_t kPHQonly)
-{
- //
- //initialize calib hist, list should not own the hist, or list->Clear/delete hist should not be called
- //
-
- DeleteCalibHist();
-
- Int_t nbin[2];
- const Double_t xmin[2]={0, 0};
- Double_t xmax[2];
-
- nbin[0]=NTRDtimebin(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=20;
- for(Int_t iter=0; iter<8; iter++){
- const TString hn = GetPHQName(0, iter);
- THnF *hi = new THnF(hn, "", 2, nbin, xmin, xmax);
- //fgHistPHQ owns the hists
- fgHistPHQ->AddAt(hi, iter);
- list->Add(hi);
- }
-
- if(kPHQonly)
- return;
-
- nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=20; fgHistGain = new THnF("TRDCalibHistGain", "", 2, nbin, xmin, xmax);
- nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=AliTRDseedV1::kNtb; fgHistT0 = new THnF("TRDCalibHistT0", "", 2, nbin, xmin, xmax);
- nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=AliTRDseedV1::kNtb; fgHistVd = new THnF("TRDCalibHistVd", "", 2, nbin, xmin, xmax);
-
- list->Add(fgHistGain);
- list->Add(fgHistT0);
- list->Add(fgHistVd);
-}
-
-Bool_t AliTRDdEdxUtils::ReadCalibHist(const TString filename, const TString listname)
-{
- //
- //used in AliTRDPreprocessorOffline
- //read in calib hist from file, only for PHQ
- //
- DeleteCalibHist();
-
- //maybe already open by others... don't close
- TFile fcalib(filename);
-
- TObjArray * array = (TObjArray*)fcalib.Get(listname);
-
- for(Int_t iter=0; iter<8; iter++){
- const TString hn = GetPHQName(0, iter);
- THnF * tmph=0x0;
- if(array){
- tmph = (THnF *) array->FindObject(hn);
- }
- else{
- tmph = (THnF *) fcalib.Get(hn);
- }
- if(!tmph){
- printf("AliTRDdEdxUtils::ReadCalibHist warning calib hist not found! %s %s\n", filename.Data(), listname.Data());
- fcalib.ls();
- if(array){
- array->ls();
- }
- return kFALSE;
- }
- THnF *hi = (THnF*)tmph->Clone();
- fgHistPHQ->AddAt(hi, iter);
- }
-
- return kTRUE;
-}
-
-void AliTRDdEdxUtils::FillCalibHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnF * hcalib, const Double_t scale)
-{
- //
- //fill calibration hist
- //
- if(!hcalib){printf("AliTRDdEdxUtils::FillCalibHist errro hcalib null!!\n"); exit(1);}
-
- for(Int_t ii=0; ii<ncls; ii++){
- const Double_t dq = (*arrayQ)[ii];
- const Double_t xx = (*arrayX)[ii];
-
- const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
- const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
- const Double_t xmax = hcalib->GetAxis(0)->GetXmax();
-
- if(xx>=xmax || xx<xmin){
- printf("AliTRDdEdxUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(), xx, xmin, xmax); exit(1);
- }
-
- const Double_t var[]={xx, TMath::Min(dq, qmax)/scale};
- hcalib->Fill(var);
- }
-}
-
-void AliTRDdEdxUtils::FillCalibHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale)
-{
- //
- //get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled
- //
-
- THnF * hcalib = AliTRDdEdxUtils::GetHistPHQ(kinvq, mag, charge);
-
- TVectorD arrayQ(200), arrayX(200);
- const Int_t ncls = AliTRDdEdxUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
- FillCalibHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);
-
- static Int_t kprint = 100;
- if(kprint<0){
- printf("\nAliTRDdEdxUtils::FillCalibHist summary: \n");
- printf("\nkinvq= %d;\n", kinvq);
- for(Int_t iq=0; iq<ncls; iq++){
- printf("arrayX[%3d] = %15.0f; arrayQ[%3d] = %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
- }
- printf("\n");
- }
- kprint++;
-}
-
-Int_t AliTRDdEdxUtils::ApplyCalib(const Int_t nc0, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
-{
- //
- //apply calibration on arrayQ
- //
- if(!cobj){ printf("AliTRDdEdxUtils::ApplyCalib error gain array null!!\n"); exit(1);}
-
- TVectorD tmpq(arrayQ->GetNrows());
- TVectorD tmpx(arrayX->GetNrows());
- Int_t ncls = 0;
-
- const TVectorD * gain = (TVectorD*) cobj->At(0);
- for(Int_t ii=0; ii<nc0; ii++){
- const Double_t dq = (*arrayQ)[ii];
- const Int_t xx = (Int_t)(*arrayX)[ii];
- const Double_t gg = (*gain)[xx];
-
- if(gg<EPSILON){
- continue;
- }
-
- tmpq[ncls] = dq*gg;
- tmpx[ncls] = xx;
- ncls++;
- }
-
- (*arrayQ)=tmpq;
- (*arrayX)=tmpx;
-
- return ncls;
-}
-
-void AliTRDdEdxUtils::GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean)
-{
- //
- //calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output
- //
- const Int_t ndet = 540;
- TObjArray *obj=new TObjArray(ndet);
- obj->SetOwner();
- for(Int_t ii=0; ii<ndet; ii++){
- obj->Add(new TVectorD(AliTRDseedV1::kNtb));
- }
-
- //ibin = binlowedge of bin(ibin+1) = the number fills this bin
- for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
- const Double_t stat = hnor->GetBinContent(ibin+1);
- if(stat<EPSILON){
- continue;
- }
-
- const Int_t idet = ToDetector(ibin);
- const Int_t itb = ToTimeBin(ibin);
- TVectorD *vec=(TVectorD *)obj->At(idet);
- (*vec)[itb] = stat;
- }
-
- hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
- for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
- const Int_t idet = ToDetector(ibin);
- const TVectorD *vec=(TVectorD *)obj->At(idet);
-
- Int_t nonzero = 0;
- for(Int_t ii=0; ii<vec->GetNrows(); ii++){
- if((*vec)[ii]>EPSILON){
- nonzero++;
- }
- }
-
- Double_t mean = 0;
- const Double_t lowfrac = 0.6;
- //if there are too many 0's, reject this chamber by settig mean=rms=0
- if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){
- //only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large!
- mean = TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
- }
-
- hmean->SetBinContent(ibin+1, mean);
- }
-
- delete obj;
-}
-
-void AliTRDdEdxUtils::CalibOutput(const TList *lin, Int_t run)
-{
- //
- //produce calibration objects
- //
-
- TString objnames("TRDCalibHistGain TRDCalibHistT0 TRDCalibHistVd ");
- for(Int_t iter=0; iter<8; iter++){
- objnames+= GetPHQName(0, iter)+" ";
- }
-
- TList * lout = new TList;
- lout->SetOwner();
-
- TTreeSRedirector *calibStream = new TTreeSRedirector(Form("TRDCalibStream_%010d.root", run));
-
- const Int_t nh=lin->GetEntries();
- for(Int_t ii=0; ii<nh; ii++){
- const THnF *hh=(THnF*)lin->At(ii);
- const TString hname = hh->GetName();
- if(!objnames.Contains(hname))
- continue;
-
- TObjArray * cobj0 = GetCalibObj(hh, run, lout, calibStream);
- lout->Add(cobj0);
- }
-
- //lout->ls();
-
- //=============================================================
- //=============================================================
-
- TFile *fout=new TFile(Form("TRDCalibObj_%010d.root", run),"recreate");
- fout->cd();
- const Int_t nout=lout->GetEntries();
- for(Int_t ii=0; ii<nout; ii++){
- const TString oname = lout->At(ii)->GetName();
- if(oname.Contains("Obj")){
- TObjArray * cobj = (TObjArray*) lout->At(ii);
- cobj->Write(oname, TObjArray::kSingleKey);
- }
- }
- fout->Save();
- fout->Close();
- delete fout;
-
- fout=new TFile(Form("TRDCalibList_%010d.root", run),"recreate");
- fout->cd();
- lin->Write();
- lout->Write();
- fout->Save();
- fout->Close();
- delete fout;
-
- delete calibStream;
-
- /*
- http://root.cern.ch/root/html/TH1.html
- When an histogram is created, a reference to it is automatically added to the list of in-memory objects for the current file or directory. This default behaviour can be changed by:
-
- h->SetDirectory(0); for the current histogram h
- TH1::AddDirectory(kFALSE); sets a global switch disabling the reference
-
- When the histogram is deleted, the reference to it is removed from the list of objects in memory. When a file is closed, all histograms in memory associated with this file are automatically deleted.
- */
- delete lout;
-}
-
-TObjArray* AliTRDdEdxUtils::GetCalibObj(const THnF *hh, Int_t run, TList *lout, TTreeSRedirector *calibStream)
-{
- //
- //produce calibration objects
- //
-
- const TString hname = hh->GetName();
- const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;
-
- //----------------------------------------
- // Define nbin, tag, cobj0
- //----------------------------------------
- Int_t nbin =-999;
- if(hname.Contains("Gain") || hname.Contains("T0") || hname.Contains("Vd")){
- nbin = NTRDchamber();
- }
- else if(hname.Contains("PHQ")){
- nbin = NTRDtimebin();
- }
- else{
- printf("AliTRDdEdxUtils::GetCalibObj error wrong hname!! %s\n", hname.Data()); exit(1);
- }
-
- TString tag(hname);
- tag.ReplaceAll("Hist","Obj");
-
- TObjArray * cobj0 = new TObjArray(1);
- cobj0->SetOwner();
- cobj0->SetName(tag);
- cobj0->Add(new TVectorD(nbin));
-
- //----------------------------------------
- // Define lowFrac, highFrac
- //----------------------------------------
- Double_t lowFrac = -999, highFrac = -999;
- if(hname.Contains("Gain") || (hname.Contains("PHQ") && !kinvq) ){
- lowFrac = 0.01; highFrac = Q0Frac();
- }
- else if(hname.Contains("PHQ") && kinvq){
- lowFrac = Q1Frac(); highFrac = 0.99;
- }
- else{
- lowFrac = 0.01;
- highFrac = 0.99;
- }
-
- //----------------------------------------
- // Get analysis result
- //----------------------------------------
- TH1::AddDirectory(kFALSE);//important!
- TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
- TH2D *hpj = hh->Projection(1,0);
- FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
- if(hname.Contains("PHQ")){
- GetPHCountMeanRMS(hnor, htrdphmean);
- if(lout){
- lout->Add(htrdphmean);
- }
- }
- delete hpj;
-
- if(lout){
- lout->Add(hnor);
- lout->Add(hmpv);
- lout->Add(hwid);
- lout->Add(hres);
- }
-
- //----------------------------------------
- // Define Counter
- //----------------------------------------
- TVectorD *countDet=0x0;
- TObjArray *countSSL=0x0;
-
- if(hname.Contains("PHQ") && !kinvq){
- countDet=new TVectorD(540);
- countSSL=new TObjArray(90);//SectorStackLayer
- countSSL->SetOwner();
- for(Int_t ii=0; ii<90; ii++){
- countSSL->Add(new TVectorD(6));
- }
- }
-
- //----------------------------------------
- // Fill cobj0
- //----------------------------------------
-
- //ibin = binlowedge of bin(ibin+1) = the number fills this bin
- for(Int_t ibin=0; ibin<nbin; ibin++){
- Double_t gnor = hnor->GetBinContent(ibin+1);
- Double_t gmpv = hmpv->GetBinContent(ibin+1);
- Double_t gwid = hwid->GetBinContent(ibin+1);
- Double_t gres = hres->GetBinContent(ibin+1);
-
- //--- set additional cut by kpass
- Bool_t kpass = kTRUE;
- Double_t gtrdphmean = -999;
- if(htrdphmean){
- gtrdphmean = htrdphmean->GetBinContent(ibin+1);
- //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
- if(gtrdphmean<EPSILON){
- kpass = kFALSE;
- }
- if(gnor<TimeBinCountCut()*gtrdphmean){
- kpass = kFALSE;
- }
- }
-
- //--- set calibration constant p0
- Double_t p0= 0;
-
- //reason for gmpv=0:
- //1)gnor<=3; truncation in hist: (0, 0.6*ntot=1.8 with ntot=3]={1}, in hist entries can pile up so that ntot=2, or 3, and (ntot>nstart && ntot<=nstop) is skipped;
- //2)TruncatedMean(hist) out of range (only for Q0, not Q1).
-
- if(gmpv>EPSILON && kpass){
- if(tag.Contains("T0")){
- p0 = gmpv;
- }
- else{
- p0 = 1/gmpv;
- }
- //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
- }
-
- (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0;
-
- //--- save optional record
- if(p0>EPSILON && countDet && countSSL){
- const Int_t idet = ToDetector(ibin);
- (*countDet)[idet]=1;
-
- const Int_t isector = ToSector(ibin);
- const Int_t istack = ToStack(ibin);
- const Int_t ilayer = ToLayer(ibin);
- TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector);
- (*vecsectorstack)[ilayer]=1;
- }
-
- if(calibStream){
- (*calibStream)<<tag<<
- "run="<<run<<
- "p0="<<p0<<
-
- "nor="<<gnor<<
- "mpv="<<gmpv<<
- "wid="<<gwid<<
- "res="<<gres<<
- "gtrdphmean="<<gtrdphmean<<
-
- "ibin="<<ibin<<
- "\n";
- }
- }
-
- //----------------------------------------
- // Status Report
- //----------------------------------------
- if(countDet && countSSL){
- TVectorD count2Dstack(90);
- for(Int_t ii=0; ii<90; ii++){
- TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
- const Int_t nlayer = (Int_t)vecsectorstack->Sum();
- if(nlayer==6){
- count2Dstack[ii]=1;
- }
- }
-
- printf("\nAliTRDdEdxUtils::GetCalibObj Summary run: %d name: %s entries: %.0f ndetector: %03.0f n2dstack %02.0f\n\n", run, hname.Data(), hh->GetEntries(), countDet->Sum(), count2Dstack.Sum());
- }
-
- //----------------------------------------
- // Clean Up
- //----------------------------------------
-
- TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
- const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
- for(Int_t ihh=0; ihh<nhh; ihh++){
- if(!lout){
- delete (*hhs[ihh]);
- }
- }
-
- delete countDet;
- delete countSSL;
-
- //----------------------------------------
-
- return cobj0;
-}
-
-//===================================================================================
-// dEdx calculation
-//===================================================================================
-Double_t AliTRDdEdxUtils::ToyCook(const Bool_t kinvq, Int_t &ncluster, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
-{
- //
- //template for cookdedx
- //
- if(cobj){
- if(arrayQ && arrayX){
- ncluster = ApplyCalib(ncluster, arrayQ, arrayX, cobj);
- }
- else{
- printf("AliTRDdEdxUtils::ToyCook arrayQ arrayX null, applycalib can not be applied!\n"); exit(1);
- }
- }
-
- Double_t lowFrac =-999, highFrac = -999;
- if(kinvq){
- lowFrac = Q1Frac(); highFrac = 0.99;
- }
- else{
- lowFrac = 0.01; highFrac = Q0Frac();
- }
-
- Double_t meanQ = TruncatedMean(ncluster, arrayQ->GetMatrixArray(), lowFrac, highFrac);
- if(kinvq){
- if(meanQ>EPSILON){
- meanQ = 1/meanQ;
- }
- }
-
- return meanQ;
-}
-
-Double_t AliTRDdEdxUtils::CombineddEdx(const Bool_t kinvq, Int_t &concls, TVectorD *coarrayQ, TVectorD *coarrayX, const Int_t tpcncls, const TVectorD *tpcarrayQ, const TVectorD *tpcarrayX, const Int_t trdncls, const TVectorD *trdarrayQ, const TVectorD *trdarrayX)
-{
- //
- //combine tpc and trd dedx
- //
-
- for(Int_t iq=0; iq<tpcncls; iq++){
- (*coarrayQ)[iq]=(*tpcarrayQ)[iq];
- if(tpcarrayX && trdarrayX && coarrayX){
- (*coarrayX)[iq]=(*tpcarrayX)[iq];
- }
- }
- for(Int_t iq=0; iq<trdncls; iq++){
- (*coarrayQ)[tpcncls+iq]=(*trdarrayQ)[iq];
- if(tpcarrayX && trdarrayX && coarrayX){
- (*coarrayX)[tpcncls+iq]=159+(*trdarrayX)[iq];
- }
- }
-
- concls=trdncls+tpcncls;
-
- const Double_t coQ = ToyCook(kinvq, concls, coarrayQ, coarrayX);
-
- return coQ;
-}
-
-
-//===================================================================================
-// dEdx Getter and Setter
-//===================================================================================
-Double_t AliTRDdEdxUtils::GetAngularCorrection(const AliTRDseedV1 *seed)
-{
- //
- //return angular normalization factor
- //
-
- return TMath::Sqrt(1+seed->GetYref(1)*seed->GetYref(1)+seed->GetZref(1)*seed->GetZref(1));
-}
-
-Double_t AliTRDdEdxUtils::GetPadGain(const Int_t det, const Int_t icol, const Int_t irow)
-{
- //
- //get pad calibration
- //
- AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
- if (!calibration) {
- printf("AliTRDdEdxUtils::GetPadCalib No AliTRDcalibDB instance available\n"); exit(1);
- }
- AliTRDCalROC * calGainFactorROC = calibration->GetGainFactorROC(det);
- if(!calGainFactorROC){
- printf("AliTRDdEdxUtils::GetPadCalib no calGainFactorROC!\n"); exit(1);
- }
-
- Double_t padgain = -999;
- if( icol >= 0 &&
- icol < calGainFactorROC->GetNcols() &&
- irow >=0 &&
- irow < calGainFactorROC->GetNrows()){
- padgain = calGainFactorROC->GetValue(icol, irow);
- if(padgain<EPSILON){
- printf("AliTRDdEdxUtils::GetPadGain padgain 0! %f %f -- %d %d %d -- %d %d\n", padgain, EPSILON, det, icol, irow, calGainFactorROC->GetNcols(), calGainFactorROC->GetNrows()); exit(1);
- }
- }
- else{
- //printf("\nAliTRDdEdxUtils::GetPadGain warning!! indices out of range %d %d %d -- %d %d\n\n", det, icol, irow, calGainFactorROC->GetNcols(), calGainFactorROC->GetNrows() );
- }
-
- return padgain;
-}
-
-Double_t AliTRDdEdxUtils::GetRNDClusterQ(AliTRDcluster *cl)
-{
- //
- //get cluter q from GetRawQ, apply baseline and Kr pad-calibration
- //
-
- const Int_t det = cl->GetDetector();
- const Int_t pad3col = cl->GetPadCol();
- const Int_t padrow = cl->GetPadRow();
-
- const Double_t baseline = 10;
-
- Double_t rndqsum = 0;
- for(Int_t ii=0; ii<7; ii++){
- if(cl->GetSignals()[ii] < EPSILON){//bad pad marked by electronics
- continue;
- }
-
- const Int_t icol = pad3col+(ii-3);
- const Double_t padgain = GetPadGain(det, icol, padrow);
- if(padgain<0){//indices out of range, pad3col near boundary case
- continue;
- }
-
- const Double_t rndsignal = (cl->GetSignals()[ii] - baseline)/(IsPadGainOn()? padgain : 1);
-
- //sum it anyway even if signal below baseline, as long as the total is positive
- rndqsum += rndsignal;
- }
-
- return rndqsum;
-}
-
-Double_t AliTRDdEdxUtils::GetClusterQ(const Bool_t kinvq, const AliTRDseedV1 * seed, const Int_t itb)
-{
- //
- //get cluster charge
- //
- Double_t dq = 0;
- AliTRDcluster *cl = 0x0;
-
- //GetRNDClusterQ(cl)>0 ensures that the total sum of q is above baseline*NsignalPhysical.
- cl = seed->GetClusters(itb); if(cl && GetRNDClusterQ(cl)>0 ) dq+= GetRNDClusterQ(cl);//cl->GetRawQ();
- cl = seed->GetClusters(itb+AliTRDseedV1::kNtb); if(cl && GetRNDClusterQ(cl)>0 ) dq+= GetRNDClusterQ(cl);//cl->GetRawQ();
-
- dq /= GetAngularCorrection(seed);
-
- dq /= QScale();
-
- if(kinvq){
- if(dq>EPSILON){
- dq = 1/dq;
- }
- }
-
- return dq;
-}
-
-Int_t AliTRDdEdxUtils::GetArrayClusterQ(const Bool_t kinvq, TVectorD *arrayQ, TVectorD *arrayX, const AliTRDtrackV1 *trdtrack, Int_t timeBin0, Int_t timeBin1, Int_t tstep)
-{
- //
- //return nclustter
- //(if kinvq, return 1/q array), size of array must be larger than 31*6
- //
- if(!arrayQ || arrayQ->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
- printf("AliTRDdEdxUtils::GetArrayClusterQ error arrayQ null or size too small! %d\n", arrayQ? arrayQ->GetNrows() : -999); exit(1);
- }
- if(!arrayX || arrayX->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
- printf("AliTRDdEdxUtils::GetArrayClusterQ error arrayX null or size too small! %d\n", arrayX? arrayX->GetNrows() : -999); exit(1);
- }
-
- const Int_t mintb = 0;
- const Int_t maxtb = AliTRDseedV1::kNtb-1;
- if(timeBin0<mintb) timeBin0=mintb;
- if(timeBin1>maxtb) timeBin1=maxtb;
- if(tstep<=0) tstep=1;
-
- //============
- Int_t tbN=0;
- Double_t tbQ[200];
- Int_t tbBin[200];
-
- for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
- const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
- if(!seed)
- continue;
-
- const Int_t det = seed->GetDetector();
-
- for(Int_t itb=timeBin0; itb<=timeBin1; itb+=tstep){
- const Double_t dq = GetClusterQ(kinvq, seed, itb);
- if(dq<EPSILON)
- continue;
-
- const Int_t gtb = det * AliTRDseedV1::kNtb + itb;
-
- tbQ[tbN]=dq;
- tbBin[tbN]=gtb;
- tbN++;
- }
- }
-
- Int_t ncls = 0;
- for(Int_t iq=0; iq<tbN; iq++){
- if(tbQ[iq]<EPSILON)
- continue;
-
- (*arrayQ)[ncls] = tbQ[iq];
- (*arrayX)[ncls] = tbBin[iq];
-
- ncls++;
- }
-
- static Int_t kprint = 100;
- if(kprint<0){
- printf("\nAliTRDdEdxUtils::GetArrayClusterQ raw cluster-Q\n");
- for(Int_t iq=0; iq<ncls; iq++){
- const Int_t ichamber = ToLayer((*arrayX)[iq]);
- const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
- if(!seed){
- printf("error seed null!!\n"); exit(1);
- }
- const Double_t rawq = (*arrayQ)[iq] * 45. * GetAngularCorrection(seed);
- printf("esdid=%d; chamber=%d; timebin=%d; rawq= %.3f; myq[%d]= %e;\n", trdtrack->GetESDid(), ichamber, ToTimeBin((*arrayX)[iq]), rawq, iq, (*arrayQ)[iq]);
- }
- printf("\n");
- }
- kprint++;
-
- return ncls;
-}
-
-Int_t AliTRDdEdxUtils::UpdateArrayX(const Int_t ncls, TVectorD* arrayX)
-{
- //
- //arrayX det*Ntb+itb -> itb
- //
-
- TVectorD countChamber(6);
- for(Int_t ii = 0; ii<ncls; ii++){
- const Int_t xx = (Int_t)(*arrayX)[ii];
- const Int_t idet = ToDetector(xx);
-
- const Double_t ich = AliTRDgeometry::GetLayer(idet);
- const Double_t itb = ToTimeBin(xx);
- (*arrayX)[ii] = ich*AliTRDseedV1::kNtb+itb;
-
- countChamber[ich] = 1;
- }
-
- const Double_t nch = countChamber.Sum();
- return (Int_t) nch;
-}
-
-void AliTRDdEdxUtils::SetChamberQT(const AliTRDtrackV1 *trdtrack, const Int_t kcalib, THnF * hgain, THnF * ht0, THnF * hvd)
-{
- //
- //CookdEdx at TRD track level, use chamber info, related calibrations: chamber-gain; T0, Vd based on raw PH distribution
- //
-
- static Int_t kprint = 100;
-
- fgNchamber = 0;
- for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
- //initialize output, default values: 0, so that summation and weighting will automatically discard default quantities
- fgChamberQ[ichamber] = fgChamberTmean[ichamber] = 0;
-
- const AliTRDseedV1 *seed = trdtrack->GetTracklet(ichamber);
- if (!seed)
- continue;
-
- const Int_t idet = seed->GetDetector();
-
- //-------------------------------------------------------------------------
-
- Double_t qsum = 0, qtsum = 0, w2sum = 0;
- for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){
- const Double_t dq = GetClusterQ(0, seed, itb);
- if(dq<EPSILON)
- continue;
-
- qsum += dq;
- qtsum += dq*itb;
- w2sum += dq*itb*itb;
- }
- if(qsum<EPSILON)
- continue;
-
- //-------------------------------------------------------------------------
-
- Double_t tbm, tbr = 0;
- tbm = GetMeanRMS(qsum, qtsum, w2sum, &tbr);
-
- qsum /= 1.25e3/45.;
-
- if(hgain){
- const Double_t var[]={idet, qsum};
- hgain->Fill(var);
- }
- if(ht0){
- const Double_t var[]={idet, tbm};
- ht0->Fill(var);
- }
- if(hvd){
- const Double_t var[]={idet, tbr};
- hvd->Fill(var);
- }
-
- Double_t gain = 1, t0 = 0, vd = 1;
- if(kcalib){
- if(!fgObjGain) {printf("AliTRDdEdxUtils::SetChamberQT error Gain array null!!\n"); exit(1);}
- if(! fgObjT0) {printf("AliTRDdEdxUtils::SetChamberQT error T0 array null!!\n"); exit(1);}
- if(! fgObjVd) {printf("AliTRDdEdxUtils::SetChamberQT error Vd array null!!\n"); exit(1);}
-
- const TVectorD * gainvec = (TVectorD*) fgObjGain->At(0); gain = (*gainvec)[idet];
- const TVectorD * t0vec = (TVectorD*) fgObjT0->At(0); t0 = (* t0vec)[idet];
- const TVectorD * vdvec = (TVectorD*) fgObjVd->At(0); vd = (* vdvec)[idet];
- }
- if(kprint<0){
- printf("\nAliTRDdEdxUtils::CookdEdxV2\n");
- printf("idet = %d;\n", idet);
- printf("gain = %15f; t0 = %15f; vd = %15f;\n", gain, t0, vd);
- printf("\n");
- }
-
- qsum *= gain;
- tbm = (tbm-t0)*vd;
-
- if(qsum<EPSILON)
- continue;
-
- //-------------------------------------------------------------------------
-
- //should have non-zero value, initialized with default 0 (except for calibrated tbm, may be very close to 0)
- fgChamberQ[ichamber] = qsum;
- fgChamberTmean[ichamber] = tbm;
- fgNchamber++;
- }
-
- if(kprint<0){
- printf("\nAliTRDdEdxUtils::CookdEdxV2 summary:\n");
-
- printf("\nfgNchamber = %d\n", fgNchamber);
- for(Int_t ich=0; ich<AliTRDtrackV1::kNplane; ich++){
- printf("fgChamberTmean[%d] = %15f; fgChamberQ[%d] = %15f;\n", ich, fgChamberTmean[ich], ich, fgChamberQ[ich]);
- }
- }
-
- fgTrackTmean = -999;
- if(fgNchamber){
- fgTrackTmean = 0;
- for(Int_t ich=0; ich<AliTRDtrackV1::kNplane; ich++){
- fgTrackTmean += fgChamberTmean[ich];
- }
- fgTrackTmean /= fgNchamber;
- }
-
- if(kprint<0){
- printf("\nAliTRDdEdxUtils::CookdEdxV2\n");
- printf("GetTrackTmean() %15f\n", GetTrackTmean());
- printf("\n");
- }
- kprint++;
-
- return;
-}
-
-
-//===================================================================================
-// dEdx Parameterization
-//===================================================================================
-
-Double_t AliTRDdEdxUtils::Q0MeanTRDpp(const Double_t bg)
-{
- //
- //truncated Mean Q_{xx} in TRD
- //
-
- Double_t par[8];
- //03132012161150
- //opt: ppQ0
-par[0]= 2.397001e-01;
-par[1]= 1.334697e+00;
-par[2]= 6.967470e+00;
-par[3]= 9.055289e-02;
-par[4]= 9.388760e+00;
-par[5]= 9.452742e-04;
-par[6]= -1.866091e+00;
-par[7]= 1.403545e+00;
-
- ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see2.log:hhtype2Q0b2c2 scale 0.428543 at ltbg 0.650000
- return 0.428543*MeandEdxTR(&bg, par);
-}
-
-Double_t AliTRDdEdxUtils::Q0MeanTRDPbPb(const Double_t bg)
-{
- //
- //truncated Mean Q_{xx} in TRD
- //
-
- Double_t par[8];
-
- //03132012161259
- //opt: PbPbQ0
-par[0]= 1.844912e-01;
-par[1]= 2.509702e+00;
-par[2]= 6.744031e+00;
-par[3]= 7.355123e-02;
-par[4]= 1.166023e+01;
-par[5]= 1.736186e-04;
-par[6]= -1.716063e+00;
-par[7]= 1.611366e+00;
-
- ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale 0.460994 at ltbg 0.650000
- return 0.460994*MeandEdxTR(&bg, par);
-}
-
-Double_t AliTRDdEdxUtils::Q1MeanTRDpp(const Double_t bg)
-{
- //
- //truncated Mean 1/(1/Q)_{xx} in TRD
- //
-
- Double_t par[8];
-
- //So 4. Mär 13:30:51 CET 2012
- //opt: trdppQ1
- par[0]= 2.434646e-01;
- par[1]= 1.400211e+00;
- par[2]= 6.937471e+00;
- par[3]= 7.758118e-02;
- par[4]= 1.097372e+01;
- par[5]= 4.297518e-04;
- par[6]= -1.806266e+00;
- par[7]= 1.543811e+00;
-
- //hhtype2Q1b2c2 scale 0.418629 at ltbg 0.650000
-
- return 0.418629*MeandEdxTR(&bg, par);
-}
-
-Double_t AliTRDdEdxUtils::Q1MeanTRDPbPb(const Double_t bg)
-{
- //
- //truncated Mean 1/(1/Q)_{xx} in TRD
- //
-
- Double_t par[8];
-
- //So 4. Mär 13:30:52 CET 2012
- //opt: trdPbPbQ1
- par[0]= 2.193660e-01;
- par[1]= 2.051864e+00;
- par[2]= 6.825112e+00;
- par[3]= 6.151693e-02;
- par[4]= 1.390343e+01;
- par[5]= 6.010032e-05;
- par[6]= -1.676324e+00;
- par[7]= 1.838873e+00;
-
- //hhtype4Q1b2c2 scale 0.457988 at ltbg 0.650000
-
- return 0.457988*MeandEdxTR(&bg, par);
-}
-
-Double_t AliTRDdEdxUtils::QMeanTPC(const Double_t bg)
-{
- //
- //bethe bloch in TPC
- //
-
- Double_t par[5];
- //Mi 15. Feb 14:48:05 CET 2012
- //train_2012-02-13_1214.12001, tpcsig
- par[0]= 4.401269;
- par[1]= 9.725370;
- par[2]= 0.000178;
- par[3]= 1.904962;
- par[4]= 1.426576;
-
- return MeandEdx(&bg, par);
-}
-
-Double_t AliTRDdEdxUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
-{
- //
- //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
- //npar = 8 = 3+5
- //
- Double_t ptr[4]={0};
- for(int ii=0; ii<3; ii++){
- ptr[ii+1]=pin[ii];
- }
- return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
-}
-
-Double_t AliTRDdEdxUtils::MeanTR(const Double_t * xx, const Double_t * par)
-{
- //
- //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
- //npar = 4
- //
-
- const Double_t bg = xx[0];
- const Double_t gamma = sqrt(1+bg*bg);
-
- const Double_t p0 = TMath::Abs(par[1]);
- const Double_t p1 = TMath::Abs(par[2]);
- const Double_t p2 = TMath::Abs(par[3]);
-
- const Double_t zz = TMath::Log(gamma);
- const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
-
- return par[0]+tryield;
-}
-
-Double_t AliTRDdEdxUtils::MeandEdx(const Double_t * xx, const Double_t * par)
-{
- //
- //ALEPH parametrization for dEdx
- //npar = 5
- //
-
- const Double_t bg = xx[0];
- const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
-
- const Double_t p0 = TMath::Abs(par[0]);
- const Double_t p1 = TMath::Abs(par[1]);
- const Double_t p2 = TMath::Abs(par[2]);
- const Double_t p3 = TMath::Abs(par[3]);
- const Double_t p4 = TMath::Abs(par[4]);
-
- const Double_t aa = TMath::Power(beta, p3);
- const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
-
- //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
-
- return (p1-aa-bb)*p0/aa;
-}
-
-Double_t AliTRDdEdxUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
-{
- //
- //f(x)-> f(y) with y=log10(x)
- //
- const Double_t x2[]={TMath::Power(10, xx[0])};
- return func(x2, par);
-}
-
-//===================================================================================
-// Detector, Data and Control Constant
-//===================================================================================
-Int_t AliTRDdEdxUtils::ToDetector(const Int_t gtb)
-{
- //
- //gtb = det*Ntb+itb
- //
- return gtb/AliTRDseedV1::kNtb;
-}
-
-Int_t AliTRDdEdxUtils::ToTimeBin(const Int_t gtb)
-{
- //
- //gtb = det*Ntb+itb
- //
- return gtb%AliTRDseedV1::kNtb;
-}
-
-Int_t AliTRDdEdxUtils::ToSector(const Int_t gtb)
-{
- //
- //return sector
- //
- return AliTRDgeometry::GetSector(ToDetector(gtb));
-}
-
-Int_t AliTRDdEdxUtils::ToStack(const Int_t gtb)
-{
- //
- //return stack
- //
- return AliTRDgeometry::GetStack(ToDetector(gtb));
-}
-
-Int_t AliTRDdEdxUtils::ToLayer(const Int_t gtb)
-{
- //
- //return layer
- //
- return AliTRDgeometry::GetLayer(ToDetector(gtb));
-}
-
-TString AliTRDdEdxUtils::GetRunType(const Int_t run)
-{
- //
- //return run type
- //
-
- TString type;
- if(run>=121527 && run<= 126460)//LHC10d
- type="pp2010LHC10d";
- else if(run>=126461 && run<= 130930)//LHC10e
- type="pp2010LHC10e";
- else if(run>=136782 && run <= 139846)//LHC10h
- type="PbPb2010LHC10h";
- else if(run>= 143224 && run<= 143237)//2011Feb
- type="cosmic2011Feb";
- else if(run>= 150587 && run<= 154930){
- type="cosmic2011MayJun";
-
- TString runstr(Form("%d",run));
- const TString listrun1kg("154601 154602 154629 154634 154636 154639 154643");
- if(listrun1kg.Contains(runstr)){
- type+="1kG";
- }
- else{
- type+="5kG";
- }
- }
- else{
- type="unknown";
- }
-
- type.ToUpper();
- return type;
-}
-
-void AliTRDdEdxUtils::PrintControl()
-{
- //
- //print out control variable
- //
- printf("\nAliTRDdEdxUtils::PrintControl Q0Frac %.1f, Q1Frac %.1f, TimeBinCountCut %.2f, CalibTPCnclsCut %d, IsExBOn %d, IsPadGainOn %d, QScale %.2f\n\n", Q0Frac(), Q1Frac(), TimeBinCountCut(), CalibTPCnclsCut(), IsExBOn(), IsPadGainOn(), QScale());
-}
-//===================================================================================
-//===================================================================================
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-//
-//
-//
-// Xianguo Lu
-// lu@physi.uni-heidelberg.de
-// Xianguo.Lu@cern.ch
-//
-/*
-grep " AliTRDdEdxUtils::" AliTRDdEdxUtils.cxx | grep "=" -v | grep -v "[6]" | grep -v printf |wc
-grep "(" AliTRDdEdxUtils.h | grep ";" | grep -v grep | grep -v ClassDef | grep -v "{" | grep -v typedef | wc
-*/
-
-
-#ifndef ALITRDDEDXUTILS_H
-#define ALITRDDEDXUTILS_H
-
-#ifndef TVECTORD_H
-#include "TVectorD.h"
-#endif
-
-#ifndef THNSPARSE_H
-#include "THn.h"
-#include "THnSparse.h"
-#endif
-
-#ifndef TTREESTREAM_H
-#include "TTreeStream.h"
-#endif
-
-class TH1D;
-class TH2D;
-class TObjArray;
-
-class AliESDEvent;
-class AliESDtrack;
-class AliTRDcluster;
-class AliTRDtrackV1;
-class AliTRDseedV1;
-
-class AliTRDdEdxUtils
-{
- public:
-
- //===================================================================================
- // Math and Histogram
- //===================================================================================
- static void GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres);
- static Double_t GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms=0x0, Double_t * gerr=0x0);
- static Double_t TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms=0x0, Double_t * gerr=0x0, Double_t *wws=0x0);
- static Double_t TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms=0x0, Double_t * gerr=0x0);
- static void FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, const Double_t thres, const Double_t lowfrac, const Double_t highfrac);
-
- //===================================================================================
- // TRD Analysis Fast Tool
- //===================================================================================
- static Int_t GetNtracklet(const AliESDEvent *esd);
- static AliTRDtrackV1 * GetTRDtrackV1(const AliESDtrack * esdtrack);
- static Bool_t IsInSameStack(const AliTRDtrackV1 *trdtrack);
- static Bool_t GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom);
-
- //===================================================================================
- // Calibration
- //===================================================================================
- static void SetCalibFile(const TString file){fgCalibFile = file;}
- static void DeleteCalibObj();
- static void IniCalibObj();
-
- static void SetObjPHQ(TObjArray * obj){fgObjPHQ = obj;}
- static Bool_t GenerateDefaultPHQOCDB(const TString path="local://./");
-
- static Double_t GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig);
-
- static void DeleteCalibHist();
- static void IniCalibHist(TList *list, const Bool_t kPHQonly=kFALSE);
- static Bool_t ReadCalibHist(const TString filename, const TString listname);
-
- static TObjArray * GetObjPHQ();
- static TObjArray * GetHistPHQ(){return fgHistPHQ;}
- static TObjArray * GetObjPHQ(const Bool_t kinvq, const Double_t mag, const Int_t charge);
- static THnF * GetHistPHQ(const Bool_t kinvq, const Double_t mag, const Int_t charge);
-
- static void FillCalibHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnF * hcalib, const Double_t scale);
- static void FillCalibHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale);
-
- static void CalibOutput(const TList *l0, Int_t run);
- static TObjArray* GetCalibObj(const THnF *hh, Int_t run=-999, TList *lout=0x0, TTreeSRedirector *calibStream=0x0);
-
- static THnF * GetHistGain(){return fgHistGain;}
- static THnF * GetHistT0(){return fgHistT0;}
- static THnF * GetHistVd(){return fgHistVd;}
-
- //===================================================================================
- // dEdx calculation
- //===================================================================================
- static Double_t ToyCook(const Bool_t kinvq, Int_t &ncluster, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj=0x0);
- static Double_t CombineddEdx(const Bool_t kinvq, Int_t &concls, TVectorD *coarrayQ, TVectorD *coarrayX, const Int_t tpcncls, const TVectorD *tpcarrayQ, const TVectorD *tpcarrayX, const Int_t trdncls, const TVectorD *trdarrayQ, const TVectorD *trdarrayX);
-
- //===================================================================================
- // dEdx Getter and Setter
- //===================================================================================
- static Int_t GetArrayClusterQ(const Bool_t kinvq, TVectorD *arrayQ, TVectorD *arrayX, const AliTRDtrackV1 *trdtrack, Int_t timeBin0=-1, Int_t timeBin1=1000, Int_t tstep=1);
- static Int_t UpdateArrayX(const Int_t ncls, TVectorD* arrayX);
- static void SetChamberQT(const AliTRDtrackV1 *trdtrack, const Int_t kcalib, THnF * hgain=0x0, THnF * ht0=0x0, THnF * hvd=0x0);
-
- static Int_t GetNchamber() {return fgNchamber;}
- static Double_t GetChamberQ(const Int_t ich) {return fgChamberQ[ich];}
- static Double_t GetChamberTmean(const Int_t ich) {return fgChamberTmean[ich];}
- static Double_t GetTrackTmean() {return fgTrackTmean;}
-
- //===================================================================================
- // dEdx Parameterization
- //===================================================================================
-
- static Double_t MeandEdx(const Double_t * xx, const Double_t * par);
- static Double_t MeanTR(const Double_t * xx, const Double_t * par);
- static Double_t MeandEdxTR(const Double_t * xx, const Double_t * par);
-
- static Double_t QMeanTPC(const Double_t bg);
- static Double_t Q0MeanTRDpp(const Double_t bg);
- static Double_t Q1MeanTRDpp(const Double_t bg);
- static Double_t Q0MeanTRDPbPb(const Double_t bg);
- static Double_t Q1MeanTRDPbPb(const Double_t bg);
-
- typedef Double_t (*FFunc)(const Double_t *xx, const Double_t *par);
-
- static Double_t MeandEdxLogx(const Double_t * xx, const Double_t * par){return ToLogx(MeandEdx, xx, par);}
- static Double_t MeanTRLogx(const Double_t * xx, const Double_t * par){return ToLogx(MeanTR, xx, par);}
- static Double_t MeandEdxTRLogx(const Double_t * xx, const Double_t * par){return ToLogx(MeandEdxTR, xx, par);}
-
- //===================================================================================
- // Detector, Data and Control Constant
- //===================================================================================
- static Int_t NTRDchamber(){return 18*5*6;} //540
- static Int_t NTRDtimebin(){return NTRDchamber()*31;} //16740
- static Int_t ToDetector(const Int_t gtb);
- static Int_t ToTimeBin(const Int_t gtb);
- static Int_t ToSector(const Int_t gtb);
- static Int_t ToStack(const Int_t gtb);
- static Int_t ToLayer(const Int_t gtb);
-
- static TString GetRunType(const Int_t run);
-
- static void SetPadGainOn(const Bool_t kon){ fgPadGainOn = kon; }
- static void SetExBOn(const Bool_t kon){ fgExBOn = kon; }
- static void SetQScale(const Double_t scale){ fgQScale = scale; }
- static void SetQ0Frac(const Double_t q0){ fgQ0Frac = q0; }
- static void SetQ1Frac(const Double_t q1){ fgQ1Frac = q1; }
- static void SetTimeBinCountCut(const Double_t tbc){ fgTimeBinCountCut = tbc; }
- static void SetCalibTPCnclsCut(const Int_t tpc){ fgCalibTPCnclsCut = tpc; }
-
- static Bool_t IsPadGainOn(){return fgPadGainOn;}
- static Bool_t IsExBOn(){return fgExBOn;}
- static Double_t QScale(){return fgQScale;}
- static Double_t Q0Frac(){return fgQ0Frac;}
- static Double_t Q1Frac(){return fgQ1Frac;}
- static Double_t TimeBinCountCut(){return fgTimeBinCountCut;}
- static Int_t CalibTPCnclsCut(){return fgCalibTPCnclsCut;}
-
- static void PrintControl();
- //===================================================================================
- //===================================================================================
-
- private:
-
- //dEdx Getter and Setter
- static Double_t GetAngularCorrection(const AliTRDseedV1 *seed);
- static Double_t GetPadGain(const Int_t det, const Int_t icol, const Int_t irow);
- static Double_t GetRNDClusterQ(AliTRDcluster *cl);
- static Double_t GetClusterQ(const Bool_t kinvq, const AliTRDseedV1 * seed, const Int_t itb);
-
- //dEdx Parameterization
- static Double_t ToLogx(FFunc func, const Double_t * xx, const Double_t * par);
-
- //Calibration
- static Int_t ApplyCalib(const Int_t nc0, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj);
- static void GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean);
- static Int_t GetPHQIterator(const Bool_t kinvq, const Double_t mag, const Int_t charge);
- static TString GetPHQName(const Bool_t kobj, const Int_t iter);
-
- //Detector, Data and Control Constant
- static THnF *fgHistGain;//PH hist
- static THnF *fgHistT0;//PH hist
- static THnF *fgHistVd;//PH hist
- static TObjArray * fgHistPHQ;//array containing 8 THnF!
-
- static TString fgCalibFile; //private path for calibration object
-
- static TObjArray * fgObjGain;//chamber gain obj
- static TObjArray * fgObjT0;//t0 obj
- static TObjArray * fgObjVd;//Vd obj
- static TObjArray * fgObjPHQ;//array containing 8 TObjArray!
-
- static Int_t fgNchamber; //number of chamber used in cookdedx
- static Double_t fgChamberQ[6]; //dqdl in chamber [i]
- static Double_t fgChamberTmean[6]; //Q-weighted timebin \sum Q*T / \sum Q
-
- static Double_t fgTrackTmean; //mean timebin over track
-
- static Bool_t fgPadGainOn; //pad gain
- static Bool_t fgExBOn; //exbon
- static Double_t fgQScale; //Qscale
- static Double_t fgQ0Frac; //q0frac
- static Double_t fgQ1Frac; //q1frac
- static Double_t fgTimeBinCountCut; //tbcut
- static Int_t fgCalibTPCnclsCut; //tpccut
-};
-
-#endif
#include "AliTRDReconstructor.h"
#include "AliTRDPIDResponse.h"
#include "AliTRDrecoParam.h"
-#include "AliTRDdEdxUtils.h"
+#include "AliTRDdEdxBaseUtils.h"
+#include "AliTRDdEdxCalibHistArray.h"
+#include "AliTRDdEdxCalibUtils.h"
+#include "AliTRDdEdxReconUtils.h"
ClassImp(AliTRDtrackV1)
static Int_t nprint = 0;
if(!nprint){
- AliTRDdEdxUtils::PrintControl();
+ AliTRDdEdxBaseUtils::PrintControl();
nprint++;
}
//do truncated mean
- AliTRDdEdxUtils::SetObjPHQ(AliTRDcalibDB::Instance()->GetPHQ());
- const Double_t mag = AliTRDdEdxUtils::IsExBOn() ? GetBz() : -1;
- const Double_t charge = AliTRDdEdxUtils::IsExBOn() ? Charge() : -1;
+ AliTRDdEdxCalibUtils::SetObjArray(AliTRDcalibDB::Instance()->GetPHQ());
+ const Double_t mag = AliTRDdEdxBaseUtils::IsExBOn() ? GetBz() : -1;
+ const Double_t charge = AliTRDdEdxBaseUtils::IsExBOn() ? Charge() : -1;
fTruncatedMean = CookTruncatedMean(0, mag, charge, kTRUE, fNchamberdEdx, fNclusterdEdx);
return kTRUE;
//
TVectorD arrayQ(200), arrayX(200);
- ncls = AliTRDdEdxUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, this, timeBin0, timeBin1, tstep);
+ ncls = AliTRDdEdxReconUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, this, timeBin0, timeBin1, tstep);
- const TObjArray *cobj = kcalib ? AliTRDdEdxUtils::GetObjPHQ(kinvq, mag, charge) : NULL;
+ const TObjArray *cobj = kcalib ? AliTRDdEdxCalibUtils::GetObj(kinvq, mag, charge) : NULL;
- const Double_t tmean = AliTRDdEdxUtils::ToyCook(kinvq, ncls, &arrayQ, &arrayX, cobj);
+ const Double_t tmean = AliTRDdEdxReconUtils::ToyCook(kinvq, ncls, &arrayQ, &arrayX, cobj);
- nch = AliTRDdEdxUtils::UpdateArrayX(ncls, &arrayX);
+ nch = AliTRDdEdxReconUtils::UpdateArrayX(ncls, &arrayX);
if(Qs && Xs){
(*Qs)=arrayQ;
AliTRDtransform.cxx
AliTRDtrackletOflHelper.cxx
AliTRDpidUtil.cxx
- AliTRDdEdxUtils.cxx
+ AliTRDdEdxCalibHistArray.cxx
+ AliTRDdEdxBaseUtils.cxx
+ AliTRDdEdxCalibUtils.cxx
+ AliTRDdEdxReconUtils.cxx
AliTRDpidESD.cxx
AliTRDReconstructor.cxx
AliTRDseedV1.cxx AliTRDtrackV1.cxx
#pragma link C++ class AliTRDpidESD+;
#pragma link C++ class AliTRDpidUtil+;
-#pragma link C++ class AliTRDdEdxUtils+;
+#pragma link C++ class AliTRDdEdxCalibHistArray+;
+#pragma link C++ class AliTRDdEdxBaseUtils+;
+#pragma link C++ class AliTRDdEdxCalibUtils+;
+#pragma link C++ class AliTRDdEdxReconUtils+;
#pragma link C++ class AliTRDReconstructor+;
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff