]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/macrosSDD/CalibrateSDD.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update master to aliroot
[u/mrichter/AliRoot.git] / ITS / macrosSDD / CalibrateSDD.C
CommitLineData
ba2089c4 1#if !defined(__CINT__) || defined(__MAKECINT__)
2#include <TMath.h>
3#include <TH1.h>
4#include <TH2.h>
5#include <TF1.h>
6#include <TProfile.h>
7#include <TSpline.h>
8#include <TCanvas.h>
9#include <TStyle.h>
10#include <TFile.h>
11#include <TGrid.h>
12#include <TGraphErrors.h>
13#include <TSystem.h>
14#include <TLatex.h>
3379b34a 15#include <TMinuit.h>
ba2089c4 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"
25#endif
26
3379b34a 27Bool_t verbose = kFALSE;
ba2089c4 28
29/*
30 // master equations
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
33 //
34
35 where
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
41 //
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)
46 T0_true : real T0
47 deltaT0 : error in T0
48*/
49
50//------------------- Details of CDB objects to create ---------------
51int firstRun = -1;
52int lastRun = -1;
53TString cdbComment = "";
54//--------------------------------------------------------------------
55
56const Int_t kSDDMin=240,kSDDMax=499,kNSDD=kSDDMax-kSDDMin+1;
3379b34a 57enum {kSXvsX,kSXvsZ,kSTDvsX,kSTDvsZ,kSDXvsXclean,kSVDvsX,kSCorrMapX,kSXvsXCorr,kSXvsZCorr,kSDVvsZ,kSDVvsZOrig,kSGrDxx,kSGrTDx,kSGrTXCorr, kNStore};
ba2089c4 58const Int_t maxVZOrd = 4; // max order of polinomial for VD vs Z correction (normal sensors)
59const Double_t chi2BadVZ = 2.9; // treat as bad module if chi2/ndf is worse
60//
61enum {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
3379b34a 64const 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
69};
70const Char_t* kLRNames[] = {"0","1"}; // identifiers for left/right sides in histo names
71/*
ba2089c4 72const 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
77};
78const Char_t* kLRNames[] = {"left","right"}; // identifiers for left/right sides in histo names
3379b34a 79*/
ba2089c4 80//
3379b34a 81const double kMaxChi2SimpleMap = 1.0; // if polinomial fit gives chi2 below this value, make map out of it
ba2089c4 82double minBEntryX = 10; // min entries in the bin to account (vs Xdrift)
83double minBEntryZ = 20; // min entries in the bin to account (vs Z)
84double skipDXEdge = 1000.; // microns to skip from each side of XDrift when smoothing
85double edgeSmearMinErr = 60.; // minimal error for extrapolation
86double edgeSmearMaxErr = 600.; // maximal error for extrapolation
87double wDXEdge = 5000.; // width in microns for XDrift edges to average
88double wVDEdge = 2500; // width in microns for VDrift averaging region on the edges
89double threshClean= 0.3; // clean edge bins with occupancy smaller than in the neighbourhood
90int nVDEdge = 20; // min width in n bins for VDrift averaging region on the edges
91int smoothLevel= 5; // smoothing level
92int minBinsEdge = 5; // min number of X>0 good bins to fit the edge
93Bool_t useSplieForR2V = kFALSE; // if true, extract Vd profile using Derivateve of spline (may be inadequate for sharp turns)
94//
c3793e48 95
96Bool_t userLeftRightFromTASK = kTRUE; // VERY IMPORTANT: histos produced by AliAnalysisTaskITSAlignQA have inverted left/right side definitions
3379b34a 97Bool_t forceT0CorrTo0 = kFALSE;//kTRUE;
ba2089c4 98Bool_t forceRPhiCorrTo0 = kTRUE;
3379b34a 99Bool_t userDummyCorrMap = kFALSE;//kTRUE;
100Bool_t userDummyt0Corr = kFALSE;//kTRUE;
101Bool_t userDummyxlCorr = kFALSE;//kTRUE;
102Int_t userRebinX = 1;
103Int_t userRebinZ = 2;
ba2089c4 104//
105//
106AliITSsegmentationSDD sddSeg;
107//
108//----------------- working variables --------------------------------
109Int_t currSDDId=-1,currSDDSide=-1,currMod=-1;
110TProfile* currHDXvsX = 0; // raw DX vs X histo from QA
111TProfile* currHDXvsZ = 0; // raw DX vs Z histo from QA
112TProfile* currHTDvsX = 0; // raw DriftTime-T0 vs X histo from QA
113TProfile* currHTDvsZ = 0; // raw DriftTime-T0 vs Z histo from QA
114TH1* currHDXvsXclean = 0; // smoothed DX vs X histo
3379b34a 115TH1* currHVDvsX = 0; // extracted VDrift vs X profile
ba2089c4 116TH1* currHCorrMapX = 0; // extracted correction map
117TH1* currHDXvsXCorr = 0; // DX vs X histo after application of the map
118TH1* currHDXvsZCorrMT0 = 0; // DX vs X histo after application of the map and t0 correction
119TH1* currHDVvsZ = 0; // correction for the VDrift vs Z
120TH1* currHDVvsZOrig = 0; // correction for the VDrift vs Z (before error processing)
121TGraphErrors* currGrDxx = 0; // DX final vs XDrift
122TGraphErrors* currGrTDx = 0; // DX final vs TDrift
123TGraphErrors* currGrTXCorr = 0; // TDrift vs XDritf
124TCollection* qaHistos = 0; // input QA collection
125TObjArray procHistos; // processed histos buffer
126AliITSresponseSDD* sddResp=0; // SDD response, updated
127TObjArray *vdarrayOld = 0; // olf VDrifts array
128TObjArray* vDriftArr=0; // VDrifts array, updated
129TObjArray* corrMaps=0; // holder for correction maps
130//
131TH1* resOffsDXraw[2]={0}; // Dx vs X offset at Xdrift=0 raw
132TH1* resOffsDXAMap[2]={0}; // Dx vs X offset (mean) after correction by map
133TH1* resOffsDX[2]={0}; // Dx vs X offset at Xdrift=0 after correction by map (for each side)
134TH1* resVDCorr[2]={0}; // VDrift correction (for each side)
135TH1* resVDMean[2]={0}; // average VDrift (for each side)
136TH1* resVDCorrZ[2]={0}; // mean VDrift vs Z corrections (for each side)
137TH1* resT0Corr = 0; // correction to modules T0
138TH1* resXLocCorr = 0; // correction to module location
139TCanvas* sddCanv = 0; // report canvas
140//
141Bool_t sddRespUpdT0OK = kFALSE; // flag that SDDresponse object T0 was updated
142Bool_t sddRespUpdVDOK = kFALSE; // flag that SDDresponse object VDrift correction was updated
143Bool_t sddVDriftUpdOK = kFALSE; // flag that SDD VDrift object was updated
144//
145TString pathSDDRespOld = "";
146TString pathSDDVDriftOld = "";
147TString pathSDDCorrMapOld = "";
148//
149//--------------------------------------------------------------------
150void Process(const char* pathSDDResp=0,const char* pathSDDVDrift=0, const char* pathSDDCorrMap=0);
151void Process(TCollection* qa, const char* pathSDDResp=0,const char* pathSDDVDrift=0, const char* pathSDDCorrMap=0);
152Bool_t ProcSDDSideDXvsX();
153Bool_t ProcSDDSideDXvsZ();
154TProfile* GetQAHisto(Int_t qaType);
155TProfile* H2Profile(TH2* h2);
156TH1* H2ProfileAsTH1(TH2* h2);
3379b34a 157TH1* GetProfHEntries(TProfile* prof);
ba2089c4 158TH1* ProfileAsTH1(TProfile* prof, const char* addName);
159TH1* CleanProfile(TProfile* profR);
160TH1* Vdrift2Resid(TH1* vd);
161TH1* Resid2Vdrift(TH1* res);
162void RedoProfileErrors(TH1* profH1,TProfile* prof);
163void SetModuleID(Int_t mdID,Int_t side=-1);
164void PrepareModuleHistos(Int_t mdID, Int_t side, Bool_t vsZ);
165void CheckResultDX();
166void CalcDXCorrections();
167double GetVOld(double z);
168Int_t GetStoreID(int type, int imd=-1,int side=-1);
169void CleanPrev();
170void StoreCurrent();
171double ftVdZ(double *x, double *par);
172void PlotReport(const char* psname="repSDDQA.ps");
173TH1* GetPadBaseHisto(TPad* pad);
174Bool_t PlotHisto(TH1* h, Option_t* opt="", int mrkStyle=20,int mrkCol=kBlack, double mrkSize=1.);
175TLatex* AddPadLabel(const char*txt,float x=0.1,float y=0.9,int color=kBlack,float size=0.04);
176void GetHistoMinMaxInRange(TH1* h, double &mn,double &mx);
177double edgeLow(double* x, double *par);
178double edgeHigh(double* x, double *par);
179void CureEdges(TH1* prof);
3379b34a 180void SafeRebin(TProfile* prof, Int_t factor, Bool_t xprof);
181TH1* FitDXEdges(TProfile* prof);
182Bool_t TestMapFunction(TH1* smap, TF1* fun, double lft, double rgt);
183double ftPolComb(double* x, double *par);
184TH1* SimpleMap(TH1* prof);
ba2089c4 185//
186Bool_t LoadSDDVDrift(TString& path, TObjArray *&arr);
187Bool_t LoadSDDResponse(TString& path, AliITSresponseSDD *&resp);
188Bool_t LoadSDDCorrMap(TString& path, TObjArray *&maps);
189AliCDBEntry* GetCDBEntry(const char* path);
190//
191void UpdateSDDResponse(AliITSresponseSDD *resp, Bool_t t0=kTRUE, Bool_t vdrift=kTRUE);
192void UpdateSDDVDrift(AliITSDriftSpeedArraySDD* vdArr, int imd, int side);
193TObjArray* UpdateSDDVDrift();
194TObjArray* CreateCorrMaps();
195AliITSresponseSDD* UpdateSDDResponse(Bool_t t0=kTRUE, Bool_t vdrift=kTRUE);
196AliITSCorrMap1DSDD* CreateCorrMap(TH1* mapHisto, int imd, int side, AliITSCorrMap1DSDD* updateMap=0);
197void PrepCDBObj(TObject *obj,const char* path,int firstrun=0,int lastrun=999999999,const char* comment="");
198
199//-------------------------------------------------------------------
200
201//_________________________________________________________________________
202void Process(TCollection* qa, const char* pathSDDResp, const char* pathSDDVDrift, const char* pathSDDCorrMap)
203{
204 // process all
205 qaHistos = qa;
206 Process(pathSDDResp,pathSDDVDrift,pathSDDCorrMap);
207}
208
209//_________________________________________________________________________
210void Process(const char* pathSDDResp, const char* pathSDDVDrift, const char* pathSDDCorrMap)
211{
212 // process all
213 procHistos.Clear();
214 pathSDDRespOld = pathSDDResp;
215 pathSDDVDriftOld = pathSDDVDrift;
216 pathSDDCorrMapOld = pathSDDCorrMap;
217 sddRespUpdT0OK = sddRespUpdVDOK = sddVDriftUpdOK = kFALSE;
218 //
219 if (pathSDDVDriftOld.IsNull()) printf("Attention: Old VDrift is missing!\n");
220 //
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
227 StoreCurrent();
228 }
229 CalcDXCorrections(); // calculate delta's
230 //
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
235 StoreCurrent();
236 }
237 }
238 //
239 corrMaps = CreateCorrMaps(); // create correction maps
240 PrepCDBObj(corrMaps,"ITS/Calib/MapsTimeSDD",firstRun,lastRun,cdbComment.Data());
241 //
242 if (!pathSDDVDriftOld.IsNull()) {
243 vDriftArr = UpdateSDDVDrift();
244 PrepCDBObj(vDriftArr,"ITS/Calib/DriftSpeedSDD",firstRun,lastRun,cdbComment.Data());
245 }
246 //
247 if (!pathSDDRespOld.IsNull()) {
248 sddResp = UpdateSDDResponse();
249 PrepCDBObj(sddResp,"ITS/Calib/RespSDD",firstRun,lastRun,cdbComment.Data());
250 }
251 //
252 PlotReport();
253}
254
255//_________________________________________________________________________
256Bool_t ProcSDDSideDXvsX()
257{
258 // process one side of the SDD module
259 currHDXvsXCorr = ProfileAsTH1(currHDXvsX,"corrCheck");
260 RedoProfileErrors(currHDXvsXCorr,currHDXvsX);
261 //
262 if ( (currHDXvsXclean=CleanProfile(currHDXvsX)) ) {
263 //
264 delete currHDXvsXCorr;
265 currHDXvsXCorr = (TH1*)currHDXvsXclean->Clone( Form("%s_corrCheck",currHDXvsX->GetName()) );
266 //
3379b34a 267 // try simple solution
268 if (!(currHCorrMapX=SimpleMap(currHDXvsXclean))) {
269 currHVDvsX = Resid2Vdrift(currHDXvsXclean);
270 currHCorrMapX = Vdrift2Resid(currHVDvsX);
271 }
ba2089c4 272 //
273 currHDXvsXCorr->Add(currHCorrMapX,-1);
274 }
275 // check results
276 CheckResultDX();
277 //
278 return kTRUE;
279}
280
281//_________________________________________________________________________________
282Bool_t ProcSDDSideDXvsZ()
283{
284 // extract correction for Vdrift vs Z from DXvsZ profile and mean_drift_time vs Z
285 //
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.;
291 if (!iniDone) {
292 iniDone = kTRUE;
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);
296 }
297 }
298 //
299 int nb = currHDXvsZ->GetNbinsX();
300 if (currHDXvsZ->GetEntries()<minBEntryZ*nb) return kFALSE;
301 //
302 currHDVvsZ = ProfileAsTH1(currHDXvsZ,"_vzcorr");
303 currHDVvsZ->Reset();
304 //
305 currHDXvsZCorrMT0 = ProfileAsTH1(currHDXvsZ,"_zcorrMapT0");
306 currHDXvsZCorrMT0->Reset();
307 //
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;
312 //
313 for (int ib=ib0;ib<=ib1;ib++) {
314 double ne = currHTDvsZ->GetBinEntries(ib);
315 if (ne<1) continue;
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);
325 }
326 currHDXvsZCorrMT0->Fit(fitP0,"q0","");
327 double pedestal = fitP0->GetParameter(0);
328 //
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);
336 if (t<1) continue;
337 // corrections
338 // account for the corrections leading the mean DX to 0
339 dx -= pedestal;
340 double v = dx/t;
341 double ve = TMath::Sqrt(dxe*dxe + v*v*te*te)/t;
342 //
343 nbUse++;
344 vmeanE += 1./(ve*ve);
345 vmean += v/(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;
348 statZB[ibDest] = ne;
349 currHDVvsZ->SetBinContent(ibDest,v);
350 currHDVvsZ->SetBinError(ibDest,ve);
351 //
352 }
353 //
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
357 vmeanE *= nbUse;
358 /*
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);
364 }
365 }
366 */
367 // save original
368 currHDVvsZOrig = (TH1*) currHDVvsZ->Clone(Form("%s_orig",currHDVvsZ->GetName()));
369 // find leftmost good bin
370 int bL,bR;
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);
376 break;
377 }
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);
382 break;
383 }
384 // change bad bins
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));
390 }
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));
396 }
397 //
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;
403 }
404 //
405 // analyse chi2's
406 int bestOrd = 0;
407 double bestChi2 = 9999;
408 for (int iOrd=0;iOrd<=maxVZOrd;iOrd++) {
409 if (chi2s[iOrd]<bestChi2-0.5) {bestChi2 = chi2s[iOrd]; bestOrd = iOrd;}
410 }
411 TF1* fitVvsZ = fitVvsZs[bestOrd];
412 currHDVvsZ->Fit(fitVvsZ,"q");
413 //
414 // extract mean correction neglecting the constant term
415 vmean = 0;
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
420 norm += statZB[ib];
421 }
422 //
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);
427 }
428 //
429 // if (currSDDSide) vmean = -vmean;
430 resVDCorrZ[currSDDSide]->SetBinContent(currMod+1, norm>0 ? vmean/norm : 0);
431 delete statZB;
432 return kTRUE;
433}
434
435//_________________________________________________________________________
436void PrepareModuleHistos(Int_t mdID, Int_t side, Bool_t vsX)
437{
438 // retrieve QA histos, if needed, convert to TH1
439 SetModuleID(mdID, side);
440 if (vsX) {
441 currHDXvsX = GetQAHisto(kDxX);
442 currHTDvsX = GetQAHisto(kTDX);
443 }
444 else {
445 currHTDvsZ = GetQAHisto(kTDZ);
446 currHDXvsZ = GetQAHisto(kDxZ);
447 }
448}
449
450//_________________________________________________________________________
451void SetModuleID(Int_t mdID,Int_t side)
452{
453 // set current module ID,side
454 currSDDId = mdID;
455 currSDDSide = side;
456 currMod = currSDDId - kSDDMin;
457}
458
459//_________________________________________________________________________
460TProfile* GetQAHisto(Int_t qaType)
461{
462 // retrieve needed QA histo
463 if (!qaHistos) {printf("QA Histos collection is not set\n"); exit(1);}
464 //
465 if (currSDDId<kSDDMin || currSDDId>kSDDMax || currSDDSide<0 || currSDDSide>2)
466 {printf("Illegal SDD module ID/Side: %d/%d\n",currSDDId,currSDDSide); exit(1);}
467 //
468 if (qaType<0 || qaType>=kNQATypes)
469 {printf("Illegal QA Histo Type %d\n",qaType); exit(1);}
470 //
c3793e48 471 int trueSide = userLeftRightFromTASK ? (1-currSDDSide) : currSDDSide;
472 const char* hname = Form("%s%d_%s",kQAHistoNames[qaType],currSDDId,kLRNames[trueSide]);
ba2089c4 473 TH1* h = (TH1*) qaHistos->FindObject(hname);
474 if (!h) {printf("Did not find histo %s in QA output collection\n",hname); exit(1);}
475 //
476 if (h->InheritsFrom(TH2::Class())) h = H2Profile((TH2*)h); // make sure it is not TH2 histo
477 //
3379b34a 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);
480 //
ba2089c4 481 return (TProfile*)h;
482}
483
484//_________________________________________________________________________
485TH1* H2ProfileAsTH1(TH2* h2)
486{
487 // extract profile as TH1 histo
488 TProfile* prf = h2->ProfileX();
489 TH1* prof = ProfileAsTH1(prf, "_profH1");
490 delete prf;
491 return prof;
492}
493
494//_________________________________________________________________________
495TH1* ProfileAsTH1(TProfile* prf, const char* addName)
496{
497 // convert profile to TH1 histo
498 TString nm = prf->GetName(); nm += addName;
499 TAxis* xax = prf->GetXaxis();
3379b34a 500 TH1* prof = 0;
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());
504 }
505 else {
506 prof = new TH1F(nm.Data(),nm.Data(),xax->GetNbins(),bins->GetArray());
507 }
ba2089c4 508 for (int i=1;i<=prof->GetNbinsX();i++) {
509 prof->SetBinContent(i, prf->GetBinContent(i));
510 prof->SetBinError(i, prf->GetBinError(i));
511 }
512 return prof;
513}
514
515//_________________________________________________________________________
516TProfile* H2Profile(TH2* h2)
517{
518 // extract profile
519 TString nm = h2->GetName(); nm += "_prof";
520 TProfile* prf = h2->ProfileX(nm.Data());
521 return prf;
522}
523
524//_________________________________________________________________________________
525TH1* CleanProfile(TProfile* profR)
526{
527 // clean profile copy
3379b34a 528 TH1* prof = FitDXEdges(profR); // cure edges
529
ba2089c4 530 int ib0 = profR->FindBin(1), ib1 = profR->FindBin(sddSeg.Dx()-1);
531 int nbn = profR->GetNbinsX();
3379b34a 532 int nSkip = skipDXEdge/profR->GetBinWidth(nbn/2);
533 int nMean = wDXEdge/profR->GetBinWidth(nbn/2);
ba2089c4 534 //
ba2089c4 535 // get mean occupancy skipping first and last nSkip occupied bins
536 //
537 int nbCntL=0, nbCntR=0,nskipL=0,nskipR=0;
538 double smL=0,smR=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;}
543 smL += vl;
544 nbCntL++;
545 }
546 if (nbCntL<1) return 0;
547 smL /= nbCntL;
548 //
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;}
553 smR += vl;
554 nbCntR++;
555 }
556 //
557 if (nbCntR<1) return 0;
558 smR /= nbCntR;
559 //
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;
566 //
567 return prof;
568}
569
570
571//_________________________________________________________________________________
572TH1* Resid2Vdrift(TH1* res)
573{
574 // extract Vdrift profile from residuals profile
575 //
576 TSpline3 spl(res);
577 TString nm = res->GetName(); nm += "_vd";
578 TH1* vd = (TH1*) res->Clone(nm.Data());
579 vd->Reset();
580 int nb = vd->GetNbinsX();
581 //
3379b34a 582 int nbcnt=0, nbfl = wVDEdge/vd->GetBinWidth(nb/2);
ba2089c4 583 if (nbfl<nVDEdge) nbfl=nVDEdge;
584 //
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
591 int nbd = 1;
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;
596 err2 += err*err;
597 nbd++;
598 }
599 if (nbd) err2/=nbd;
600 if (err2>0 && v>0 && v<3) {vAv += v/err2; eAv += 1./err2;}
3379b34a 601 err2 = TMath::Sqrt(err2)/res->GetBinWidth(nb/2);
ba2089c4 602 vd->SetBinError(i, err2 );
603 nbcnt++;
604 }
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);
609 //
610 int ib = 0;
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);
618 }
619 //
620 double vL=0,eL=0,vR=0,eR=0;
621 // get the mean of leftmost stable vdrift
622 int lastCounted = 0;
623 nbcnt = 0;
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;
628 eL += 1./berr/berr;
629 nbcnt++;
630 lastCounted = i;
631 }
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 !!
638 //
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);
646 }
647 nbcnt= 0;
648 lastCounted = 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;
653 eR += 1./berr/berr;
654 nbcnt++;
655 lastCounted = i;
656 }
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)));
666 //
667 return vd;
668}
669
670
671//_________________________________________________________________________________
672TH1* Vdrift2Resid(TH1* vd)
673{
674 // convert Vdrift profile to residuals profile (final correction map)
675 //
676 TString nm = vd->GetName(); nm += "_vd2res";
677 TH1* res = (TH1*) vd->Clone(nm.Data());
678 res->Reset();
3379b34a 679 if (userDummyCorrMap) return res;
ba2089c4 680 int ib0 = res->FindBin(1);
681 int ib1 = res->FindBin(sddSeg.Dx()-1);
682 double resv=0;
683 int lastB=0;
684 for (lastB=ib1+1;lastB--;) if (vd->GetBinError(lastB)>1e-9) break; // find last good bin
685 double lastX = vd->GetBinCenter(lastB);
686 // extend by 1mm
687 lastX += 1000;
688 if (lastX>sddSeg.Dx()) lastX = sddSeg.Dx();
689 lastB = vd->FindBin(lastX);
690 //
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);
696 }
697 double vcorr = (resv)/lastX;
698 //
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);
706 }
707 //
708 return res;
709}
710
711//__________________________________________________________________________________
712void CheckResultDX()
713{
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);
718 //
719 if (currHDXvsX->GetEntries()<minBEntryX*currHDXvsX->GetNbinsX()) {
720 if (currHCorrMapX) currHCorrMapX->Reset();
721 return;
722 }
723 //
724 // vdrift correction
725 int b0 = currHDXvsXCorr->FindBin(1),b1 = currHDXvsXCorr->FindBin(sddSeg.Dx()-1);
726 int nb = 0;
727 for (int ib=b0;ib<b1;ib++) if (currHTDvsX->GetBinEntries(ib)>=minBEntryX) nb++;
3379b34a 728 currHDXvsX->Fit(fitP0,"q0N","");
ba2089c4 729 double offsRaw = fitP0->GetParameter(0);
3379b34a 730 currHDXvsXCorr->Fit(fitP0,"q0N","");
ba2089c4 731 double offsAMap = fitP0->GetParameter(0);
732
733 //
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;
738 int ip = 0;
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));
749 //
750 currGrTDx->SetPoint(ip,t, currHDXvsXCorr->GetBinContent(i));
751 currGrTDx->SetPointError(ip, currHTDvsX->GetBinError(i), currHDXvsXCorr->GetBinError(i));
752 //
753 currGrTXCorr->SetPoint(ip,t, currHTDvsX->GetBinCenter(i));
754 currGrTXCorr->SetPointError(ip,currHTDvsX->GetBinError(i), currHTDvsX->GetBinWidth(i));
755 //
756 ip++;
757 }
758 double del = tmax-tmin;
759 tmin += kFOffs*del;
760 tmax -= kFOffs*del;
761 del = xmax - xmin;
762 xmin += kFOffs*del;
763 xmax -= kFOffs*del;
764 //
c3793e48 765 fitP1->SetParameters(0,0);
ba2089c4 766 currGrDxx->Fit(fitP1,"q","",xmin, xmax);
767 double offs = fitP1->GetParameter(0); // offset of correction line at Xdrift=0
768 //
769 // printf("Fitting VD correction in the range %.1f:%.1f\n",tmin,tmax);
c3793e48 770 fitP1->SetParameters(0,0);
ba2089c4 771 currGrTDx->Fit(fitP1,"q","",tmin,tmax);
772 double vcor = fitP1->GetParameter(1);
773 //
c3793e48 774 fitP1->SetParameters(0,0);
ba2089c4 775 currGrTXCorr->Fit(fitP1,"q","",tmin,tmax);
776 double vav = fitP1->GetParameter(1);
777 //
778 // store results
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);
783 }
784 //
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);
789 }
790 //
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);
795 }
796 //
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);
801 }
802 //
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);
807 }
808 //
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);
814 //
815}
816
817//__________________________________________________________________________
818void CalcDXCorrections()
819{
820 // estimate time0 and alignment correction for the whole module
821 if (!resT0Corr) {
822 resT0Corr = new TH1F("T0Corr","T0 Correction",kNSDD,kSDDMin-0.5,kSDDMax+0.5);
823 resT0Corr->SetMarkerColor(2);
824 resT0Corr->SetMarkerStyle(20);
825 }
826 //
827 if (!resXLocCorr) {
828 resXLocCorr = new TH1F("XLocCorr","XLoc Correction",kNSDD,kSDDMin-0.5,kSDDMax+0.5);
829 resXLocCorr->SetMarkerColor(2);
830 resXLocCorr->SetMarkerStyle(20);
831 }
832 //
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);
839 //
840 double vsum=0,t0Corr=0,xlCorr=0;
841 if (vL>1 && vR>1) { // both sides available
842 vsum = vL + vR;
843 t0Corr = -(offsL+offsR)/vsum;
844 xlCorr = -(offsL*vR - offsR*vL)/vsum;
845 }
846 /*
847 else if (vL>1) t0Corr = -offsL/vL; // only one side is available
848 else if (vR>1) t0Corr = -offsR/vR;
849 */
850 else if (vL>1) xlCorr = -offsL; // only one side is available
851 else if (vR>1) xlCorr = offsR;
852 //
3379b34a 853 if (userDummyt0Corr) t0Corr = 0;
854 if (userDummyxlCorr) xlCorr = 0;
ba2089c4 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
858 //
859 double addMap[2]={0,0};
860 Bool_t redoMaps = kFALSE;
861 //
862 if (forceT0CorrTo0) { // T0 correction was forced to be 0, attribute this to map
863 addMap[0] -= vL*t0Corr;
864 addMap[1] -= vR*t0Corr;
865 redoMaps = kTRUE;
866 }
867 if (forceRPhiCorrTo0) { // alignment correction was forced to be 0, attribute this to map
868 addMap[0] -= xlCorr;
869 addMap[1] -= -xlCorr;
870 redoMaps = kTRUE;
871 }
872 //
873 if (redoMaps) {
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]);
883 }
884 }
885 }
886
887 //
888}
889
890//______________________________________________________________
891Int_t GetStoreID(int type, int imd,int side)
892{
893 // entry of the histo/graph of type in the procHistos array
894 //
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);
899 exit(1);
900 }
901 return (2*(imd-kSDDMin)+side)*kNStore + type;
902}
903
904//______________________________________________________________
905void CleanPrev()
906{
907 // clean "current" objects from last event
908 currHDXvsX = 0;
909 currHDXvsZ = 0;
910 currHTDvsX = 0;
911 currHTDvsZ = 0;
912 currHDXvsXclean = 0;
3379b34a 913 currHVDvsX = 0;
ba2089c4 914 currHCorrMapX = 0;
915 currHDXvsXCorr = 0;
916 currHDVvsZ = 0;
917 currGrDxx = 0;
918 currGrTDx = 0;
919 currGrTXCorr = 0;
920 //
921}
922
923//______________________________________________________________
924void StoreCurrent()
925{
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));
3379b34a 932 if (currHVDvsX) procHistos.AddAtAndExpand(currHVDvsX, GetStoreID(kSVDvsX));
ba2089c4 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));
941 //
942}
943
944//_________________________________________________________________________________
945TObjArray* CreateCorrMaps()
946{
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());
953 exit(1);
954 }
955 //
956 dest->Clear();
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);
965 }
966 }
967 //
968 return dest;
969}
970
971//_________________________________________________________________________________
972AliITSCorrMap1DSDD* CreateCorrMap(TH1* mapHisto, int imd, int side, AliITSCorrMap1DSDD* updateMap)
973{
974 // create or update correction map from histo
975 int nbCorr = 1, nbOld = 0;
976 int b0=0,b1=0;
977 if (mapHisto) {
978 b0 = mapHisto->FindBin(1);
979 b1 = mapHisto->FindBin(sddSeg.Dx()-1);
980 }
981 nbCorr = b1-b0+1;
982 AliITSCorrMap1DSDD* mpCorr = 0;
983 //
977bf0d4 984 // check if the updateMap is meaningful
985 if (updateMap && updateMap->GetNBinsDrift()>2 && nbCorr>1) {
ba2089c4 986 if (mapHisto) {
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));
993 }
994 }
995 mpCorr = updateMap;
996 }
997 else {
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));
1001 }
1002 //
1003 return mpCorr;
1004}
1005
1006//_________________________________________________________________________________
1007TObjArray* UpdateSDDVDrift()
1008{
1009 // retrieve SDD VDrift object and update it
1010 if (!vdarrayOld && !LoadSDDVDrift(pathSDDVDriftOld,vdarrayOld)) return 0;
1011 TObjArray *vdarrayUp = new TObjArray(2*kNSDD);
1012 //
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);
1023 }
1024 }
1025 //
1026 sddVDriftUpdOK = kTRUE;
1027 return vdarrayUp;
1028}
1029
1030
1031//_________________________________________________________________________________
1032void UpdateSDDVDrift(AliITSDriftSpeedArraySDD* vdArr, int imd, int side)
1033{
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) );
1038 if (!vdh) {
1039 //printf("VDrift vs Z correction is not processed for module %d/%d\n",imd,side);
1040 return;
1041 }
1042 TF1* fp = vdh->GetFunction("fitVvsZ");
1043 if (!fp) {printf("VDrift vs Z correction fit is missing SDD%d/%d\n",imd,side); return;}
1044 //
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);
1051 }
1052 //
1053}
1054
1055//_________________________________________________________________________________
1056AliITSresponseSDD* UpdateSDDResponse(Bool_t t0, Bool_t vdrift)
1057{
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;
1064 //
1065 return resp;
1066}
1067
1068//_________________________________________________________________________________
1069void UpdateSDDResponse(AliITSresponseSDD *resp, Bool_t t0, Bool_t vdrift)
1070{
1071 // update the map with extracted values
1072 printf("Updating RespSDD object: T0:%s VDrift:%s\n",t0?"ON":"OFF",vdrift?"ON":"OFF");
1073 //
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);}
1078 //
1079 for (int imd=kSDDMin;imd<=kSDDMax;imd++) {
1080 if (t0 && !forceT0CorrTo0) resp->SetModuleTimeZero(imd, resp->GetTimeZero(imd) - resT0Corr->GetBinContent(imd-kSDDMin+1));
1081 if (vdrift) {
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);
1086 }
1087 }
1088 }
1089 //
1090}
1091
1092//___________________________________________________________________
1093double GetVOld(double z)
1094{
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);
1100 return v;
1101}
1102
1103//___________________________________________________________________
1104double ftVdZ(double *x, double *par)
1105{
1106 // function to fit the vdrift dependence on Z
1107 //
1108 // convert Z to anode
1109 double z = x[0];
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();
1114 //
1115 int ord = int(par[0]);
1116 double v = par[ord+1];
1117 for (int i=ord;i--;) v = par[i+1]+ian*v;
1118 return v;
1119}
1120
1121//________________________________________________________________________________________________________
1122Bool_t LoadSDDVDrift(TString& path, TObjArray *&arr)
1123{
1124 // load VDrift object
1125 if (path.IsNull()) return kFALSE;
1126 printf("Loading SDD VDrift from %s\n",path.Data());
1127 //
1128 AliCDBEntry *entry = 0;
1129 delete arr;
1130 arr = 0;
1131 while(1) {
1132 if (path.BeginsWith("path: ")) { // must load from OCDB
1133 entry = GetCDBEntry(path.Data());
1134 if (!entry) break;
1135 arr = (TObjArray*) entry->GetObject();
1136 entry->SetObject(NULL);
1137 entry->SetOwner(kTRUE);
1138 break;
1139 }
1140 //
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);
1147 else arr = 0;
1148 entry->SetObject(NULL);
1149 entry->SetOwner(kTRUE);
1150 delete entry;
1151 }
1152 //
1153 precf->Close();
1154 delete precf;
1155 break;
1156 }
1157 //
1158 if (!arr) {printf("Failed to load SDD vdrift from %s\n",path.Data()); return kFALSE;}
1159 arr->SetOwner(kTRUE);
1160 return kTRUE;
1161}
1162
1163//________________________________________________________________________________________________________
1164Bool_t LoadSDDResponse(TString& path, AliITSresponseSDD *&resp)
1165{
1166 // load SDD response
1167 if (path.IsNull()) return kFALSE;
1168 printf("Loading SDD response from %s\n",path.Data());
1169 //
1170 AliCDBEntry *entry = 0;
1171 delete resp;
1172 resp = 0;
1173 //
1174 while(1) {
1175 if (path.BeginsWith("path: ")) { // must load from OCDB
1176 entry = GetCDBEntry(path.Data());
1177 if (!entry) break;
1178 resp = (AliITSresponseSDD*) entry->GetObject();
1179 entry->SetObject(NULL);
1180 entry->SetOwner(kTRUE);
1181 break;
1182 }
1183 //
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);
1190 else resp = 0;
1191 entry->SetObject(NULL);
1192 entry->SetOwner(kTRUE);
1193 delete entry;
1194 }
1195 //
1196 precf->Close();
1197 delete precf;
1198 break;
1199 }
1200 //
1201 if (!resp) {printf("Error: Failed to load SDD response from %s\n",path.Data()); return kFALSE;}
1202 return kTRUE;
1203}
1204
1205//________________________________________________________________________________________________________
1206Bool_t LoadSDDCorrMap(TString& path, TObjArray *&maps)
1207{
1208 // Load SDD correction map
1209 //
1210 if (path.IsNull()) return kFALSE;
1211 printf("Loading SDD Correction Maps from %s\n",path.Data());
1212 //
1213 AliCDBEntry *entry = 0;
1214 delete maps;
1215 maps = 0;
1216 while(1) {
1217 if (path.BeginsWith("path: ")) { // must load from OCDB
1218 entry = GetCDBEntry(path.Data());
1219 if (!entry) break;
1220 maps = (TObjArray*) entry->GetObject();
1221 entry->SetObject(NULL);
1222 entry->SetOwner(kTRUE);
1223 break;
1224 }
1225 //
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);
1232 else maps = 0;
1233 entry->SetObject(NULL);
1234 entry->SetOwner(kTRUE);
1235 delete entry;
1236 }
1237 //
1238 precf->Close();
1239 delete precf;
1240 break;
1241 }
1242 //
1243 if (!maps) {printf("Failed to load SDD Correction Map from %s\n",path.Data()); return kFALSE;}
1244
1245 return kTRUE;
1246}
1247
1248//_______________________________________________________________________________________
1249AliCDBEntry* GetCDBEntry(const char* path)
1250{
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);
1256 if (!cdbId) {
1257 printf("Failed to create cdbId\n");
1258 return 0;
1259 }
1260 //
1261 AliCDBStorage* stor = man->GetDefaultStorage();
1262 if (!stor && !man->GetRaw()) man->SetDefaultStorage("raw://");
1263 if (man->GetRaw()) man->SetRun(cdbId->GetFirstRun());
1264 if (stor) {
1265 TString tp = stor->GetType();
1266 if (tp.Contains("alien",TString::kIgnoreCase) && !gGrid) TGrid::Connect("alien:");
1267 }
1268 entry = man->Get(cdbId->GetPath(),cdbId->GetFirstRun(),cdbId->GetVersion(),cdbId->GetSubVersion());
1269 // entry = man->Get( *cdbId );
1270 man->ClearCache();
1271 //
1272 delete cdbId;
1273 return entry;
1274 //
1275}
1276//
1277
1278//_______________________________________________________________________________________
1279Bool_t PlotHisto(TH1* h, Option_t* opt, int mrkStyle,int mrkCol, double mrkSize)
1280{
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);
1288 }
1289 h->SetLineColor(mrkCol);
1290 h->Draw(opt);
1291 //
1292 h->SetMinimum(); h->SetMaximum();
1293 double hmn=h->GetMinimum(),hmx=h->GetMaximum(); // new histo min/max
1294 //
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);
1302 }
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 );
1308 gPad->Update();
1309 return kTRUE;
1310}
1311
1312//_______________________________________________________________________________________
1313void GetHistoMinMaxInRange(TH1* h, double &mn,double &mx)
1314{
1315 // compute min/max of histo in the user range
1316 mn = 1e50;
1317 mx =-1e50;
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;
1325 }
1326}
1327
1328//_______________________________________________________________________________________
1329void PlotReport(const char* psname)
1330{
1331 // report results
1332 sddCanv = new TCanvas("sddCanv","sddCanv",700,900);
1333 //
1334 gStyle->SetOptStat(0);
1335 gStyle->SetOptTitle(0);
1336 //
1337 TString psnm1 = psname;
1338 if (psnm1.IsNull()) psnm1 = "sddQAreport.ps";
1339 TString psnm0 = psnm1 + "[";
1340 TString psnm2 = psnm1 + "]";
1341 sddCanv->Print(psnm0.Data());
1342 //
1343 // mean corrections per module/side
1344 sddCanv->Clear();
1345 sddCanv->Divide(2,3);
1346 int cntPad = 0;
1347 //
1348 for (int ix=0;ix<2;ix++) { // mean residuals before/after correction
1349 sddCanv->cd(++cntPad);
3379b34a 1350 PlotHisto(resOffsDXraw[ix],"p" ,7,kBlack,0.5);
ba2089c4 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);
1354 }
1355 //
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);
1362 }
1363 //
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);
1368 //
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);
1373 //
1374 sddCanv->cd();
1375 sddCanv->Print(psnm1.Data());
1376 //
1377 //-------------------------------------------------------------------
1378 TH1* hsdd = 0;
1379 //
1380 cntPad = 999;
1381 int nModPerPage = 3;
1382 int nRowPerMod = 2;
1383 Bool_t saved = kFALSE;
1384 int ib0=1,ib1=999;
1385 //
1386 for (int imd=kSDDMin;imd<=kSDDMax;imd++) {
1387 if (cntPad>=2*nModPerPage*nRowPerMod) {
1388 sddCanv->cd();
1389 if (imd!=kSDDMin) sddCanv->Print(psnm1.Data());
1390 sddCanv->Clear();
1391 sddCanv->Divide(2,nModPerPage*nRowPerMod);
1392 cntPad = 0;
1393 saved = kTRUE;
1394 }
1395 for (int ix=0;ix<2;ix++) {
1396 sddCanv->cd(++cntPad);
3379b34a 1397 TH1* hsddcl = (TH1*)procHistos.At(GetStoreID(kSDXvsXclean,imd,ix)); // raw residuals
1398 if (hsddcl) {
1399 ib0 = hsddcl->FindBin(1);
1400 ib1 = hsddcl->FindBin(sddSeg.Dx()-1);
1401 hsddcl->GetXaxis()->SetRange(ib0,ib1);
1402 }
1403 PlotHisto(hsddcl,"p",24,kGreen+2,0.5);
1404 //
ba2089c4 1405 hsdd = (TH1*)procHistos.At(GetStoreID(kSXvsX,imd,ix)); // raw residuals
1406 if (hsdd) {
1407 ib0 = hsdd->FindBin(1);
1408 ib1 = hsdd->FindBin(sddSeg.Dx()-1);
1409 hsdd->GetXaxis()->SetRange(ib0,ib1);
1410 }
3379b34a 1411 PlotHisto(hsdd,"p same",20,kBlack,0.4);
1412 PlotHisto(hsddcl,"p sames",24,kGreen+2,0.5);
1413 //
1414 //
ba2089c4 1415 hsdd = (TH1*)procHistos.At(GetStoreID(kSCorrMapX,imd,ix)); // map
1416 if (hsdd) hsdd->GetXaxis()->SetRange(ib0,ib1);
3379b34a 1417 PlotHisto(hsdd,"same",7,kRed,0.5);
1418 hsdd = (TH1*)procHistos.At(GetStoreID(kSXvsXCorr,imd,ix));
ba2089c4 1419 if (hsdd) hsdd->GetXaxis()->SetRange(ib0,ib1);
3379b34a 1420 PlotHisto(hsdd,"histo same",7,kBlue,0.5);
ba2089c4 1421 //
1422 AddPadLabel(Form("<#DeltaX> %d %s: Raw",imd,ix?"Right":"Left"), 0.1,0.93,kBlack,0.07);
3379b34a 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);
ba2089c4 1426 //
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);
1432 //
1433 saved = kFALSE;
1434 }
1435 //
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);
1444 //
1445 AddPadLabel(Form("<#DeltaV>:%+.4f",resVDCorrZ[ix] ? resVDCorrZ[ix]->GetBinContent(imd-kSDDMin+1):0), 0.5, 0.15, kRed, 0.07);
1446 //
1447 saved = kFALSE;
1448 }
1449 //
1450 }
1451 //
1452 sddCanv->cd();
1453 if (!saved) sddCanv->Print(psnm1.Data());
1454 sddCanv->Print(psnm2.Data());
1455}
1456
1457//__________________________________
1458TH1* GetPadBaseHisto(TPad* pad)
1459{
1460 if (!pad) pad = (TPad*)gPad;
1461 if (!pad) return 0;
1462 TList* lst = pad->GetListOfPrimitives();
1463 int size = lst->GetSize();
1464 TH1* hst=0;
1465 for (int i=0;i<size;i++) {
1466 TObject* obj = lst->At(i);
1467 if (!obj) continue;
1468 if (obj->InheritsFrom("TH1")) {hst = (TH1*)obj; break;}
1469 }
1470 return hst;
1471}
1472
1473//__________________________________
1474TLatex* AddPadLabel(const char*txt,float x,float y,int color,float size)
1475{
1476 TLatex* lt = new TLatex(x,y,txt);
1477 lt->SetNDC();
1478 lt->SetTextColor(color);
1479 lt->SetTextSize(size);
1480 lt->Draw();
1481 return lt;
1482}
1483
1484//__________________________________
1485void SetCDBObjData(int firstrun,int lastrun,const char* comment)
1486{
1487 // change range and comment of the objects to store
1488 firstRun = firstrun;
1489 lastRun = lastrun;
1490 cdbComment = comment;
1491}
1492
1493//__________________________________
1494void PrepCDBObj(TObject *obj,const char* path,int firstrun,int lastrun,const char* comment)
1495{
1496 if (firstrun<0) firstrun = 0;
1497 //
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);
1507 //
1508}
1509
1510//__________________________________________________________
1511double edgeLow(double* x, double *par)
1512{
1513 // Low TDrift edge:
1514 // residuals assuming linear dependence of "natural" residual vs Xtrue and smeared
1515 // by the finite track resolution
1516 double x0 = par[0];
1517 double sigma = par[1];
1518 double offs = par[2];
1519 double slop = par[3];
1520 //
1521 if (sigma<1) return 0;
1522 if (x0<-sigma) return 0;
1523 //
1524 double xex = x[0];
1525 xex -= x0;
1526 //
1527 double arg = xex/sigma;
1528 arg *= arg/2;
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);
1535 return res;
1536 //
1537}
1538
1539//__________________________________________________________
1540double edgeHigh(double* x, double *par)
1541{
1542 // High TDrift edge
1543 // residual assuming linear dependence of "natural" residual vs Xtrue and smeared
1544 // by the finite track resolution
1545 double x0 = par[0];
1546 double sigma = par[1];
1547 double offs = par[2];
1548 double slop = par[3];
1549 double tailCorr = par[4];
1550 //
1551 if (sigma<1) return 0;
1552 if (x0<-sigma) return 0;
1553 //
1554 double xex = (x0 - x[0])*tailCorr;
1555 //
1556 double arg = xex/sigma;
1557 arg *= arg/2;
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);
1563 return res;
1564 //
1565}
1566
1567//_________________________________________________________________________
1568void RedoProfileErrors(TH1* profH, TProfile* prof)
1569{
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++;}
1576 }
1577 if (nbCnt>0) meanStat/= nbCnt; // mean occupancy
1578 //
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;
1583 wghStat += stat;
1584 }
1585 if (wghStat) meanSpread /= wghStat; // mean spread
1586 //
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));
1591 }
1592}
1593
1594//_________________________________________________________________________
1595void CureEdges(TH1* prof)
1596{
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);
1602 //
1603 int ndf,ib0,ib1,nbn = prof->GetNbinsX();
1604 double sigma,offs,slp,chi2,x0;
1605 //
1606 // LowT edge
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;
1612 //
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);
1617 //
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;
1622 //
1623 x0 = fitEdgeLow->GetParameter(0);
1624 sigma = fitEdgeLow->GetParameter(1);
1625 offs = fitEdgeLow->GetParameter(2);
1626 slp = fitEdgeLow->GetParameter(3);
1627 if ( chi2<kMaxChi2) {
1628 x0 += 3*sigma;
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);
1636 }
1637 }
1638 //
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;
1644 //
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);
1651 //
1652 chi2 = fitEdgeHigh->GetChisquare();
1653 ndf = fitEdgeHigh->GetNDF();
1654 if (ndf>0) chi2 /= ndf;
1655 //
1656 x0 = fitEdgeHigh->GetParameter(0);
1657 sigma = fitEdgeHigh->GetParameter(1);
1658 offs = fitEdgeHigh->GetParameter(2);
1659 slp = fitEdgeHigh->GetParameter(3);
1660 if ( chi2<kMaxChi2 ) {
1661 x0 -= 3*sigma;
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);
1669 }
1670 }
1671 //
1672}
3379b34a 1673
1674//_________________________________________________________________________
1675void SafeRebin(TProfile* prof, Int_t factor, Bool_t xprof)
1676{
1677 // rebin taking into account left/right margins
1678 const int minBins = 5;
1679 int bmn,bmx;
1680 if (factor<1) return;
1681 Bool_t firstX=kTRUE,firstZ=kTRUE;
1682 //
1683 if (xprof) { // drift profiles
1684 bmn = prof->FindBin(1);
1685 bmx = prof->FindBin(sddSeg.Dx()-1);
1686 }
1687 else { // Z profile
1688 double zrange = sddSeg.Dz()/2 - 1.;
1689 bmn = prof->FindBin(-zrange);
1690 bmx = prof->FindBin( zrange);
1691 }
1692 int nbTot = prof->GetNbinsX();
1693 int nbUse = bmx - bmn + 1;
1694 int edge = bmn-1; // number of edge bins from each side
1695 //
1696 // find closest divisor
1697 int fCClose = 2;
1698 int dst = 9999;
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;}
1703 }
1704 }
1705 if (dst==9999) {
1706 printf("Could find good rebinning factor\n"); exit(1);
1707 }
1708 //
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);
1711 factor = fCClose;
1712 }
1713 //
1714 int nbUseNew = nbUse/factor;
1715 if (nbUseNew<minBins) {
1716 factor = nbUse/5;
1717 if (factor<2) factor=1;
1718 nbUseNew = nbUse/factor;
1719 }
1720 //
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;
1723 //
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);
1730 }
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;
1735 }
1736 TProfile *profNew = (TProfile*)prof->Rebin(nbTotNew,"rbTMPprof$",xnew);
1737 profNew->SetName(prof->GetName());
1738 profNew->SetTitle(prof->GetTitle());
1739 *prof = *profNew;
1740 delete profNew;
1741 //
1742 if (xprof&&firstX) firstX = kFALSE;
1743 else if (firstZ) firstZ = kFALSE;
1744}
1745
1746//____________________________________________________
1747Double_t EdgeFun(double *x, double *par)
1748{
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}];
1753 // with
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.
1757 //
1758 const double sqrt2 = 1.41421356237309515e+00;
1759 const double sqrtPi = 1.77245385090551588e+00;
1760 double px = x[0];
1761 //
1762 // edge parameters
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];
1772 //
1773 if (sig<1e-6) return 1e6;
1774 //
1775 // slope parameters
1776 double offset = par[6];
1777 double slope = par[7];
1778 double curve = par[8];
1779 //
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;
1783 //
1784 if (it==0) {
1785 tmp0 = t0;
1786 tmp1 = tc0;
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
1790 a = ped-t0*b;
1791 }
1792 else if (it==1) {
1793 tmp0 = tc0;
1794 tmp1 = tc1;
1795 if (tmp0>=tmp1) continue;
1796 a = 1.; // constant occupancy between tc0 and tc1
1797 b = 0;
1798 }
1799 else {
1800 tmp0 = tc1;
1801 tmp1 = t1;
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;
1806 }
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;
1811 //
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;
1814 //
1815 double derfc = 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;
1820 //
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;
1823 top += topLoc;
1824 norm+= nrmLoc;
1825 if (verbose) {
1826 printf("it%d | a:%+e b:%+e | topLoc: %+e nrmLoc:%+e -> top: %+e norm: %+e\n",it, a,b,topLoc,nrmLoc,top,norm);
1827 }
1828 }
1829 //
1830 double res = 0;
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);
1833 }
1834 else res = top/norm;
1835 //
1836 return res + offset + slope*px + curve*px*px;
1837}
1838
1839/*
1840//____________________________________________________
1841Double_t EdgeFun(double *x, double *par)
1842{
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}];
1847 // with
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.
1851 //
1852 const double sqrt2 = 1.41421356237309515e+00;
1853 const double sqrtPi = 1.77245385090551588e+00;
1854 double px = x[0];
1855 //
1856 // edge parameters
1857 double sig = par[0];
1858 double t0 = par[1];
1859 double t1 = par[2];
1860 double a = par[3];
1861 double b = par[4];
1862 if (sig<1e-6) return 0;
1863 //
1864 // slope parameters
1865 double offset = par[5];
1866 double slope = par[6];
1867 //
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;
1872 //
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;
1875 //
1876 double derfc = 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;
1881 //
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;
1884 //
1885 if (verbose) {
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);
1888 }
1889 //
1890 double res = 0;
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);
1893 }
1894 else res = top/norm;
1895 //
1896 return res + offset + slope*px;
1897}
1898*/
1899
1900//____________________________________________________________
1901TH1* GetProfHEntries(TProfile* prof)
1902{
1903 // create histo with entries of the profile histo
1904 TH1* hent = ProfileAsTH1(prof, "_Entries");
1905 if (!hent) return 0;
1906 hent->Reset();
1907 int nb = hent->GetNbinsX()+1;
1908 for (int ib=0;ib<=nb;ib++) hent->SetBinContent(ib, prof->GetBinEntries(ib));
1909 return hent;
1910}
1911
1912//___________________________________________________________
1913TH1* FitDXEdges(TProfile* prof)
1914{
1915 //
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();
1927 TString fstatus;
1928 //
1929 TH1* histo = ProfileAsTH1(prof,"_h1");
1930 //
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;
1937 //
1938 // RedoProfileErrors(histo,prof);
1939 // find left/right significant bins
1940 int bleft = -1;
1941 for (int i=1;i<=nb;i++) {
1942 double enti = prof->GetBinEntries(i);
1943 if (enti>kMinThresh*meanbin) {bleft=i; break;}
1944 }
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;
1949 //
1950 int bright = -1;
1951 for (int i=nb;i>0;i--) {
1952 double enti = prof->GetBinEntries(i);
1953 if (enti>kMinThresh*meanbin) {bright=i; break;}
1954 }
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;
1959 //
1960 printf("Fit left edge\n");
1961 if (!flft) flft = new TF1("lftEdge",EdgeFun,
1962 prof->GetBinLowEdge(1),
1963 prof->GetBinLowEdge(nb+1),
1964 9);
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);
1978 int cntL = 0;
1979 double chiL = 0;
1980 do {
1981 prof->Fit(flft,"q+","",fitLStart,fitLEnd);
1982 fstatus = (char*)gMinuit->fCstatu.Data();
1983 if (fstatus.Contains("CONVERGED") || fstatus.Contains("SUCCESSFUL")) {
1984 cntL=100;
1985 chiL = flft->GetNDF()>0 ? flft->GetChisquare()/flft->GetNDF() : 0;
1986 }
1987 else {
1988 TF1* oldf = (TF1*)prof->GetListOfFunctions()->FindObject(flft->GetName());
1989 if (oldf) prof->GetListOfFunctions()->Remove(oldf);
1990 }
1991 } while(++cntL<3);
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);
2004 }
2005 histo->SetBinContent(0, TMath::Max(0.0, flft->GetParameter(1)-6*flft->GetParameter(0))); // effective lower sensor edge
2006 }
2007 else {
2008 printf("Left edge bad: %f %d %s\n",chiL,cntL,fstatus.Data());
2009 }
2010 //
2011 printf("Fit right edge\n");
2012 if (!frgt) frgt = new TF1("rgtEdge",EdgeFun,
2013 prof->GetBinLowEdge(1),
2014 prof->GetBinLowEdge(nb+1),
2015 9);
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);
2029 double chiR = 0;
2030 int cntR = 0;
2031 do {
2032 prof->Fit(frgt,"q+","",fitRStart,fitREnd);
2033 fstatus = (char*)gMinuit->fCstatu.Data();
2034 if (fstatus.Contains("CONVERGED") || fstatus.Contains("SUCCESSFUL")) {
2035 cntR=100;
2036 chiR = frgt->GetNDF()>0 ? frgt->GetChisquare()/frgt->GetNDF() : 0;
2037 }
2038 else {
2039 TF1* oldf = (TF1*)prof->GetListOfFunctions()->FindObject(frgt->GetName());
2040 if (oldf) prof->GetListOfFunctions()->Remove(oldf);
2041 }
2042 } while((++cntR<3));
2043 //
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);
2056 }
2057 histo->SetBinContent(nb+1, TMath::Min((double)sddSeg.Dx(), frgt->GetParameter(2)+6*frgt->GetParameter(0))); // effective upper sensor edge
2058 }
2059 else {
2060 printf("Right edge bad: %f %d %s\n",chiR,cntR,fstatus.Data());
2061 }
2062 //
2063 return histo;
2064}
2065
2066double ftPolComb(double* x, double *par)
2067{
2068 // fit with combination of 2 polinomials of order par[0]
2069 int ord = int(par[0]);
2070 int npars = ord+1;
2071 double brk = par[1];
2072 //
2073 double px = x[0];
2074 double res = 0;
2075 int start = 2;
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];
2080 }
2081 return res;
2082}
2083
2084//______________________________________________________________
2085TH1* SimpleMap(TH1* prof)
2086{
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);
2091 //
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());
2096 //
2097 while(1) {
2098 //
2099 TF1* smapf = 0;
2100 double mean = 0;
2101 // 1) try pol1
2102 smapf = new TF1("smapf","pol1",prof->GetBinLowEdge(0),prof->GetBinLowEdge(nb+1));
2103 if (TestMapFunction(prof,smapf,lft,rgt)) {delete smapf; break;}
2104 else {
2105 mean = smapf->GetParameter(0);
2106 delete smapf;
2107 }
2108 //
2109 // 2) try pol2
2110 smapf = new TF1("smapf","pol2",prof->GetBinLowEdge(0),prof->GetBinLowEdge(nb+1));
2111 if (TestMapFunction(prof,smapf,lft,rgt)) {delete smapf; break;}
2112 else delete smapf;
2113 //
2114 /*
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);
2121 //
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
2126 //
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;}
2133 else delete smapf;
2134 //
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;}
2141 else delete smapf;
2142 //
2143 */
2144 break;
2145 }
2146 TF1* fnsel = (TF1*) prof->GetListOfFunctions()->FindObject("smapf");
2147 if (!fnsel) {delete smap; return 0;} // no simple solution
2148 //
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));
2155 smap->Reset();
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);
2160 }
2161 return smap;
2162}
2163
2164
2165//______________________________________________________________________
2166Bool_t TestMapFunction(TH1* smap, TF1* fun, double lft, double rgt)
2167{
2168 // test if fun describes the shape
2169 TString fstatus;
2170 double chi2;
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);
2179 return kTRUE;
2180 }
2181 }
2182 return kFALSE;
2183}