1 #if !defined(__CINT__) || defined(__MAKECINT__)
12 #include <TGraphErrors.h>
16 #include "AliCDBManager.h"
17 #include "AliCDBEntry.h"
18 #include "AliCDBStorage.h"
19 #include "AliITSsegmentationSDD.h"
20 #include "AliITSCorrMapSDD.h"
21 #include "AliITSCorrMap1DSDD.h"
22 #include "AliITSDriftSpeedArraySDD.h"
23 #include "AliITSDriftSpeedSDD.h"
24 #include "AliITSresponseSDD.h"
27 Bool_t verbose = kFALSE;
31 DeltaX_L = -deltaX + deltaV_L*(T_L - T0_true - deltaT0) - V_Ltrue*deltaT0 - corrMapL
32 DeltaX_R = deltaX + deltaV_R*(T_R - T0_true - deltaT0) - V_Rtrue*deltaT0 - corrMapR
36 corrMapL(R) : Vdrift vs X nonuniformity correction maps
37 V_L(R)true: true drift speed
38 deltaV_L(R) : error in drift speed
39 V_Ltrue = V_Lassumed - deltaV_L
40 V_Rtrue = V_Rassumed - deltaV_R
42 DeltaX : measured residuals (track extrapolation - measurement) for left (L)
43 and right (R) sides (right side multiplied by -1!!!)
44 deltaX : misalignment of module (Xtrue = Xassumed - deltaX)
45 T_L(R) : measured drift time (w/o T0 subtraction)
50 //------------------- Details of CDB objects to create ---------------
53 TString cdbComment = "";
54 //--------------------------------------------------------------------
56 const Int_t kSDDMin=240,kSDDMax=499,kNSDD=kSDDMax-kSDDMin+1;
57 enum {kSXvsX,kSXvsZ,kSTDvsX,kSTDvsZ,kSDXvsXclean,kSVDvsX,kSCorrMapX,kSXvsXCorr,kSXvsZCorr,kSDVvsZ,kSDVvsZOrig,kSGrDxx,kSGrTDx,kSGrTXCorr, kNStore};
58 const Int_t maxVZOrd = 4; // max order of polinomial for VD vs Z correction (normal sensors)
59 const Double_t chi2BadVZ = 2.9; // treat as bad module if chi2/ndf is worse
61 enum {kDxX,kDxZ,kTDX,kTDZ,kNQATypes};
62 // Prefixes for residuals histos prepared by PrepSDDResidMap.C from MP2 output
63 // For QA histograms the name should be changed to "hpSDDResXvsXD","hpSDDResXvsZ","hpSDDDrTimevsXD","hpSDDDrTimevsZ" respectively
64 const Char_t* kQAHistoNames[] = { //
65 "hpSDDResXvsXD", // DX vs Xdrift
66 "hpSDDResXvsZ", // DX vs local Z
67 "hpSDDDrTimevsXD", // TDrift-T0 vs Xdrift
68 "hpSDDDrTimevsZ" // TDrift-T0 vs local Z
70 const Char_t* kLRNames[] = {"0","1"}; // identifiers for left/right sides in histo names
72 const Char_t* kQAHistoNames[] = { //
73 "Xdrift_Mod_", // DX vs Xdrift
74 "ZAnode_Mod_", // DX vs local Z
75 "TD_vs_Xdrift_Mod_", // TDrift-T0 vs Xdrift
76 "TD_vs_ZAnode_Mod_" // TDrift-T0 vs local Z
78 const Char_t* kLRNames[] = {"left","right"}; // identifiers for left/right sides in histo names
81 const double kMaxChi2SimpleMap = 1.0; // if polinomial fit gives chi2 below this value, make map out of it
82 double minBEntryX = 10; // min entries in the bin to account (vs Xdrift)
83 double minBEntryZ = 20; // min entries in the bin to account (vs Z)
84 double skipDXEdge = 1000.; // microns to skip from each side of XDrift when smoothing
85 double edgeSmearMinErr = 60.; // minimal error for extrapolation
86 double edgeSmearMaxErr = 600.; // maximal error for extrapolation
87 double wDXEdge = 5000.; // width in microns for XDrift edges to average
88 double wVDEdge = 2500; // width in microns for VDrift averaging region on the edges
89 double threshClean= 0.3; // clean edge bins with occupancy smaller than in the neighbourhood
90 int nVDEdge = 20; // min width in n bins for VDrift averaging region on the edges
91 int smoothLevel= 5; // smoothing level
92 int minBinsEdge = 5; // min number of X>0 good bins to fit the edge
93 Bool_t useSplieForR2V = kFALSE; // if true, extract Vd profile using Derivateve of spline (may be inadequate for sharp turns)
96 Bool_t userLeftRightFromTASK = kTRUE; // VERY IMPORTANT: histos produced by AliAnalysisTaskITSAlignQA have inverted left/right side definitions
97 Bool_t forceT0CorrTo0 = kFALSE;//kTRUE;
98 Bool_t forceRPhiCorrTo0 = kTRUE;
99 Bool_t userDummyCorrMap = kFALSE;//kTRUE;
100 Bool_t userDummyt0Corr = kFALSE;//kTRUE;
101 Bool_t userDummyxlCorr = kFALSE;//kTRUE;
102 Int_t userRebinX = 1;
103 Int_t userRebinZ = 2;
106 AliITSsegmentationSDD sddSeg;
108 //----------------- working variables --------------------------------
109 Int_t currSDDId=-1,currSDDSide=-1,currMod=-1;
110 TProfile* currHDXvsX = 0; // raw DX vs X histo from QA
111 TProfile* currHDXvsZ = 0; // raw DX vs Z histo from QA
112 TProfile* currHTDvsX = 0; // raw DriftTime-T0 vs X histo from QA
113 TProfile* currHTDvsZ = 0; // raw DriftTime-T0 vs Z histo from QA
114 TH1* currHDXvsXclean = 0; // smoothed DX vs X histo
115 TH1* currHVDvsX = 0; // extracted VDrift vs X profile
116 TH1* currHCorrMapX = 0; // extracted correction map
117 TH1* currHDXvsXCorr = 0; // DX vs X histo after application of the map
118 TH1* currHDXvsZCorrMT0 = 0; // DX vs X histo after application of the map and t0 correction
119 TH1* currHDVvsZ = 0; // correction for the VDrift vs Z
120 TH1* currHDVvsZOrig = 0; // correction for the VDrift vs Z (before error processing)
121 TGraphErrors* currGrDxx = 0; // DX final vs XDrift
122 TGraphErrors* currGrTDx = 0; // DX final vs TDrift
123 TGraphErrors* currGrTXCorr = 0; // TDrift vs XDritf
124 TCollection* qaHistos = 0; // input QA collection
125 TObjArray procHistos; // processed histos buffer
126 AliITSresponseSDD* sddResp=0; // SDD response, updated
127 TObjArray *vdarrayOld = 0; // olf VDrifts array
128 TObjArray* vDriftArr=0; // VDrifts array, updated
129 TObjArray* corrMaps=0; // holder for correction maps
131 TH1* resOffsDXraw[2]={0}; // Dx vs X offset at Xdrift=0 raw
132 TH1* resOffsDXAMap[2]={0}; // Dx vs X offset (mean) after correction by map
133 TH1* resOffsDX[2]={0}; // Dx vs X offset at Xdrift=0 after correction by map (for each side)
134 TH1* resVDCorr[2]={0}; // VDrift correction (for each side)
135 TH1* resVDMean[2]={0}; // average VDrift (for each side)
136 TH1* resVDCorrZ[2]={0}; // mean VDrift vs Z corrections (for each side)
137 TH1* resT0Corr = 0; // correction to modules T0
138 TH1* resXLocCorr = 0; // correction to module location
139 TCanvas* sddCanv = 0; // report canvas
141 Bool_t sddRespUpdT0OK = kFALSE; // flag that SDDresponse object T0 was updated
142 Bool_t sddRespUpdVDOK = kFALSE; // flag that SDDresponse object VDrift correction was updated
143 Bool_t sddVDriftUpdOK = kFALSE; // flag that SDD VDrift object was updated
145 TString pathSDDRespOld = "";
146 TString pathSDDVDriftOld = "";
147 TString pathSDDCorrMapOld = "";
149 //--------------------------------------------------------------------
150 void Process(const char* pathSDDResp=0,const char* pathSDDVDrift=0, const char* pathSDDCorrMap=0);
151 void Process(TCollection* qa, const char* pathSDDResp=0,const char* pathSDDVDrift=0, const char* pathSDDCorrMap=0);
152 Bool_t ProcSDDSideDXvsX();
153 Bool_t ProcSDDSideDXvsZ();
154 TProfile* GetQAHisto(Int_t qaType);
155 TProfile* H2Profile(TH2* h2);
156 TH1* H2ProfileAsTH1(TH2* h2);
157 TH1* GetProfHEntries(TProfile* prof);
158 TH1* ProfileAsTH1(TProfile* prof, const char* addName);
159 TH1* CleanProfile(TProfile* profR);
160 TH1* Vdrift2Resid(TH1* vd);
161 TH1* Resid2Vdrift(TH1* res);
162 void RedoProfileErrors(TH1* profH1,TProfile* prof);
163 void SetModuleID(Int_t mdID,Int_t side=-1);
164 void PrepareModuleHistos(Int_t mdID, Int_t side, Bool_t vsZ);
165 void CheckResultDX();
166 void CalcDXCorrections();
167 double GetVOld(double z);
168 Int_t GetStoreID(int type, int imd=-1,int side=-1);
171 double ftVdZ(double *x, double *par);
172 void PlotReport(const char* psname="repSDDQA.ps");
173 TH1* GetPadBaseHisto(TPad* pad);
174 Bool_t PlotHisto(TH1* h, Option_t* opt="", int mrkStyle=20,int mrkCol=kBlack, double mrkSize=1.);
175 TLatex* AddPadLabel(const char*txt,float x=0.1,float y=0.9,int color=kBlack,float size=0.04);
176 void GetHistoMinMaxInRange(TH1* h, double &mn,double &mx);
177 double edgeLow(double* x, double *par);
178 double edgeHigh(double* x, double *par);
179 void CureEdges(TH1* prof);
180 void SafeRebin(TProfile* prof, Int_t factor, Bool_t xprof);
181 TH1* FitDXEdges(TProfile* prof);
182 Bool_t TestMapFunction(TH1* smap, TF1* fun, double lft, double rgt);
183 double ftPolComb(double* x, double *par);
184 TH1* SimpleMap(TH1* prof);
186 Bool_t LoadSDDVDrift(TString& path, TObjArray *&arr);
187 Bool_t LoadSDDResponse(TString& path, AliITSresponseSDD *&resp);
188 Bool_t LoadSDDCorrMap(TString& path, TObjArray *&maps);
189 AliCDBEntry* GetCDBEntry(const char* path);
191 void UpdateSDDResponse(AliITSresponseSDD *resp, Bool_t t0=kTRUE, Bool_t vdrift=kTRUE);
192 void UpdateSDDVDrift(AliITSDriftSpeedArraySDD* vdArr, int imd, int side);
193 TObjArray* UpdateSDDVDrift();
194 TObjArray* CreateCorrMaps();
195 AliITSresponseSDD* UpdateSDDResponse(Bool_t t0=kTRUE, Bool_t vdrift=kTRUE);
196 AliITSCorrMap1DSDD* CreateCorrMap(TH1* mapHisto, int imd, int side, AliITSCorrMap1DSDD* updateMap=0);
197 void PrepCDBObj(TObject *obj,const char* path,int firstrun=0,int lastrun=999999999,const char* comment="");
199 //-------------------------------------------------------------------
201 //_________________________________________________________________________
202 void Process(TCollection* qa, const char* pathSDDResp, const char* pathSDDVDrift, const char* pathSDDCorrMap)
206 Process(pathSDDResp,pathSDDVDrift,pathSDDCorrMap);
209 //_________________________________________________________________________
210 void Process(const char* pathSDDResp, const char* pathSDDVDrift, const char* pathSDDCorrMap)
214 pathSDDRespOld = pathSDDResp;
215 pathSDDVDriftOld = pathSDDVDrift;
216 pathSDDCorrMapOld = pathSDDCorrMap;
217 sddRespUpdT0OK = sddRespUpdVDOK = sddVDriftUpdOK = kFALSE;
219 if (pathSDDVDriftOld.IsNull()) printf("Attention: Old VDrift is missing!\n");
221 for (int imd=kSDDMin;imd<=kSDDMax;imd++) {
222 for (int ix=0;ix<2;ix++) {
223 CleanPrev(); // clean data from previous module
224 printf("Processing %d/%d\n",imd,ix);
225 PrepareModuleHistos(imd,ix, kTRUE);
226 ProcSDDSideDXvsX(); // process DX vs X histos to get corr.maps, deltaT0, deltaX0, deltaV_mean
229 CalcDXCorrections(); // calculate delta's
231 for (int ix=0;ix<2;ix++) {
232 CleanPrev(); // clean data from previous module
233 PrepareModuleHistos(imd,ix, kFALSE);
234 ProcSDDSideDXvsZ(); // process deltaV vs anode profile
239 corrMaps = CreateCorrMaps(); // create correction maps
240 PrepCDBObj(corrMaps,"ITS/Calib/MapsTimeSDD",firstRun,lastRun,cdbComment.Data());
242 if (!pathSDDVDriftOld.IsNull()) {
243 vDriftArr = UpdateSDDVDrift();
244 PrepCDBObj(vDriftArr,"ITS/Calib/DriftSpeedSDD",firstRun,lastRun,cdbComment.Data());
247 if (!pathSDDRespOld.IsNull()) {
248 sddResp = UpdateSDDResponse();
249 PrepCDBObj(sddResp,"ITS/Calib/RespSDD",firstRun,lastRun,cdbComment.Data());
255 //_________________________________________________________________________
256 Bool_t ProcSDDSideDXvsX()
258 // process one side of the SDD module
259 currHDXvsXCorr = ProfileAsTH1(currHDXvsX,"corrCheck");
260 RedoProfileErrors(currHDXvsXCorr,currHDXvsX);
262 if ( (currHDXvsXclean=CleanProfile(currHDXvsX)) ) {
264 delete currHDXvsXCorr;
265 currHDXvsXCorr = (TH1*)currHDXvsXclean->Clone( Form("%s_corrCheck",currHDXvsX->GetName()) );
267 // try simple solution
268 if (!(currHCorrMapX=SimpleMap(currHDXvsXclean))) {
269 currHVDvsX = Resid2Vdrift(currHDXvsXclean);
270 currHCorrMapX = Vdrift2Resid(currHVDvsX);
273 currHDXvsXCorr->Add(currHCorrMapX,-1);
281 //_________________________________________________________________________________
282 Bool_t ProcSDDSideDXvsZ()
284 // extract correction for Vdrift vs Z from DXvsZ profile and mean_drift_time vs Z
286 static TF1* fitP0 = new TF1("fitP0","pol0",-4000,4000);
287 double chi2s[maxVZOrd+1] = {0};
288 static TF1* fitVvsZs[maxVZOrd+1] = {0};
289 static Bool_t iniDone = kFALSE;
290 double zrange = sddSeg.Dz()/2 - 1.;
293 for (int iord=0;iord<=maxVZOrd;iord++) {
294 fitVvsZs[iord] = new TF1("fitVvsZ",ftVdZ, -zrange, zrange ,iord+2);
295 fitVvsZs[iord]->FixParameter(0, iord+0.1);
299 int nb = currHDXvsZ->GetNbinsX();
300 if (currHDXvsZ->GetEntries()<minBEntryZ*nb) return kFALSE;
302 currHDVvsZ = ProfileAsTH1(currHDXvsZ,"_vzcorr");
305 currHDXvsZCorrMT0 = ProfileAsTH1(currHDXvsZ,"_zcorrMapT0");
306 currHDXvsZCorrMT0->Reset();
308 int nbUse=0, ib0 = currHDVvsZ->FindBin(-zrange), ib1 = currHDVvsZ->FindBin( zrange);
309 double *statZB = new double[nb+1]; // entries per Z bin
310 memset(statZB,0,sizeof(double)*(nb+1));
311 double vmean=0, vmeanE=0, norm=0;
313 for (int ib=ib0;ib<=ib1;ib++) {
314 double ne = currHTDvsZ->GetBinEntries(ib);
316 double dx = currHDXvsZ->GetBinContent(ib); // raw residual X vs Z
317 double dxe= currHDXvsZ->GetBinError(ib);
318 double t = currHTDvsZ->GetBinContent(ib); // mean assumed (TDrift - T0)
319 double te = currHTDvsZ->GetBinError(ib);
320 double vCorr = resVDCorr[currSDDSide]->GetBinContent(currMod+1);
321 dx -= vCorr*t; // subtract mean V correction
322 dxe = TMath::Sqrt(dxe*dxe + vCorr*vCorr*te*te);
323 currHDXvsZCorrMT0->SetBinContent(ib,dx);
324 currHDXvsZCorrMT0->SetBinError(ib,dxe);
326 currHDXvsZCorrMT0->Fit(fitP0,"q0","");
327 double pedestal = fitP0->GetParameter(0);
329 for (int ib=ib0;ib<=ib1;ib++) {
330 double ne = currHTDvsZ->GetBinEntries(ib);
331 if (ne<minBEntryZ) continue;
332 double dx = currHDXvsZCorrMT0->GetBinContent(ib); // raw residual X vs Z
333 double dxe= currHDXvsZCorrMT0->GetBinError(ib);
334 double t = currHTDvsZ->GetBinContent(ib); // mean assumed (TDrift - T0)
335 double te = currHTDvsZ->GetBinError(ib);
338 // account for the corrections leading the mean DX to 0
341 double ve = TMath::Sqrt(dxe*dxe + v*v*te*te)/t;
344 vmeanE += 1./(ve*ve);
346 // printf("%d %f %f | %f %f | -> %f %f\n",ib,dx,dxe,t,te,v,ve);
347 int ibDest = currSDDSide==0 ? ib : nb+1-ib;
349 currHDVvsZ->SetBinContent(ibDest,v);
350 currHDVvsZ->SetBinError(ibDest,ve);
354 if (nbUse<maxVZOrd+1) {delete statZB; return kFALSE;}
355 if (vmeanE>0) { vmean /= vmeanE; vmeanE = 1./TMath::Sqrt(vmeanE); }
356 // fill empty and bad bins by mean value with large error
359 for (int ib=ib0;ib<=ib1;ib++) {
360 //if (statZB[ib]<minBEntryZ) currHDVvsZ->SetBinContent(ib,vmean);
361 if (currHDVvsZ->GetBinError(ib)>vmeanE || statZB[ib]<minBEntryZ) {
362 currHDVvsZ->SetBinError(ib,vmeanE*2);
363 currHDVvsZ->SetBinContent(ib,vmean);
368 currHDVvsZOrig = (TH1*) currHDVvsZ->Clone(Form("%s_orig",currHDVvsZ->GetName()));
369 // find leftmost good bin
371 double errL=0,errR=0, valL=0,valR=0;
372 for (bL=ib0;bL<=ib1;bL++) {
373 if (/*currHDVvsZ->GetBinError(bL)>vmeanE ||*/ statZB[bL]<minBEntryZ) continue;
374 errL = currHDVvsZ->GetBinError(bL);
375 valL = currHDVvsZ->GetBinContent(bL);
378 for (bR=ib1;bR>=ib0;bR--) {
379 if (/*currHDVvsZ->GetBinError(bR)>vmeanE ||*/ statZB[bR]<minBEntryZ) continue;
380 errR = currHDVvsZ->GetBinError(bR);
381 valR = currHDVvsZ->GetBinContent(bR);
385 for (int ib=bL-1;ib>=ib0;ib--) {
386 double err = errL + (vmeanE-errL)*float(ib-ib0)/(bL-ib0+1);
387 double val = vmean + (valL-vmean)*float(ib-ib0)/(bL-ib0+1);
388 currHDVvsZ->SetBinError(ib,err);
389 currHDVvsZ->SetBinContent(ib,val);//currHDVvsZ->GetBinContent(ib+1));
391 for (int ib=bR+1;ib<=ib1;ib++) {
392 double err = errR + (vmeanE-errR)*float(ib1-ib)/(ib1-bR+1);
393 double val = vmean + (valR-vmean)*float(ib1-ib)/(ib1-bR+1);
394 currHDVvsZ->SetBinError(ib,err);
395 currHDVvsZ->SetBinContent(ib,val);//currHDVvsZ->GetBinContent(ib-1));
398 for (int iord=0;iord<=maxVZOrd;iord++) {
399 currHDVvsZ->Fit(fitVvsZs[iord],"q0");
400 chi2s[iord] = fitVvsZs[iord]->GetChisquare();
401 int ndf = fitVvsZs[iord]->GetNDF();
402 if (ndf>0) chi2s[iord] /= ndf;
407 double bestChi2 = 9999;
408 for (int iOrd=0;iOrd<=maxVZOrd;iOrd++) {
409 if (chi2s[iOrd]<bestChi2-0.5) {bestChi2 = chi2s[iOrd]; bestOrd = iOrd;}
411 TF1* fitVvsZ = fitVvsZs[bestOrd];
412 currHDVvsZ->Fit(fitVvsZ,"q");
414 // extract mean correction neglecting the constant term
416 double freePar = 0;//fitVvsZ->GetParameter(1);
417 for (int ib=1;ib<=nb;ib++) {
418 if (statZB[ib]<1) continue;
419 vmean += statZB[ib]*(fitVvsZ->Eval(currHDVvsZ->GetBinCenter(ib)) - freePar); // account only Z dependent part
423 if (!resVDCorrZ[currSDDSide]) {
424 resVDCorrZ[currSDDSide] = new TH1F(Form("VDcorrvsZ%d",currSDDSide),Form("mean VDrift correction vs Z, side %d",currSDDSide),kNSDD,kSDDMin-0.5,kSDDMax+0.5);
425 resVDCorrZ[currSDDSide]->SetMarkerColor(2+currSDDSide);
426 resVDCorrZ[currSDDSide]->SetMarkerStyle(20+4*currSDDSide);
429 // if (currSDDSide) vmean = -vmean;
430 resVDCorrZ[currSDDSide]->SetBinContent(currMod+1, norm>0 ? vmean/norm : 0);
435 //_________________________________________________________________________
436 void PrepareModuleHistos(Int_t mdID, Int_t side, Bool_t vsX)
438 // retrieve QA histos, if needed, convert to TH1
439 SetModuleID(mdID, side);
441 currHDXvsX = GetQAHisto(kDxX);
442 currHTDvsX = GetQAHisto(kTDX);
445 currHTDvsZ = GetQAHisto(kTDZ);
446 currHDXvsZ = GetQAHisto(kDxZ);
450 //_________________________________________________________________________
451 void SetModuleID(Int_t mdID,Int_t side)
453 // set current module ID,side
456 currMod = currSDDId - kSDDMin;
459 //_________________________________________________________________________
460 TProfile* GetQAHisto(Int_t qaType)
462 // retrieve needed QA histo
463 if (!qaHistos) {printf("QA Histos collection is not set\n"); exit(1);}
465 if (currSDDId<kSDDMin || currSDDId>kSDDMax || currSDDSide<0 || currSDDSide>2)
466 {printf("Illegal SDD module ID/Side: %d/%d\n",currSDDId,currSDDSide); exit(1);}
468 if (qaType<0 || qaType>=kNQATypes)
469 {printf("Illegal QA Histo Type %d\n",qaType); exit(1);}
471 int trueSide = userLeftRightFromTASK ? (1-currSDDSide) : currSDDSide;
472 const char* hname = Form("%s%d_%s",kQAHistoNames[qaType],currSDDId,kLRNames[trueSide]);
473 TH1* h = (TH1*) qaHistos->FindObject(hname);
474 if (!h) {printf("Did not find histo %s in QA output collection\n",hname); exit(1);}
476 if (h->InheritsFrom(TH2::Class())) h = H2Profile((TH2*)h); // make sure it is not TH2 histo
478 if ( (qaType==kDxX || qaType==kTDX) && userRebinX>1) SafeRebin((TProfile*)h,userRebinX,kTRUE);
479 if ( (qaType==kDxZ || qaType==kTDZ) && userRebinZ>1) SafeRebin((TProfile*)h,userRebinZ,kFALSE);
484 //_________________________________________________________________________
485 TH1* H2ProfileAsTH1(TH2* h2)
487 // extract profile as TH1 histo
488 TProfile* prf = h2->ProfileX();
489 TH1* prof = ProfileAsTH1(prf, "_profH1");
494 //_________________________________________________________________________
495 TH1* ProfileAsTH1(TProfile* prf, const char* addName)
497 // convert profile to TH1 histo
498 TString nm = prf->GetName(); nm += addName;
499 TAxis* xax = prf->GetXaxis();
501 const TArrayD *bins = xax->GetXbins();
502 if (bins->GetSize() == 0) {
503 prof = new TH1F(nm.Data(),nm.Data(),prf->GetNbinsX(),xax->GetXmin(),xax->GetXmax());
506 prof = new TH1F(nm.Data(),nm.Data(),xax->GetNbins(),bins->GetArray());
508 for (int i=1;i<=prof->GetNbinsX();i++) {
509 prof->SetBinContent(i, prf->GetBinContent(i));
510 prof->SetBinError(i, prf->GetBinError(i));
515 //_________________________________________________________________________
516 TProfile* H2Profile(TH2* h2)
519 TString nm = h2->GetName(); nm += "_prof";
520 TProfile* prf = h2->ProfileX(nm.Data());
524 //_________________________________________________________________________________
525 TH1* CleanProfile(TProfile* profR)
527 // clean profile copy
528 TH1* prof = FitDXEdges(profR); // cure edges
530 int ib0 = profR->FindBin(1), ib1 = profR->FindBin(sddSeg.Dx()-1);
531 int nbn = profR->GetNbinsX();
532 int nSkip = skipDXEdge/profR->GetBinWidth(nbn/2);
533 int nMean = wDXEdge/profR->GetBinWidth(nbn/2);
535 // get mean occupancy skipping first and last nSkip occupied bins
537 int nbCntL=0, nbCntR=0,nskipL=0,nskipR=0;
539 for (int ib=ib0;ib<=ib1 && nbCntL<nMean;ib++) { // skipping left, get average of good bins
540 double vl = profR->GetBinEntries(ib);
541 if (vl<minBEntryX) continue;
542 if (nskipL<nSkip) {nskipL++; continue;}
546 if (nbCntL<1) return 0;
549 for (int ib=ib1+1;ib>ib0 && nbCntR<nMean;ib--) { // skipping right
550 double vl = profR->GetBinEntries(ib);
551 if (vl<minBEntryX) continue;
552 if (nskipR<nSkip) {nskipR++; continue;}
557 if (nbCntR<1) return 0;
560 prof->GetXaxis()->SetRange(ib0,ib1);
561 prof->Smooth(smoothLevel,"r");
562 prof->GetXaxis()->SetRange(1,nbn);
563 // kill bins with small statistics, from left and right
564 for (int ib=1;ib<ib1;ib++) if (profR->GetBinEntries(ib)<threshClean*smL) {prof->SetBinContent(ib,0);prof->SetBinError(ib,0);} else break;
565 for (int ib=ib1;ib>ib0;ib--) if (profR->GetBinEntries(ib)<threshClean*smR) {prof->SetBinContent(ib,0);prof->SetBinError(ib,0);} else break;
571 //_________________________________________________________________________________
572 TH1* Resid2Vdrift(TH1* res)
574 // extract Vdrift profile from residuals profile
577 TString nm = res->GetName(); nm += "_vd";
578 TH1* vd = (TH1*) res->Clone(nm.Data());
580 int nb = vd->GetNbinsX();
582 int nbcnt=0, nbfl = wVDEdge/vd->GetBinWidth(nb/2);
583 if (nbfl<nVDEdge) nbfl=nVDEdge;
585 double vAv=0,eAv=0; // average vdrift / v0
586 for (int i=1;i<=nb;i++) {
587 double deriv = useSplieForR2V ? spl.Derivative(vd->GetBinCenter(i)) : (res->GetBinContent(i)-res->GetBinContent(i-1))/res->GetBinWidth(i);
588 double v = TMath::Abs(deriv-1)>1e-6 ? 1./(1-deriv) : 1.;
589 vd->SetBinContent(i, v);
590 double err2 = 0; // set relative error
592 for (int i1=i-1;i1<=i+1;i1++) {
593 if (i1<1 || i1>nb) continue;
594 double err = res->GetBinError(i);
595 if (err<1e-9) err = 1e4;
600 if (err2>0 && v>0 && v<3) {vAv += v/err2; eAv += 1./err2;}
601 err2 = TMath::Sqrt(err2)/res->GetBinWidth(nb/2);
602 vd->SetBinError(i, err2 );
605 vAv = eAv>0 ? vAv/eAv : 1.0;
606 // printf("mean V = %f in %d bins\n",vAv,nbcnt);
607 // cure anomalous bins
608 // for (int i=1;i<=nb;i++) if (vd->GetBinContent(i)<0.1) vd->SetBinContent(i,vAv);
611 for (ib=1;ib<nb;ib++) { // detect up to which bin the left tail is unstable
612 double x = vd->GetBinCenter(ib);
613 if (x<0 || res->GetBinError(ib)<1e-9) {vd->SetBinContent(ib,1); vd->SetBinError(ib,0); continue;}
614 double vl = vd->GetBinContent(ib);
615 if (TMath::Abs(vl-vAv)<0.2) break;
616 vd->SetBinContent(ib,1);
617 vd->SetBinError(ib,0);
620 double vL=0,eL=0,vR=0,eR=0;
621 // get the mean of leftmost stable vdrift
624 for (int i=ib;i<=nb && nbcnt<nbfl;i++) {
625 double berr = vd->GetBinError(i);
626 if (berr<1e-9 || berr>0.5 || TMath::Abs(vd->GetBinContent(i)-vAv)>0.2) continue;
627 vL += vd->GetBinContent(i)/berr/berr;
632 vL = eL>0 ? vL/eL : vAv;
633 //printf("VLeft: %f, in %d bins (%d %d)\n",vL,nbcnt,ib,lastCounted);
634 for (int i=1;i<=ib;i++) vd->SetBinContent(i, vL);
635 // for safety check if there are no outliers in first 3 "stable" bins
636 for (int i=ib+1;i<ib+3;i++) if ( vd->GetBinError(i)<1e-9 || TMath::Abs(vd->GetBinContent(i)-vL)>0.2) vd->SetBinContent(i, vL);
637 vd->SetBinContent(vd->FindBin(1),1.); // no correction at t=0 !!
639 double lmax = sddSeg.Dx()-1;
640 for (ib=nb+1;ib--;) { // detect up to which bin the right tail is unstable
641 double x = vd->GetBinCenter(ib);
642 if (x>=lmax || res->GetBinError(ib)<1e-9) {vd->SetBinContent(ib,1);vd->SetBinError(ib,0); continue;}
643 if (TMath::Abs(vd->GetBinContent(ib)-vAv)<0.4) break;
644 vd->SetBinContent(ib,vAv);
645 vd->SetBinError(ib,0);
649 for (int i=ib;i>=1 && nbcnt<nbfl;i--) {
650 double berr = vd->GetBinError(i);
651 if (berr<1e-9 || berr>0.5 || TMath::Abs(vd->GetBinContent(i)-vAv)>0.4) continue;
652 vR += vd->GetBinContent(i)/berr/berr;
657 vR = eR>0 ? vR/eR : vAv;
658 // printf("VRight: %f, in %d bins (%d %d)\n",vR,nbcnt,lastCounted,ib);
659 for (int i=ib;i<=nb;i++) vd->SetBinContent(i, vR);
660 // for safety check if there are no outliers in first 3 "stable bins
661 for (int i=(lastCounted+ib)/2;i<ib;i++)
662 if ( vd->GetBinError(i)<1e-9 || (vd->GetBinError(i)>0.3 && TMath::Abs(vd->GetBinContent(i)-vR)>0.2)) vd->SetBinContent(i, vR);
663 // fit the empty bins on the right
664 // vd->Fit(fitP1,"+","",vd->GetBinCenter(ib-nbfl),vd->GetBinCenter(nb));
665 //for (int i=ib;i<=nb;i++) vd->SetBinContent(i, fitP1->Eval(vd->GetBinCenter(i)));
671 //_________________________________________________________________________________
672 TH1* Vdrift2Resid(TH1* vd)
674 // convert Vdrift profile to residuals profile (final correction map)
676 TString nm = vd->GetName(); nm += "_vd2res";
677 TH1* res = (TH1*) vd->Clone(nm.Data());
679 if (userDummyCorrMap) return res;
680 int ib0 = res->FindBin(1);
681 int ib1 = res->FindBin(sddSeg.Dx()-1);
684 for (lastB=ib1+1;lastB--;) if (vd->GetBinError(lastB)>1e-9) break; // find last good bin
685 double lastX = vd->GetBinCenter(lastB);
688 if (lastX>sddSeg.Dx()) lastX = sddSeg.Dx();
689 lastB = vd->FindBin(lastX);
691 // 1st iteration : estimate correction at max Xdrift
692 for (int i=ib0;i<=lastB;i++) {
693 double dx = res->GetBinWidth(i);
694 double v = vd->GetBinContent(i);
695 resv += dx*(1.-1./v);
697 double vcorr = (resv)/lastX;
699 // 2nd iteration : create new residuals forcing them to be 0 at Xdrift=0 and maxXDrift
700 resv = res->GetBinWidth(ib0)*vcorr;
701 for (int i=ib0;i<=lastB;i++) {
702 double dx = res->GetBinWidth(i);
703 double v = vd->GetBinContent(i);
704 resv += dx*(1.-1./v - vcorr);
705 res->SetBinContent(i, resv);
711 //__________________________________________________________________________________
714 // check mean residuals before and after correction
715 const double kFOffs = 0.05; // skip edges
716 static TF1* fitP1 = new TF1("fitP1","pol1",-5000,40000);
717 static TF1* fitP0 = new TF1("fitP0","pol0",-5000,40000);
719 if (currHDXvsX->GetEntries()<minBEntryX*currHDXvsX->GetNbinsX()) {
720 if (currHCorrMapX) currHCorrMapX->Reset();
725 int b0 = currHDXvsXCorr->FindBin(1),b1 = currHDXvsXCorr->FindBin(sddSeg.Dx()-1);
727 for (int ib=b0;ib<b1;ib++) if (currHTDvsX->GetBinEntries(ib)>=minBEntryX) nb++;
728 currHDXvsX->Fit(fitP0,"q0N","");
729 double offsRaw = fitP0->GetParameter(0);
730 currHDXvsXCorr->Fit(fitP0,"q0N","");
731 double offsAMap = fitP0->GetParameter(0);
734 currGrDxx = new TGraphErrors(nb); // residual vs Xdrift
735 currGrTDx = new TGraphErrors(nb); // residual vs tdrift
736 currGrTXCorr = new TGraphErrors(nb); // Xdrift vs tdrift
737 double tmin = 1e6, tmax = -1e6, xmin = 1e6, xmax = -1e6;
739 for (int i=b0;i<b1;i++) {
740 if (currHTDvsX->GetBinEntries(i)<minBEntryX) continue;
741 double t = currHTDvsX->GetBinContent(i);
742 double x = currHDXvsXCorr->GetBinCenter(i);
743 if (tmin>t) tmin = t;
744 if (tmax<t) tmax = t;
745 if (xmin>x) xmin = x;
746 if (xmax<x) xmax = x;
747 currGrDxx->SetPoint(ip,x,currHDXvsXCorr->GetBinContent(i));
748 currGrDxx->SetPointError(ip,currHDXvsXCorr->GetBinWidth(i),currHDXvsXCorr->GetBinError(i));
750 currGrTDx->SetPoint(ip,t, currHDXvsXCorr->GetBinContent(i));
751 currGrTDx->SetPointError(ip, currHTDvsX->GetBinError(i), currHDXvsXCorr->GetBinError(i));
753 currGrTXCorr->SetPoint(ip,t, currHTDvsX->GetBinCenter(i));
754 currGrTXCorr->SetPointError(ip,currHTDvsX->GetBinError(i), currHTDvsX->GetBinWidth(i));
758 double del = tmax-tmin;
765 fitP1->SetParameters(0,0);
766 currGrDxx->Fit(fitP1,"q","",xmin, xmax);
767 double offs = fitP1->GetParameter(0); // offset of correction line at Xdrift=0
769 // printf("Fitting VD correction in the range %.1f:%.1f\n",tmin,tmax);
770 fitP1->SetParameters(0,0);
771 currGrTDx->Fit(fitP1,"q","",tmin,tmax);
772 double vcor = fitP1->GetParameter(1);
774 fitP1->SetParameters(0,0);
775 currGrTXCorr->Fit(fitP1,"q","",tmin,tmax);
776 double vav = fitP1->GetParameter(1);
779 if (!resOffsDXraw[currSDDSide]) {
780 resOffsDXraw[currSDDSide] = new TH1F(Form("OffsRaw%d",currSDDSide),Form("DX Raw Offset, side %d",currSDDSide),kNSDD,kSDDMin-0.5,kSDDMax+0.5);
781 resOffsDXraw[currSDDSide]->SetMarkerColor(2+currSDDSide);
782 resOffsDXraw[currSDDSide]->SetMarkerStyle(20+4*currSDDSide);
785 if (!resOffsDXAMap[currSDDSide]) {
786 resOffsDXAMap[currSDDSide] = new TH1F(Form("OffsAMap%d",currSDDSide),Form("DX Offset after Map corr. Mean, side %d",currSDDSide),kNSDD,kSDDMin-0.5,kSDDMax+0.5);
787 resOffsDXAMap[currSDDSide]->SetMarkerColor(3+currSDDSide);
788 resOffsDXAMap[currSDDSide]->SetMarkerStyle(20+4*currSDDSide);
791 if (!resOffsDX[currSDDSide]) {
792 resOffsDX[currSDDSide] = new TH1F(Form("Offs%d",currSDDSide),Form("DX Offset at TD=0 after Map corr., side %d",currSDDSide),kNSDD,kSDDMin-0.5,kSDDMax+0.5);
793 resOffsDX[currSDDSide]->SetMarkerColor(2+currSDDSide);
794 resOffsDX[currSDDSide]->SetMarkerStyle(20+4*currSDDSide);
797 if (!resVDCorr[currSDDSide]) {
798 resVDCorr[currSDDSide] = new TH1F(Form("VDcorr%d",currSDDSide),Form("VDrift correction, side %d",currSDDSide),kNSDD,kSDDMin-0.5,kSDDMax+0.5);
799 resVDCorr[currSDDSide]->SetMarkerColor(2+currSDDSide);
800 resVDCorr[currSDDSide]->SetMarkerStyle(20+4*currSDDSide);
803 if (!resVDMean[currSDDSide]) {
804 resVDMean[currSDDSide] = new TH1F(Form("VDmean%d",currSDDSide),Form("VDrift mean, side %d",currSDDSide),kNSDD,kSDDMin-0.5,kSDDMax+0.5);
805 resVDMean[currSDDSide]->SetMarkerColor(2+currSDDSide);
806 resVDMean[currSDDSide]->SetMarkerStyle(20+4*currSDDSide);
809 resOffsDXraw[currSDDSide]->SetBinContent(currMod+1, offsRaw);
810 resOffsDXAMap[currSDDSide]->SetBinContent(currMod+1, offsAMap);
811 resOffsDX[currSDDSide]->SetBinContent(currMod+1, offs);
812 resVDCorr[currSDDSide]->SetBinContent(currMod+1, vcor);
813 resVDMean[currSDDSide]->SetBinContent(currMod+1, vav);
817 //__________________________________________________________________________
818 void CalcDXCorrections()
820 // estimate time0 and alignment correction for the whole module
822 resT0Corr = new TH1F("T0Corr","T0 Correction",kNSDD,kSDDMin-0.5,kSDDMax+0.5);
823 resT0Corr->SetMarkerColor(2);
824 resT0Corr->SetMarkerStyle(20);
828 resXLocCorr = new TH1F("XLocCorr","XLoc Correction",kNSDD,kSDDMin-0.5,kSDDMax+0.5);
829 resXLocCorr->SetMarkerColor(2);
830 resXLocCorr->SetMarkerStyle(20);
833 if (!resVDMean[0] || !resVDMean[1]) return;
834 if (!resOffsDX[0] || !resOffsDX[1]) return;
835 double vL = resVDMean[0]->GetBinContent(currMod+1); // true mean VL
836 double vR = resVDMean[1]->GetBinContent(currMod+1); // true mean VR
837 double offsL = resOffsDX[0]->GetBinContent(currMod+1);
838 double offsR = resOffsDX[1]->GetBinContent(currMod+1);
840 double vsum=0,t0Corr=0,xlCorr=0;
841 if (vL>1 && vR>1) { // both sides available
843 t0Corr = -(offsL+offsR)/vsum;
844 xlCorr = -(offsL*vR - offsR*vL)/vsum;
847 else if (vL>1) t0Corr = -offsL/vL; // only one side is available
848 else if (vR>1) t0Corr = -offsR/vR;
850 else if (vL>1) xlCorr = -offsL; // only one side is available
851 else if (vR>1) xlCorr = offsR;
853 if (userDummyt0Corr) t0Corr = 0;
854 if (userDummyxlCorr) xlCorr = 0;
855 // printf("SDD%d VL:%f VR:%f offsL:%+f offsR:%+f dT:%+f dX:%+f\n",currSDDId, vL,vR, offsL,offsR, t0Corr,xlCorr);
856 resT0Corr->SetBinContent(currMod+1, t0Corr); // T0 correction
857 resXLocCorr->SetBinContent(currMod+1, xlCorr); // X alignment correction
859 double addMap[2]={0,0};
860 Bool_t redoMaps = kFALSE;
862 if (forceT0CorrTo0) { // T0 correction was forced to be 0, attribute this to map
863 addMap[0] -= vL*t0Corr;
864 addMap[1] -= vR*t0Corr;
867 if (forceRPhiCorrTo0) { // alignment correction was forced to be 0, attribute this to map
869 addMap[1] -= -xlCorr;
874 for (int ix=0;ix<2;ix++) {
875 TH1* map = (TH1*)procHistos.At( GetStoreID(kSCorrMapX, currSDDId, ix) );
876 TH1* mapc = (TH1*)procHistos.At( GetStoreID(kSXvsXCorr, currSDDId, ix) );
877 if (!map || !mapc) continue;
878 int ib0 = map->FindBin(1);
879 int ib1 = map->FindBin(sddSeg.Dx()-1);
880 for (int ib=ib0+1;ib<ib1;ib++) {
881 map->AddBinContent(ib, addMap[ix]);
882 mapc->AddBinContent(ib, -addMap[ix]);
890 //______________________________________________________________
891 Int_t GetStoreID(int type, int imd,int side)
893 // entry of the histo/graph of type in the procHistos array
895 if (imd<0) imd = currSDDId;
896 if (side<0) side = currSDDSide;
897 if (type<0||type>=kNStore || imd<kSDDMin || imd>kSDDMax || side<0 || side>1) {
898 printf("Wrong object requested: type: %d, Mod:%d/%d\n",type,imd,side);
901 return (2*(imd-kSDDMin)+side)*kNStore + type;
904 //______________________________________________________________
907 // clean "current" objects from last event
923 //______________________________________________________________
926 // store "current" objects in procHistos
927 if (currHDXvsX) procHistos.AddAtAndExpand(currHDXvsX, GetStoreID(kSXvsX));
928 if (currHDXvsZ) procHistos.AddAtAndExpand(currHDXvsZ, GetStoreID(kSXvsZ));
929 if (currHTDvsX) procHistos.AddAtAndExpand(currHTDvsX, GetStoreID(kSTDvsX));
930 if (currHTDvsZ) procHistos.AddAtAndExpand(currHTDvsZ, GetStoreID(kSTDvsZ));
931 if (currHDXvsXclean) procHistos.AddAtAndExpand(currHDXvsXclean, GetStoreID(kSDXvsXclean));
932 if (currHVDvsX) procHistos.AddAtAndExpand(currHVDvsX, GetStoreID(kSVDvsX));
933 if (currHCorrMapX) procHistos.AddAtAndExpand(currHCorrMapX, GetStoreID(kSCorrMapX));
934 if (currHDXvsXCorr) procHistos.AddAtAndExpand(currHDXvsXCorr, GetStoreID(kSXvsXCorr));
935 if (currHDXvsZCorrMT0) procHistos.AddAtAndExpand(currHDXvsZCorrMT0, GetStoreID(kSXvsZCorr));
936 if (currHDVvsZ) procHistos.AddAtAndExpand(currHDVvsZ, GetStoreID(kSDVvsZ));
937 if (currHDVvsZOrig) procHistos.AddAtAndExpand(currHDVvsZOrig, GetStoreID(kSDVvsZOrig));
938 if (currGrDxx) procHistos.AddAtAndExpand(currGrDxx, GetStoreID(kSGrDxx));
939 if (currGrTDx) procHistos.AddAtAndExpand(currGrTDx, GetStoreID(kSGrTDx));
940 if (currGrTXCorr) procHistos.AddAtAndExpand(currGrTXCorr, GetStoreID(kSGrTXCorr));
944 //_________________________________________________________________________________
945 TObjArray* CreateCorrMaps()
947 // create correction maps for all modules
948 printf("Creating correction maps (update %s)\n",pathSDDCorrMapOld.Data());
949 TObjArray *dest = new TObjArray(2*kNSDD);
950 TObjArray* update = 0;
951 if (!pathSDDCorrMapOld.IsNull() && !LoadSDDCorrMap(pathSDDCorrMapOld,update)) {
952 printf("The update of correction map was requested but the source %s is not found\n",pathSDDCorrMapOld.Data());
957 AliITSCorrMap1DSDD *updMap = 0;
958 for (int imd=kSDDMin;imd<=kSDDMax;imd++) {
959 for (int side=0;side<2;side++) {
960 TH1* mph = (TH1*)procHistos.At( GetStoreID(kSCorrMapX,imd,side) );
961 //if (!mph) printf("Correction map is missing for module %d/%d\n",imd,side);
962 if (update) updMap = (AliITSCorrMap1DSDD*)update->At(2*(imd-kSDDMin) + side);
963 AliITSCorrMap1DSDD* mp = CreateCorrMap(mph,imd,side, updMap);
964 dest->AddAtAndExpand(mp, 2*(imd-kSDDMin) + side);
971 //_________________________________________________________________________________
972 AliITSCorrMap1DSDD* CreateCorrMap(TH1* mapHisto, int imd, int side, AliITSCorrMap1DSDD* updateMap)
974 // create or update correction map from histo
975 int nbCorr = 1, nbOld = 0;
978 b0 = mapHisto->FindBin(1);
979 b1 = mapHisto->FindBin(sddSeg.Dx()-1);
982 AliITSCorrMap1DSDD* mpCorr = 0;
984 // check if the updateMap is meaningful
985 if (updateMap && updateMap->GetNBinsDrift()>2 && nbCorr>1) {
987 TSpline3 spl(mapHisto);
988 nbOld = updateMap->GetNBinsDrift();
989 double dx = sddSeg.Dx()/nbOld;
990 for (int ip=0;ip<nbOld;ip++) {
991 double x = dx*(0.5+ip);
992 updateMap->SetCellContent(0,ip,updateMap->GetCellContent(0,ip)-spl.Eval(x));
998 mpCorr = new AliITSCorrMap1DSDD(Form("DriftTimeMap_%d_%d",imd,side),nbCorr);
999 if (side==0) mpCorr->SetInversionBit(); // !!! left side should return correction*-1
1000 if (mapHisto) for (int ib=b0;ib<=b1;ib++) mpCorr->SetCellContent(0,ib-b0,-mapHisto->GetBinContent(ib));
1006 //_________________________________________________________________________________
1007 TObjArray* UpdateSDDVDrift()
1009 // retrieve SDD VDrift object and update it
1010 if (!vdarrayOld && !LoadSDDVDrift(pathSDDVDriftOld,vdarrayOld)) return 0;
1011 TObjArray *vdarrayUp = new TObjArray(2*kNSDD);
1013 for (int imd=kSDDMin;imd<=kSDDMax;imd++) {
1014 for (int side=0;side<2;side++) {
1015 int iad = 2*(imd-kSDDMin)+side;
1016 AliITSDriftSpeedArraySDD* drarr = (AliITSDriftSpeedArraySDD*) vdarrayOld->At( iad );
1017 AliITSDriftSpeedArraySDD* drarrUp = new AliITSDriftSpeedArraySDD();
1018 AliITSDriftSpeedSDD* vOr = drarr->GetDriftSpeedObject(0);
1019 AliITSDriftSpeedSDD* vUp = new AliITSDriftSpeedSDD(*vOr);
1020 drarrUp->AddDriftSpeed(vUp);
1021 vdarrayUp->AddAt(drarrUp, iad);
1022 UpdateSDDVDrift(drarrUp, imd, side);
1026 sddVDriftUpdOK = kTRUE;
1031 //_________________________________________________________________________________
1032 void UpdateSDDVDrift(AliITSDriftSpeedArraySDD* vdArr, int imd, int side)
1034 // update vdrift vs anode in the object
1035 AliITSDriftSpeedSDD* ds;
1036 if (!vdArr || !(ds=vdArr->GetDriftSpeedObject(0))) {printf("No VDrift object for module %d/%d\n",imd,side); exit(1);}
1037 TH1* vdh = (TH1*)procHistos.At( GetStoreID(kSDVvsZ,imd,side) );
1039 //printf("VDrift vs Z correction is not processed for module %d/%d\n",imd,side);
1042 TF1* fp = vdh->GetFunction("fitVvsZ");
1043 if (!fp) {printf("VDrift vs Z correction fit is missing SDD%d/%d\n",imd,side); return;}
1045 int ord = (int)fp->GetParameter(0); // 1st param is the order of poly
1046 int ordOld = ds->GetDegreeofPoly();
1047 if (ord>ordOld) ds->SetDegreeofPoly(ord);
1048 for (int ip=0;ip<ord+1;ip++) { // don't store offset (par[1])
1049 double par = ds->GetDriftSpeedParameter(ip) - fp->GetParameter(ip+1);
1050 ds->SetDriftSpeedParameter(ip, par);
1055 //_________________________________________________________________________________
1056 AliITSresponseSDD* UpdateSDDResponse(Bool_t t0, Bool_t vdrift)
1058 // retrieve RespSDD object and update it
1059 AliITSresponseSDD* resp = 0;
1060 if (!LoadSDDResponse(pathSDDRespOld, resp)) return 0;
1061 UpdateSDDResponse(resp, t0, vdrift);
1062 sddRespUpdT0OK = t0;
1063 sddRespUpdVDOK = vdrift;
1068 //_________________________________________________________________________________
1069 void UpdateSDDResponse(AliITSresponseSDD *resp, Bool_t t0, Bool_t vdrift)
1071 // update the map with extracted values
1072 printf("Updating RespSDD object: T0:%s VDrift:%s\n",t0?"ON":"OFF",vdrift?"ON":"OFF");
1074 if (t0 && !resT0Corr)
1075 {printf("T0 update is requested but corrections were not processed"); exit(1);}
1076 if (vdrift && !(resVDCorr[0] && resVDCorr[1]))
1077 {printf("VDrift update is requested but corrections were not processed"); exit(1);}
1079 for (int imd=kSDDMin;imd<=kSDDMax;imd++) {
1080 if (t0 && !forceT0CorrTo0) resp->SetModuleTimeZero(imd, resp->GetTimeZero(imd) - resT0Corr->GetBinContent(imd-kSDDMin+1));
1082 for (int ix=0;ix<2;ix++) {
1083 double vdZ = sddVDriftUpdOK&&resVDCorrZ[ix] ? resVDCorrZ[ix]->GetBinContent(imd-kSDDMin+1) : 0; // contribution from DXvsZ correction
1084 double vdX = resVDCorr[ix]->GetBinContent(imd-kSDDMin+1); // contribution from DXvsX correction
1085 resp->SetDeltaVDrift(imd, resp->GetDeltaVDrift(imd,ix) - (vdX-vdZ), ix);
1092 //___________________________________________________________________
1093 double GetVOld(double z)
1095 // return VDrift assumed in reconstruction
1096 if (!vdarrayOld && !LoadSDDVDrift(pathSDDVDriftOld,vdarrayOld)) return 0;
1097 AliITSDriftSpeedArraySDD* drarr = (AliITSDriftSpeedArraySDD*) vdarrayOld->At( 2*currMod + currSDDSide);
1098 float anode = sddSeg.GetAnodeFromLocal( currSDDSide==0 ? 1.:-1. ,z*1e-4);
1099 double v = drarr->GetDriftSpeed(0, anode);
1103 //___________________________________________________________________
1104 double ftVdZ(double *x, double *par)
1106 // function to fit the vdrift dependence on Z
1108 // convert Z to anode
1110 double ian = (z/sddSeg.Dz() + 0.5);
1111 if (ian<0) ian = 0.;
1112 else if (ian>1) ian = 1.;
1113 ian *= sddSeg.NpzHalf();
1115 int ord = int(par[0]);
1116 double v = par[ord+1];
1117 for (int i=ord;i--;) v = par[i+1]+ian*v;
1121 //________________________________________________________________________________________________________
1122 Bool_t LoadSDDVDrift(TString& path, TObjArray *&arr)
1124 // load VDrift object
1125 if (path.IsNull()) return kFALSE;
1126 printf("Loading SDD VDrift from %s\n",path.Data());
1128 AliCDBEntry *entry = 0;
1132 if (path.BeginsWith("path: ")) { // must load from OCDB
1133 entry = GetCDBEntry(path.Data());
1135 arr = (TObjArray*) entry->GetObject();
1136 entry->SetObject(NULL);
1137 entry->SetOwner(kTRUE);
1141 if (gSystem->AccessPathName(path.Data())) break;
1142 TFile* precf = TFile::Open(path.Data());
1143 if (precf->FindKey("TObjArray")) arr = (TObjArray*)precf->Get("TObjArray");
1144 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
1145 arr = (TObjArray*) entry->GetObject();
1146 if (arr && arr->InheritsFrom(TObjArray::Class())) entry->SetObject(NULL);
1148 entry->SetObject(NULL);
1149 entry->SetOwner(kTRUE);
1158 if (!arr) {printf("Failed to load SDD vdrift from %s\n",path.Data()); return kFALSE;}
1159 arr->SetOwner(kTRUE);
1163 //________________________________________________________________________________________________________
1164 Bool_t LoadSDDResponse(TString& path, AliITSresponseSDD *&resp)
1166 // load SDD response
1167 if (path.IsNull()) return kFALSE;
1168 printf("Loading SDD response from %s\n",path.Data());
1170 AliCDBEntry *entry = 0;
1175 if (path.BeginsWith("path: ")) { // must load from OCDB
1176 entry = GetCDBEntry(path.Data());
1178 resp = (AliITSresponseSDD*) entry->GetObject();
1179 entry->SetObject(NULL);
1180 entry->SetOwner(kTRUE);
1184 if (gSystem->AccessPathName(path.Data())) break;
1185 TFile* precf = TFile::Open(path.Data());
1186 if (precf->FindKey("AliITSresponseSDD")) resp = (AliITSresponseSDD*)precf->Get("AliITSresponseSDD");
1187 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
1188 resp = (AliITSresponseSDD*) entry->GetObject();
1189 if (resp && resp->InheritsFrom(AliITSresponseSDD::Class())) entry->SetObject(NULL);
1191 entry->SetObject(NULL);
1192 entry->SetOwner(kTRUE);
1201 if (!resp) {printf("Error: Failed to load SDD response from %s\n",path.Data()); return kFALSE;}
1205 //________________________________________________________________________________________________________
1206 Bool_t LoadSDDCorrMap(TString& path, TObjArray *&maps)
1208 // Load SDD correction map
1210 if (path.IsNull()) return kFALSE;
1211 printf("Loading SDD Correction Maps from %s\n",path.Data());
1213 AliCDBEntry *entry = 0;
1217 if (path.BeginsWith("path: ")) { // must load from OCDB
1218 entry = GetCDBEntry(path.Data());
1220 maps = (TObjArray*) entry->GetObject();
1221 entry->SetObject(NULL);
1222 entry->SetOwner(kTRUE);
1226 if (gSystem->AccessPathName(path.Data())) break;
1227 TFile* precf = TFile::Open(path.Data());
1228 if (precf->FindKey("TObjArray")) maps = (TObjArray*)precf->Get("TObjArray");
1229 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
1230 maps = (TObjArray*) entry->GetObject();
1231 if (maps && maps->InheritsFrom(TObjArray::Class())) entry->SetObject(NULL);
1233 entry->SetObject(NULL);
1234 entry->SetOwner(kTRUE);
1243 if (!maps) {printf("Failed to load SDD Correction Map from %s\n",path.Data()); return kFALSE;}
1248 //_______________________________________________________________________________________
1249 AliCDBEntry* GetCDBEntry(const char* path)
1251 // return object from the OCDB
1252 AliCDBEntry *entry = 0;
1253 printf("Loading object %s\n",path);
1254 AliCDBManager* man = AliCDBManager::Instance();
1255 AliCDBId* cdbId = AliCDBId::MakeFromString(path);
1257 printf("Failed to create cdbId\n");
1261 AliCDBStorage* stor = man->GetDefaultStorage();
1262 if (!stor && !man->GetRaw()) man->SetDefaultStorage("raw://");
1263 if (man->GetRaw()) man->SetRun(cdbId->GetFirstRun());
1265 TString tp = stor->GetType();
1266 if (tp.Contains("alien",TString::kIgnoreCase) && !gGrid) TGrid::Connect("alien:");
1268 entry = man->Get(cdbId->GetPath(),cdbId->GetFirstRun(),cdbId->GetVersion(),cdbId->GetSubVersion());
1269 // entry = man->Get( *cdbId );
1278 //_______________________________________________________________________________________
1279 Bool_t PlotHisto(TH1* h, Option_t* opt, int mrkStyle,int mrkCol, double mrkSize)
1281 const double kOffsH = 0.15;
1282 if (!h) return kFALSE;
1283 TString opts = opt; opts.ToLower();
1284 if (opts.Contains("p")) {
1285 h->SetMarkerStyle(mrkStyle);
1286 h->SetMarkerColor(mrkCol);
1287 h->SetMarkerSize(mrkSize);
1289 h->SetLineColor(mrkCol);
1292 h->SetMinimum(); h->SetMaximum();
1293 double hmn=h->GetMinimum(),hmx=h->GetMaximum(); // new histo min/max
1295 TH1* hbase = GetPadBaseHisto((TPad*)gPad); if (!hbase) return 0;
1296 double smn = hbase->GetMinimum(),smx = hbase->GetMaximum(); // base set min/max?
1297 hbase->SetMinimum(); hbase->SetMaximum();
1298 double omn = hbase->GetMinimum(),omx = hbase->GetMaximum(); // base real min max
1299 if (smn<omn && smx>omx) { // min/max for bas histo was set by hand: extract original min/max
1300 omx = (smn*kOffsH+smx*(1+kOffsH))/(1+2*kOffsH);
1301 omn = (smn-kOffsH*omx)/(1+kOffsH);
1303 if (hmn<omn) omn = hmn;
1304 if (hmx>omx) omx = hmx;
1305 double del = omx-omn;
1306 hbase->SetMinimum( omn - kOffsH*del );
1307 hbase->SetMaximum( omx + kOffsH*del );
1312 //_______________________________________________________________________________________
1313 void GetHistoMinMaxInRange(TH1* h, double &mn,double &mx)
1315 // compute min/max of histo in the user range
1318 int b0 = h->GetXaxis()->GetFirst(), b1 = h->GetXaxis()->GetLast();
1319 for (int i=b0;i<=b1;i++) {
1320 double e = h->GetBinError(i);
1321 if (TMath::Abs(e)<1e-9) continue;
1322 double v = h->GetBinContent(i);
1323 if (mn>v-e) mn = v-e;
1324 if (mx<v+e) mx = v+e;
1328 //_______________________________________________________________________________________
1329 void PlotReport(const char* psname)
1332 sddCanv = new TCanvas("sddCanv","sddCanv",700,900);
1334 gStyle->SetOptStat(0);
1335 gStyle->SetOptTitle(0);
1337 TString psnm1 = psname;
1338 if (psnm1.IsNull()) psnm1 = "sddQAreport.ps";
1339 TString psnm0 = psnm1 + "[";
1340 TString psnm2 = psnm1 + "]";
1341 sddCanv->Print(psnm0.Data());
1343 // mean corrections per module/side
1345 sddCanv->Divide(2,3);
1348 for (int ix=0;ix<2;ix++) { // mean residuals before/after correction
1349 sddCanv->cd(++cntPad);
1350 PlotHisto(resOffsDXraw[ix],"p" ,7,kBlack,0.5);
1351 PlotHisto(resOffsDX[ix] ,"p same",7,kRed ,1);
1352 AddPadLabel(Form("<#DeltaX> %s : Raw",ix?"Right":"Left"), 0.1,0.93,kBlack,0.05);
1353 AddPadLabel("After Map", 0.5,0.93,kRed,0.05);
1356 for (int ix=0;ix<2;ix++) { // mean residuals before/after correction
1357 sddCanv->cd(++cntPad);
1358 PlotHisto(resVDCorr[ix] ,"p" ,7,kBlack,1);
1359 PlotHisto(resVDCorrZ[ix],"p same",7,kRed ,1);
1360 AddPadLabel(Form("<#DeltaV> %s : from #DeltaX vs X",ix?"Right":"Left"), 0.1,0.93,kBlack,0.05);
1361 AddPadLabel("from #DeltaX vs Z", 0.6,0.93,kRed,0.05);
1364 sddCanv->cd(++cntPad);
1365 PlotHisto(resT0Corr,"p",7,kBlue,1);
1366 AddPadLabel("T0 Correction", 0.3,0.93,kRed,0.05);
1367 if (forceT0CorrTo0) AddPadLabel("Forced to 0 by transferring to maps", 0.15,0.88,kRed,0.05);
1369 sddCanv->cd(++cntPad);
1370 PlotHisto(resXLocCorr,"p",7,kBlack,1);
1371 AddPadLabel("#DeltaX Correction", 0.3,0.93,kRed,0.05);
1372 if (forceRPhiCorrTo0) AddPadLabel("Forced to 0 by transferring to maps", 0.15,0.88,kRed,0.05);
1375 sddCanv->Print(psnm1.Data());
1377 //-------------------------------------------------------------------
1381 int nModPerPage = 3;
1383 Bool_t saved = kFALSE;
1386 for (int imd=kSDDMin;imd<=kSDDMax;imd++) {
1387 if (cntPad>=2*nModPerPage*nRowPerMod) {
1389 if (imd!=kSDDMin) sddCanv->Print(psnm1.Data());
1391 sddCanv->Divide(2,nModPerPage*nRowPerMod);
1395 for (int ix=0;ix<2;ix++) {
1396 sddCanv->cd(++cntPad);
1397 TH1* hsddcl = (TH1*)procHistos.At(GetStoreID(kSDXvsXclean,imd,ix)); // raw residuals
1399 ib0 = hsddcl->FindBin(1);
1400 ib1 = hsddcl->FindBin(sddSeg.Dx()-1);
1401 hsddcl->GetXaxis()->SetRange(ib0,ib1);
1403 PlotHisto(hsddcl,"p",24,kGreen+2,0.5);
1405 hsdd = (TH1*)procHistos.At(GetStoreID(kSXvsX,imd,ix)); // raw residuals
1407 ib0 = hsdd->FindBin(1);
1408 ib1 = hsdd->FindBin(sddSeg.Dx()-1);
1409 hsdd->GetXaxis()->SetRange(ib0,ib1);
1411 PlotHisto(hsdd,"p same",20,kBlack,0.4);
1412 PlotHisto(hsddcl,"p sames",24,kGreen+2,0.5);
1415 hsdd = (TH1*)procHistos.At(GetStoreID(kSCorrMapX,imd,ix)); // map
1416 if (hsdd) hsdd->GetXaxis()->SetRange(ib0,ib1);
1417 PlotHisto(hsdd,"same",7,kRed,0.5);
1418 hsdd = (TH1*)procHistos.At(GetStoreID(kSXvsXCorr,imd,ix));
1419 if (hsdd) hsdd->GetXaxis()->SetRange(ib0,ib1);
1420 PlotHisto(hsdd,"histo same",7,kBlue,0.5);
1422 AddPadLabel(Form("<#DeltaX> %d %s: Raw",imd,ix?"Right":"Left"), 0.1,0.93,kBlack,0.07);
1423 AddPadLabel("Clean", 0.35,0.93,kGreen+2,0.07);
1424 AddPadLabel("Map", 0.42,0.93,kRed,0.07);
1425 AddPadLabel("+Map", 0.5,0.93,kBlue,0.07);
1427 AddPadLabel(Form("#DeltaV:%+.4f | #DeltaT0:%+5.0f | #DeltaX:%+4.0f",
1428 resVDCorr[ix] ? resVDCorr[ix]->GetBinContent(imd-kSDDMin+1):0,
1429 resT0Corr ? resT0Corr->GetBinContent(imd-kSDDMin+1):0,
1430 resXLocCorr ? resXLocCorr->GetBinContent(imd-kSDDMin+1):0),
1431 0.5, 0.15, kRed, 0.07);
1436 for (int ix=0;ix<2;ix++) {
1437 sddCanv->cd(++cntPad);
1438 hsdd = (TH1*)procHistos.At(GetStoreID(kSDVvsZ,imd,ix)); // correction Vd vs Z
1439 PlotHisto(hsdd," ",7,kBlack,1);
1440 hsdd = (TH1*)procHistos.At(GetStoreID(kSDVvsZOrig,imd,ix)); // correction Vd vs Z
1441 PlotHisto(hsdd,"same",24,kBlue,1);
1442 AddPadLabel(Form("<#DeltaV> vs Z %d %s | Stat:%.2e",imd,ix?"Right":"Left",
1443 ((TH1*)procHistos.At(GetStoreID(kSXvsX,imd,ix)))->GetEntries()), 0.1,0.93,kBlack,0.07);
1445 AddPadLabel(Form("<#DeltaV>:%+.4f",resVDCorrZ[ix] ? resVDCorrZ[ix]->GetBinContent(imd-kSDDMin+1):0), 0.5, 0.15, kRed, 0.07);
1453 if (!saved) sddCanv->Print(psnm1.Data());
1454 sddCanv->Print(psnm2.Data());
1457 //__________________________________
1458 TH1* GetPadBaseHisto(TPad* pad)
1460 if (!pad) pad = (TPad*)gPad;
1462 TList* lst = pad->GetListOfPrimitives();
1463 int size = lst->GetSize();
1465 for (int i=0;i<size;i++) {
1466 TObject* obj = lst->At(i);
1468 if (obj->InheritsFrom("TH1")) {hst = (TH1*)obj; break;}
1473 //__________________________________
1474 TLatex* AddPadLabel(const char*txt,float x,float y,int color,float size)
1476 TLatex* lt = new TLatex(x,y,txt);
1478 lt->SetTextColor(color);
1479 lt->SetTextSize(size);
1484 //__________________________________
1485 void SetCDBObjData(int firstrun,int lastrun,const char* comment)
1487 // change range and comment of the objects to store
1488 firstRun = firstrun;
1490 cdbComment = comment;
1493 //__________________________________
1494 void PrepCDBObj(TObject *obj,const char* path,int firstrun,int lastrun,const char* comment)
1496 if (firstrun<0) firstrun = 0;
1498 AliCDBManager* man = AliCDBManager::Instance();
1499 man->UnsetDefaultStorage();
1500 man->SetDefaultStorage("local://");
1501 AliCDBMetaData* md = new AliCDBMetaData();
1502 md->SetResponsible("Ruben Shahoyan");
1503 md->SetComment(comment);
1504 AliCDBId id(path,firstrun,lastrun<=0 ? (AliCDBRunRange::Infinity()) : lastrun);
1505 //AliCDBStorage* st = man->GetStorage("local//.");
1506 man->Put(obj,id,md);
1510 //__________________________________________________________
1511 double edgeLow(double* x, double *par)
1514 // residuals assuming linear dependence of "natural" residual vs Xtrue and smeared
1515 // by the finite track resolution
1517 double sigma = par[1];
1518 double offs = par[2];
1519 double slop = par[3];
1521 if (sigma<1) return 0;
1522 if (x0<-sigma) return 0;
1527 double arg = xex/sigma;
1529 double res = arg<50 ? slop*sigma*TMath::Exp(-arg)/TMath::Sqrt(2*TMath::Pi()) : 0;
1530 double erftrm = 1.+TMath::Erf(xex/sigma/TMath::Sqrt(2));
1531 //printf("%+e %+e %+e\n",x[0],x0,erftrm);
1532 if (xex<0 && TMath::Abs(erftrm)<1e-10) res = -xex*slop;
1533 else res /= erftrm/2.;
1534 res += (offs + (slop-1.)*xex);
1539 //__________________________________________________________
1540 double edgeHigh(double* x, double *par)
1543 // residual assuming linear dependence of "natural" residual vs Xtrue and smeared
1544 // by the finite track resolution
1546 double sigma = par[1];
1547 double offs = par[2];
1548 double slop = par[3];
1549 double tailCorr = par[4];
1551 if (sigma<1) return 0;
1552 if (x0<-sigma) return 0;
1554 double xex = (x0 - x[0])*tailCorr;
1556 double arg = xex/sigma;
1558 double res = arg<50 ? slop*sigma*TMath::Exp(-arg)/TMath::Sqrt(2*TMath::Pi()) : 0;
1559 double erftrm = 1.+TMath::Erf(xex/sigma/TMath::Sqrt(2));
1560 if (xex<0 && TMath::Abs(erftrm)<1e-10) res = xex*slop;
1561 else res /= -erftrm/2.;
1562 res += (offs + (slop-1.)*xex);
1567 //_________________________________________________________________________
1568 void RedoProfileErrors(TH1* profH, TProfile* prof)
1570 // cure errors of low.stat bins
1571 int nbCnt = 0, nbn = prof->GetNbinsX();
1572 double meanStat = 0, meanSpread = 0, wghStat = 0;
1573 for (int i=1;i<=nbn;i++) {
1574 double stat = prof->GetBinEntries(i);
1575 if (stat>0) {meanStat+=stat; nbCnt++;}
1577 if (nbCnt>0) meanStat/= nbCnt; // mean occupancy
1579 for (int i=1;i<=nbn;i++) {
1580 double stat = prof->GetBinEntries(i);
1581 if (stat<meanStat/2) continue;
1582 meanSpread += prof->GetBinError(i)*TMath::Sqrt(stat)*stat;
1585 if (wghStat) meanSpread /= wghStat; // mean spread
1587 for (int i=1;i<=nbn;i++) { // assign error acording to occupancy
1588 double stat = prof->GetBinEntries(i);
1589 if (stat>meanStat/2 || stat<1) continue;
1590 profH->SetBinError(i, meanSpread/TMath::Sqrt(stat));
1594 //_________________________________________________________________________
1595 void CureEdges(TH1* prof)
1597 // cure edges of the profile histo
1598 const double kMaxChi2 = 20.;
1599 const double kSlpDf = 0.05;
1600 static TF1* fitEdgeLow = new TF1("fitEdgeLow" ,edgeLow , -5000, sddSeg.Dx()+5000,4);
1601 static TF1* fitEdgeHigh = new TF1("fitEdgeHigh",edgeHigh, -5000, sddSeg.Dx()+5000,5);
1603 int ndf,ib0,ib1,nbn = prof->GetNbinsX();
1604 double sigma,offs,slp,chi2,x0;
1607 // find 1st non-empty bin
1608 for (ib0=1;ib0<=nbn;ib0++) if (prof->GetBinError(ib0)>1e-9) break;
1609 x0 = prof->GetBinCenter(ib0);
1610 ib1 = prof->FindBin(wDXEdge);
1611 if (ib1-ib0<minBinsEdge) ib1 = ib0+minBinsEdge;
1613 fitEdgeLow->SetParameters(100,100,0,1);
1614 fitEdgeLow->SetParLimits(0,0, sddSeg.Dx());
1615 fitEdgeLow->SetParLimits(1,edgeSmearMinErr, edgeSmearMaxErr);
1616 fitEdgeLow->SetParLimits(3,1.-kSlpDf, 1.+kSlpDf);
1618 prof->Fit(fitEdgeLow,"q","",prof->GetBinLowEdge(ib0)+1, prof->GetBinCenter(ib1+1)-1);
1619 chi2 = fitEdgeLow->GetChisquare();
1620 ndf = fitEdgeLow->GetNDF();
1621 if (ndf>0) chi2 /= ndf;
1623 x0 = fitEdgeLow->GetParameter(0);
1624 sigma = fitEdgeLow->GetParameter(1);
1625 offs = fitEdgeLow->GetParameter(2);
1626 slp = fitEdgeLow->GetParameter(3);
1627 if ( chi2<kMaxChi2) {
1629 ib1 = prof->FindBin(x0);
1630 for (int i=ib0;i<=ib1;i++) {
1631 if (prof->GetBinError(i)<1e-9) continue;
1632 double xb = prof->GetBinCenter(i);
1633 double polval = offs+(slp-1.)*(xb-x0);
1634 if (xb>0) prof->AddBinContent(i, polval - fitEdgeLow->Eval( xb ) );
1635 else prof->SetBinContent(i,polval);
1639 // find last non-empty bin
1640 for (ib1=nbn;ib1>=1;ib1--) if (prof->GetBinError(ib1)>1e-9) break;
1641 x0 = prof->GetBinCenter(ib1);
1642 ib0 = prof->FindBin(sddSeg.Dx() - wDXEdge);
1643 if (ib1-ib0<minBinsEdge) ib0 = ib1-minBinsEdge;
1645 fitEdgeHigh->SetParameters(prof->GetBinCenter(ib0)+wDXEdge-100,100,0,1,1.);
1646 fitEdgeHigh->SetParLimits(0,0, sddSeg.Dx()+150);
1647 fitEdgeHigh->SetParLimits(1,edgeSmearMinErr, edgeSmearMaxErr);
1648 fitEdgeHigh->SetParLimits(3,1.-kSlpDf, 1.+kSlpDf);
1649 fitEdgeHigh->SetParLimits(4,0.3, 3.);
1650 prof->Fit(fitEdgeHigh,"q+","",prof->GetBinLowEdge(ib0)+1, prof->GetBinCenter(ib1+1)+1);
1652 chi2 = fitEdgeHigh->GetChisquare();
1653 ndf = fitEdgeHigh->GetNDF();
1654 if (ndf>0) chi2 /= ndf;
1656 x0 = fitEdgeHigh->GetParameter(0);
1657 sigma = fitEdgeHigh->GetParameter(1);
1658 offs = fitEdgeHigh->GetParameter(2);
1659 slp = fitEdgeHigh->GetParameter(3);
1660 if ( chi2<kMaxChi2 ) {
1662 ib0 = prof->FindBin(x0);
1663 for (int i=ib0;i<=ib1;i++) {
1664 if (prof->GetBinError(i)<1e-9) continue;
1665 double xb = prof->GetBinCenter(i);
1666 double polval = offs+(slp-1.)*(xb-x0);
1667 if (xb<sddSeg.Dx()) prof->AddBinContent(i, polval - fitEdgeHigh->Eval( xb ) );
1668 else prof->SetBinContent(i,polval);
1674 //_________________________________________________________________________
1675 void SafeRebin(TProfile* prof, Int_t factor, Bool_t xprof)
1677 // rebin taking into account left/right margins
1678 const int minBins = 5;
1680 if (factor<1) return;
1681 Bool_t firstX=kTRUE,firstZ=kTRUE;
1683 if (xprof) { // drift profiles
1684 bmn = prof->FindBin(1);
1685 bmx = prof->FindBin(sddSeg.Dx()-1);
1688 double zrange = sddSeg.Dz()/2 - 1.;
1689 bmn = prof->FindBin(-zrange);
1690 bmx = prof->FindBin( zrange);
1692 int nbTot = prof->GetNbinsX();
1693 int nbUse = bmx - bmn + 1;
1694 int edge = bmn-1; // number of edge bins from each side
1696 // find closest divisor
1699 for (int i=2;i<nbUse;i++) {
1700 if ((nbUse%i)==0 && (nbUse/fCClose)>=minBins) {
1701 int dsti = TMath::Abs(factor-i);
1702 if (dsti<dst) {fCClose=i; dst=dsti;}
1706 printf("Could find good rebinning factor\n"); exit(1);
1709 if (fCClose!=factor) {
1710 if ( (xprof&&firstX) || ((!xprof)&&firstZ) ) printf("Rebin%c: For roundness will use factor %d instead of %d\n",xprof ? 'X':'Z',fCClose,factor);
1714 int nbUseNew = nbUse/factor;
1715 if (nbUseNew<minBins) {
1717 if (factor<2) factor=1;
1718 nbUseNew = nbUse/factor;
1721 if ( (xprof&&firstX) || ((!xprof)&&firstZ) ) printf("Rebin%c: Will rebin %d to %d\n",xprof ? 'X':'Z',nbUse,nbUseNew);
1722 if (factor<2) return;
1724 int nbTotNew = nbUseNew + 2*edge;
1725 double *xnew = new double[nbTotNew+1];
1726 TAxis* xax = prof->GetXaxis();
1727 for (int i=1;i<=edge;i++) { // edges are not rebined
1728 xnew[i-1] = xax->GetBinLowEdge(i);
1729 xnew[nbTotNew-i+1] = xax->GetBinLowEdge(nbTot+2-i);
1731 int cnt = 0, bcnt = edge+1;
1732 for (int i=edge+1;i<=nbTot-edge+1;i++) {
1733 if (cnt==0) {xnew[bcnt-1] = xax->GetBinLowEdge(i); bcnt++;}
1734 if (++cnt>=factor) cnt = 0;
1736 TProfile *profNew = (TProfile*)prof->Rebin(nbTotNew,"rbTMPprof$",xnew);
1737 profNew->SetName(prof->GetName());
1738 profNew->SetTitle(prof->GetTitle());
1742 if (xprof&&firstX) firstX = kFALSE;
1743 else if (firstZ) firstZ = kFALSE;
1746 //____________________________________________________
1747 Double_t EdgeFun(double *x, double *par)
1749 // Function to fit Xresiduals vs X, accounting for the limited sensor size and finite gaussian track resolution
1750 // Assumes that the hits density along sensor X has linear dependence rho(x) = a+b*x, and the track have resolution N(x,sig)
1751 // Then, the <residual> seen at coordinate X will be
1752 // R(x) = Integrate[(y-x)*F,{y,t0,t1}]/Integrate[F,{y,t0,t1}];
1754 // F = (a+b*y)*Exp(-(y-x)^2/2/sig^2)
1755 // where t0 and t1 are start and end coordinates of the physical module (in fact, only the coordinate
1756 // of the fitted edge is relevant.
1758 const double sqrt2 = 1.41421356237309515e+00;
1759 const double sqrtPi = 1.77245385090551588e+00;
1763 double sig = par[0]; // resolution
1764 double t0 = par[1]; // active left edge
1765 double t1 = par[2]; // active right edge
1766 double tc0 = t0+par[3]; // constant occupancy left edge
1767 double tc1 = t1-par[4]; // constant occupancy right edge
1768 if (t0>t1) return 1e6;
1769 if (tc0<t0) tc0 = t0;
1770 if (tc1>t1) tc1 = t1;
1771 double ped = par[5];
1773 if (sig<1e-6) return 1e6;
1776 double offset = par[6];
1777 double slope = par[7];
1778 double curve = par[8];
1780 double top=0,norm=0;
1781 for (int it=0;it<3;it++) { // assume three regions of occupance: rise, const, fall
1782 double tmp0,tmp1,a,b;
1787 double rise = tmp1-tmp0;
1788 if (rise<1e-9) continue;
1789 b = (1.-ped)/rise; // linear occupancy rise from ped at t0 to 1 at tc0
1795 if (tmp0>=tmp1) continue;
1796 a = 1.; // constant occupancy between tc0 and tc1
1802 double rise = tmp1-tmp0;
1803 if (rise<1e-9) continue;
1804 b = -(1.-ped)/rise; // linear occupancy fall from 1 at tc1 to ped at t1
1805 a = ped+(1.-ped)*tmp1/rise;
1807 double q0 = (tmp0-px)/sig/sqrt2;
1808 double q1 = (tmp1-px)/sig/sqrt2;
1809 double expq0 = TMath::Abs(q0)<27. ? TMath::Exp(-q0*q0) : 0;
1810 double expq1 = TMath::Abs(q1)<27. ? TMath::Exp(-q1*q1) : 0;
1812 double erfcq0 = TMath::Abs(q0)<27. ? TMath::Erfc(TMath::Abs(q0)) : 0;
1813 double erfcq1 = TMath::Abs(q1)<27. ? TMath::Erfc(TMath::Abs(q1)) : 0;
1816 if (q0>0 && q1>0) derfc = erfcq0 - erfcq1;
1817 else if (q0<0 && q1>0) derfc = (2.-erfcq0) - erfcq1;
1818 else if (q0>0 && q1<0) derfc = erfcq0 - (2.-erfcq1);
1819 else if (q0<0 && q1<0) derfc = erfcq1 - erfcq0;
1821 double topLoc = sig*(expq0*(a+b*tmp0) - expq1*(a+b*tmp1) + b*sqrtPi/sqrt2*sig*derfc);
1822 double nrmLoc = b*(expq0-expq1)*sig + sqrtPi/sqrt2*(a+b*px)*derfc;
1826 printf("it%d | a:%+e b:%+e | topLoc: %+e nrmLoc:%+e -> top: %+e norm: %+e\n",it, a,b,topLoc,nrmLoc,top,norm);
1831 if (TMath::Abs(norm)==0) {
1832 printf("!! x= %f sig=%+e t0=%+e t1=%+e | Top=%e Norm=%e -> %+e\n",px,sig,t0,t1,top,norm,res);
1834 else res = top/norm;
1836 return res + offset + slope*px + curve*px*px;
1840 //____________________________________________________
1841 Double_t EdgeFun(double *x, double *par)
1843 // Function to fit Xresiduals vs X, accounting for the limited sensor size and finite gaussian track resolution
1844 // Assumes that the hits density along sensor X has linear dependence rho(x) = a+b*x, and the track have resolution N(x,sig)
1845 // Then, the <residual> seen at coordinate X will be
1846 // R(x) = Integrate[(y-x)*F,{y,t0,t1}]/Integrate[F,{y,t0,t1}];
1848 // F = (a+b*y)*Exp(-(y-x)^2/2/sig^2)
1849 // where t0 and t1 are start and end coordinates of the physical module (in fact, only the coordinate
1850 // of the fitted edge is relevant.
1852 const double sqrt2 = 1.41421356237309515e+00;
1853 const double sqrtPi = 1.77245385090551588e+00;
1857 double sig = par[0];
1862 if (sig<1e-6) return 0;
1865 double offset = par[5];
1866 double slope = par[6];
1868 double q0 = (t0-px)/sig/sqrt2;
1869 double q1 = (t1-px)/sig/sqrt2;
1870 double expq0 = TMath::Abs(q0)<27. ? TMath::Exp(-q0*q0) : 0;
1871 double expq1 = TMath::Abs(q1)<27. ? TMath::Exp(-q1*q1) : 0;
1873 double erfcq0 = TMath::Abs(q0)<27. ? TMath::Erfc(TMath::Abs(q0)) : 0;
1874 double erfcq1 = TMath::Abs(q1)<27. ? TMath::Erfc(TMath::Abs(q1)) : 0;
1877 if (q0>=0 && q1>=0) derfc = erfcq0 - erfcq1;
1878 else if (q0<=0 && q1>=0) derfc = (2.-erfcq0) - erfcq1;
1879 else if (q0>=0 && q1<=0) derfc = erfcq0 - (2.-erfcq1);
1880 else if (q0<=0 && q1<=0) derfc = erfcq1 - erfcq0;
1882 double top = sig*(expq0*(a+b*t0) - expq1*(a+b*t1) + b*sqrtPi/sqrt2*sig*derfc);
1883 double norm= b*(expq0-expq1)*sig + sqrtPi/sqrt2*(a+b*px)*derfc;
1886 printf("x=%+.2f | q0:%+e q1:%+e exp0:%+e exp1:%+e erf0:%+e erf1:%+e (->%+.e)-> %+e/%+e\n",
1887 px,q0,q1,expq0,expq1,erfcq0,erfcq1, derfc, top,norm);
1891 if (TMath::Abs(norm)<1e-300) {
1892 printf("!! x= %f sig=%+e t0=%+e t1=%+e a=%+e b=%+e | Top=%e Norm=%e -> %+e\n",px,sig,t0,t1,a,b,top,norm,res);
1894 else res = top/norm;
1896 return res + offset + slope*px;
1900 //____________________________________________________________
1901 TH1* GetProfHEntries(TProfile* prof)
1903 // create histo with entries of the profile histo
1904 TH1* hent = ProfileAsTH1(prof, "_Entries");
1905 if (!hent) return 0;
1907 int nb = hent->GetNbinsX()+1;
1908 for (int ib=0;ib<=nb;ib++) hent->SetBinContent(ib, prof->GetBinEntries(ib));
1912 //___________________________________________________________
1913 TH1* FitDXEdges(TProfile* prof)
1916 static TF1 *flft=0,*frgt=0;
1917 const double kMinTot = 1000;
1918 const double kMinThresh = 0.15;
1919 const double kEdgeTol = 2000.;
1920 const double kFitLgt = 10000.;
1921 const double kMinRes = 300.;
1922 const double kMaxRes = 600.;
1923 const double kMaxDX = 5000.;
1924 const double kMaxRiseRange = 10000;
1925 const double kMaxChi2 = 6;
1926 double xspan = sddSeg.Dx();
1929 TH1* histo = ProfileAsTH1(prof,"_h1");
1931 int nb = prof->GetNbinsX();
1932 double tot = prof->GetEntries();
1933 histo->SetBinContent(0,0); // lower active edge
1934 histo->SetBinContent(nb+1,xspan); // upper active edge
1935 if (tot<kMinTot) return histo;
1936 double meanbin = tot/nb;
1938 // RedoProfileErrors(histo,prof);
1939 // find left/right significant bins
1941 for (int i=1;i<=nb;i++) {
1942 double enti = prof->GetBinEntries(i);
1943 if (enti>kMinThresh*meanbin) {bleft=i; break;}
1945 if (bleft<0) return kFALSE;
1946 if (bleft<prof->FindBin(1)) bleft = 1;
1947 double lEdgeMn = TMath::Max(0.,prof->GetBinLowEdge(bleft)-kEdgeTol/2);
1948 double lEdgeMx = lEdgeMn + kEdgeTol;
1951 for (int i=nb;i>0;i--) {
1952 double enti = prof->GetBinEntries(i);
1953 if (enti>kMinThresh*meanbin) {bright=i; break;}
1955 if (bright<0) return kFALSE;
1956 if (bright>prof->FindBin(xspan-1)) bleft = prof->FindBin(xspan-1);
1957 double rEdgeMx = TMath::Min(xspan,prof->GetBinLowEdge(bright)+kEdgeTol/2);
1958 double rEdgeMn = rEdgeMx - kEdgeTol;
1960 printf("Fit left edge\n");
1961 if (!flft) flft = new TF1("lftEdge",EdgeFun,
1962 prof->GetBinLowEdge(1),
1963 prof->GetBinLowEdge(nb+1),
1965 flft->SetParameters((kMaxRes+kMinRes)/2, (lEdgeMn+lEdgeMx)/2,(rEdgeMn+rEdgeMx)/2,lEdgeMx,0,0.5);
1966 flft->SetParLimits(0,kMinRes,kMaxRes);
1967 flft->SetParLimits(1,lEdgeMn,lEdgeMx);
1968 flft->SetParLimits(2,rEdgeMn,rEdgeMx);
1969 flft->SetParLimits(3,0,kMaxRiseRange);
1970 flft->FixParameter(4,0);
1971 flft->SetParLimits(5,0.,1.);
1972 flft->SetParLimits(6,-kMaxDX,kMaxDX);
1973 flft->SetParLimits(7,-kMaxDX/xspan,kMaxDX/xspan);
1974 flft->SetParLimits(8,-kMaxDX/xspan/xspan,kMaxDX/xspan/xspan);
1975 flft->SetLineWidth(1);
1976 double fitLStart = TMath::Max(prof->GetBinLowEdge(1),lEdgeMn-5.*(kMinRes+kMaxRes)/2.);
1977 double fitLEnd = TMath::Min(fitLStart+kFitLgt,rEdgeMn);
1981 prof->Fit(flft,"q+","",fitLStart,fitLEnd);
1982 fstatus = (char*)gMinuit->fCstatu.Data();
1983 if (fstatus.Contains("CONVERGED") || fstatus.Contains("SUCCESSFUL")) {
1985 chiL = flft->GetNDF()>0 ? flft->GetChisquare()/flft->GetNDF() : 0;
1988 TF1* oldf = (TF1*)prof->GetListOfFunctions()->FindObject(flft->GetName());
1989 if (oldf) prof->GetListOfFunctions()->Remove(oldf);
1992 if (chiL<kMaxChi2 && cntL>=100) {
1993 // flft->SetParameter(6,0.);
1994 // flft->SetParameter(7,0.);
1995 // flft->SetParameter(8,0.);
1996 // int mxbin = histo->FindBin(flft->GetParameter(1)+6*flft->GetParameter(0));
1997 double edgL = TMath::Max(skipDXEdge,flft->GetParameter(1)+3*flft->GetParameter(0));
1998 printf("Smoothing left edge up to %.4f\n",edgL);
1999 int mxbin = histo->FindBin(edgL);
2000 for (int i=1;i<=mxbin;i++) {
2001 double x = histo->GetBinCenter(i);
2002 double val = flft->GetParameter(6)+x*(flft->GetParameter(7)+x*flft->GetParameter(8));
2003 histo->SetBinContent(i,val);
2005 histo->SetBinContent(0, TMath::Max(0.0, flft->GetParameter(1)-6*flft->GetParameter(0))); // effective lower sensor edge
2008 printf("Left edge bad: %f %d %s\n",chiL,cntL,fstatus.Data());
2011 printf("Fit right edge\n");
2012 if (!frgt) frgt = new TF1("rgtEdge",EdgeFun,
2013 prof->GetBinLowEdge(1),
2014 prof->GetBinLowEdge(nb+1),
2016 frgt->SetParameters((kMaxRes+kMinRes)/2, (lEdgeMn+lEdgeMx)/2,(rEdgeMn+rEdgeMx)/2,0,(rEdgeMx-rEdgeMn)/2.,0.5);
2017 frgt->SetParLimits(0,kMinRes,kMaxRes);
2018 frgt->SetParLimits(1,lEdgeMn,lEdgeMx);
2019 frgt->SetParLimits(2,rEdgeMn,rEdgeMx);
2020 frgt->FixParameter(3,0);
2021 frgt->SetParLimits(4,0,kMaxRiseRange);
2022 frgt->SetParLimits(5,0.,1.);
2023 frgt->SetParLimits(6,-kMaxDX,kMaxDX);
2024 frgt->SetParLimits(7,-kMaxDX/xspan,kMaxDX/xspan);
2025 frgt->SetParLimits(8,-kMaxDX/xspan/xspan,kMaxDX/xspan/xspan);
2026 frgt->SetLineWidth(1);
2027 double fitREnd = TMath::Min(prof->GetBinLowEdge(nb+1),rEdgeMx+5.*(kMinRes+kMaxRes)/2.);
2028 double fitRStart = TMath::Max(fitREnd-kFitLgt,lEdgeMx);
2032 prof->Fit(frgt,"q+","",fitRStart,fitREnd);
2033 fstatus = (char*)gMinuit->fCstatu.Data();
2034 if (fstatus.Contains("CONVERGED") || fstatus.Contains("SUCCESSFUL")) {
2036 chiR = frgt->GetNDF()>0 ? frgt->GetChisquare()/frgt->GetNDF() : 0;
2039 TF1* oldf = (TF1*)prof->GetListOfFunctions()->FindObject(frgt->GetName());
2040 if (oldf) prof->GetListOfFunctions()->Remove(oldf);
2042 } while((++cntR<3));
2044 if (chiR<kMaxChi2 && cntR>=100) {
2045 // frgt->SetParameter(6,0.);
2046 // frgt->SetParameter(7,0.);
2047 // frgt->SetParameter(8,0.);
2048 // int mnbin = histo->FindBin(frgt->GetParameter(2)-6*frgt->GetParameter(0));
2049 double edgR = TMath::Min(xspan-skipDXEdge,frgt->GetParameter(2)-3*frgt->GetParameter(0));
2050 printf("Smoothing right edge from %.4f\n",edgR);
2051 int mnbin = histo->FindBin(edgR);
2052 for (int i=mnbin;i<=nb;i++) {
2053 double x = histo->GetBinCenter(i);
2054 double val = frgt->GetParameter(6)+x*(frgt->GetParameter(7)+x*frgt->GetParameter(8));
2055 histo->SetBinContent(i,val);
2057 histo->SetBinContent(nb+1, TMath::Min((double)sddSeg.Dx(), frgt->GetParameter(2)+6*frgt->GetParameter(0))); // effective upper sensor edge
2060 printf("Right edge bad: %f %d %s\n",chiR,cntR,fstatus.Data());
2066 double ftPolComb(double* x, double *par)
2068 // fit with combination of 2 polinomials of order par[0]
2069 int ord = int(par[0]);
2071 double brk = par[1];
2076 if (px<=brk) start += npars;
2077 for (int i=npars;i--;) {
2078 // printf("%.1f (%+.1f) | %d %d %e\n",px,brk,start,start+i,par[start+i]);
2079 res = px*res+par[start+i];
2084 //______________________________________________________________
2085 TH1* SimpleMap(TH1* prof)
2087 // get limits as over/under flows
2088 int nb = prof->GetNbinsX();
2089 double lft = prof->GetBinContent(0);
2090 double rgt = prof->GetBinContent(nb+1);
2092 int b0 = prof->FindBin(lft+1);
2093 int b1 = prof->FindBin(rgt-1);
2094 TString nm = prof->GetName(); nm += "_map";
2095 TH1* smap = (TH1*)prof->Clone(nm.Data());
2102 smapf = new TF1("smapf","pol1",prof->GetBinLowEdge(0),prof->GetBinLowEdge(nb+1));
2103 if (TestMapFunction(prof,smapf,lft,rgt)) {delete smapf; break;}
2105 mean = smapf->GetParameter(0);
2110 smapf = new TF1("smapf","pol2",prof->GetBinLowEdge(0),prof->GetBinLowEdge(nb+1));
2111 if (TestMapFunction(prof,smapf,lft,rgt)) {delete smapf; break;}
2115 double middle = (lft+rgt)/2;
2116 int bmid = prof->FindBin(middle);
2117 double a0 = prof->GetBinContent(b0);
2118 double a1 = prof->GetBinContent(bmid+1);
2119 double s0 = prof->GetBinContent(bmid-1);
2120 double s1 = prof->GetBinContent(b1);
2122 s0 = (s0-a0)/((rgt-lft)/2); // slope for 1st part
2123 a0 -= b0*lft; // offset for 1st part
2124 s1 = (s1-a1)/((rgt-lft)/2); // slope for 2nd part
2125 a1 -= s1*middle; // offset for 2nd part
2127 // 3) try pol1 + pol1
2128 smapf = new TF1("smapf",ftPolComb,prof->GetBinLowEdge(0),prof->GetBinLowEdge(nb+1),2+2*2);
2129 smapf->SetParameters(1,middle,a0,b0,a1,b1);
2130 smapf->FixParameter(0,1);
2131 smapf->SetParLimits(1,lft+4*prof->GetBinWidth(nb/2),rgt-4*prof->GetBinWidth(nb/2));
2132 if (TestMapFunction(prof,smapf,lft,rgt)) {delete smapf; break;}
2135 // 3) try pol2 + pol2
2136 smapf = new TF1("smapf",ftPolComb,prof->GetBinLowEdge(0),prof->GetBinLowEdge(nb+1),2+2*3);
2137 smapf->SetParameters(2,middle,a0,b0,0,a1,b1,0);
2138 smapf->FixParameter(0,2);
2139 smapf->SetParLimits(1,lft+4*prof->GetBinWidth(nb/2),rgt-4*prof->GetBinWidth(nb/2));
2140 if (TestMapFunction(prof,smapf,lft,rgt)) {delete smapf; break;}
2146 TF1* fnsel = (TF1*) prof->GetListOfFunctions()->FindObject("smapf");
2147 if (!fnsel) {delete smap; return 0;} // no simple solution
2149 // function is ok, set edges to 0
2150 b0 = smap->FindBin(1);
2151 b1 = smap->FindBin(sddSeg.Dx()-1);
2152 double f0 = fnsel->Eval( smap->GetBinCenter(b0) );
2153 double f1 = fnsel->Eval( smap->GetBinCenter(b1) );
2154 double slp = (f1-f0)/(smap->GetBinCenter(b1) - smap->GetBinCenter(b0));
2156 for (int ib=b0+1;ib<b1;ib++) {
2157 double x = smap->GetBinCenter(ib);
2158 double diff = fnsel->Eval(x) - (f0+slp*x);
2159 smap->SetBinContent(ib, diff);
2165 //______________________________________________________________________
2166 Bool_t TestMapFunction(TH1* smap, TF1* fun, double lft, double rgt)
2168 // test if fun describes the shape
2171 smap->Fit(fun,"0qN","",lft,rgt);
2172 fstatus = (char*)gMinuit->fCstatu.Data();
2173 if ( fstatus.Contains("CONVERGED") || fstatus.Contains("SUCCESSFUL")) {
2174 chi2 = fun->GetNDF()>0 ? fun->GetChisquare()/fun->GetNDF() : 0;
2175 if (chi2<=kMaxChi2SimpleMap) {
2176 fun->SetLineWidth(1);
2177 fun->SetLineStyle(2);
2178 smap->Fit(fun,"q","",lft,rgt);