]> 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
fHits size reduced to previous assignment
[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>
15#include "AliCDBManager.h"
16#include "AliCDBEntry.h"
17#include "AliCDBStorage.h"
18#include "AliITSsegmentationSDD.h"
19#include "AliITSCorrMapSDD.h"
20#include "AliITSCorrMap1DSDD.h"
21#include "AliITSDriftSpeedArraySDD.h"
22#include "AliITSDriftSpeedSDD.h"
23#include "AliITSresponseSDD.h"
24#endif
25
26/* $Id$ */
27
28/*
29 // master equations
30 DeltaX_L = -deltaX + deltaV_L*(T_L - T0_true - deltaT0) - V_Ltrue*deltaT0 - corrMapL
31 DeltaX_R = deltaX + deltaV_R*(T_R - T0_true - deltaT0) - V_Rtrue*deltaT0 - corrMapR
32 //
33
34 where
35 corrMapL(R) : Vdrift vs X nonuniformity correction maps
36 V_L(R)true: true drift speed
37 deltaV_L(R) : error in drift speed
38 V_Ltrue = V_Lassumed - deltaV_L
39 V_Rtrue = V_Rassumed - deltaV_R
40 //
41 DeltaX : measured residuals (track extrapolation - measurement) for left (L)
42 and right (R) sides (right side multiplied by -1!!!)
43 deltaX : misalignment of module (Xtrue = Xassumed - deltaX)
44 T_L(R) : measured drift time (w/o T0 subtraction)
45 T0_true : real T0
46 deltaT0 : error in T0
47*/
48
49//------------------- Details of CDB objects to create ---------------
50int firstRun = -1;
51int lastRun = -1;
52TString cdbComment = "";
53//--------------------------------------------------------------------
54
55const Int_t kSDDMin=240,kSDDMax=499,kNSDD=kSDDMax-kSDDMin+1;
56enum {kSXvsX,kSXvsZ,kSTDvsX,kSTDvsZ,kSDXvsXclean,kSVDvxX,kSCorrMapX,kSXvsXCorr,kSXvsZCorr,kSDVvsZ,kSDVvsZOrig,kSGrDxx,kSGrTDx,kSGrTXCorr, kNStore};
57const Int_t maxVZOrd = 4; // max order of polinomial for VD vs Z correction (normal sensors)
58const Double_t chi2BadVZ = 2.9; // treat as bad module if chi2/ndf is worse
59//
60enum {kDxX,kDxZ,kTDX,kTDZ,kNQATypes};
61// Prefixes for residuals histos prepared by PrepSDDResidMap.C from MP2 output
62// For QA histograms the name should be changed to "hpSDDResXvsXD","hpSDDResXvsZ","hpSDDDrTimevsXD","hpSDDDrTimevsZ" respectively
63const Char_t* kQAHistoNames[] = { //
64 "Xdrift_Mod_", // DX vs Xdrift
65 "ZAnode_Mod_", // DX vs local Z
66 "TD_vs_Xdrift_Mod_", // TDrift-T0 vs Xdrift
67 "TD_vs_ZAnode_Mod_" // TDrift-T0 vs local Z
68};
69const Char_t* kLRNames[] = {"left","right"}; // identifiers for left/right sides in histo names
70//
71double minBEntryX = 10; // min entries in the bin to account (vs Xdrift)
72double minBEntryZ = 20; // min entries in the bin to account (vs Z)
73double skipDXEdge = 1000.; // microns to skip from each side of XDrift when smoothing
74double edgeSmearMinErr = 60.; // minimal error for extrapolation
75double edgeSmearMaxErr = 600.; // maximal error for extrapolation
76double wDXEdge = 5000.; // width in microns for XDrift edges to average
77double wVDEdge = 2500; // width in microns for VDrift averaging region on the edges
78double threshClean= 0.3; // clean edge bins with occupancy smaller than in the neighbourhood
79int nVDEdge = 20; // min width in n bins for VDrift averaging region on the edges
80int smoothLevel= 5; // smoothing level
81int minBinsEdge = 5; // min number of X>0 good bins to fit the edge
82Bool_t useSplieForR2V = kFALSE; // if true, extract Vd profile using Derivateve of spline (may be inadequate for sharp turns)
83//
84Bool_t forceT0CorrTo0 = kTRUE;
85Bool_t forceRPhiCorrTo0 = kTRUE;
86//
87//
88AliITSsegmentationSDD sddSeg;
89//
90//----------------- working variables --------------------------------
91Int_t currSDDId=-1,currSDDSide=-1,currMod=-1;
92TProfile* currHDXvsX = 0; // raw DX vs X histo from QA
93TProfile* currHDXvsZ = 0; // raw DX vs Z histo from QA
94TProfile* currHTDvsX = 0; // raw DriftTime-T0 vs X histo from QA
95TProfile* currHTDvsZ = 0; // raw DriftTime-T0 vs Z histo from QA
96TH1* currHDXvsXclean = 0; // smoothed DX vs X histo
97TH1* currHVDvxX = 0; // extracted VDrift vs X profile
98TH1* currHCorrMapX = 0; // extracted correction map
99TH1* currHDXvsXCorr = 0; // DX vs X histo after application of the map
100TH1* currHDXvsZCorrMT0 = 0; // DX vs X histo after application of the map and t0 correction
101TH1* currHDVvsZ = 0; // correction for the VDrift vs Z
102TH1* currHDVvsZOrig = 0; // correction for the VDrift vs Z (before error processing)
103TGraphErrors* currGrDxx = 0; // DX final vs XDrift
104TGraphErrors* currGrTDx = 0; // DX final vs TDrift
105TGraphErrors* currGrTXCorr = 0; // TDrift vs XDritf
106TCollection* qaHistos = 0; // input QA collection
107TObjArray procHistos; // processed histos buffer
108AliITSresponseSDD* sddResp=0; // SDD response, updated
109TObjArray *vdarrayOld = 0; // olf VDrifts array
110TObjArray* vDriftArr=0; // VDrifts array, updated
111TObjArray* corrMaps=0; // holder for correction maps
112//
113TH1* resOffsDXraw[2]={0}; // Dx vs X offset at Xdrift=0 raw
114TH1* resOffsDXAMap[2]={0}; // Dx vs X offset (mean) after correction by map
115TH1* resOffsDX[2]={0}; // Dx vs X offset at Xdrift=0 after correction by map (for each side)
116TH1* resVDCorr[2]={0}; // VDrift correction (for each side)
117TH1* resVDMean[2]={0}; // average VDrift (for each side)
118TH1* resVDCorrZ[2]={0}; // mean VDrift vs Z corrections (for each side)
119TH1* resT0Corr = 0; // correction to modules T0
120TH1* resXLocCorr = 0; // correction to module location
121TCanvas* sddCanv = 0; // report canvas
122//
123Bool_t sddRespUpdT0OK = kFALSE; // flag that SDDresponse object T0 was updated
124Bool_t sddRespUpdVDOK = kFALSE; // flag that SDDresponse object VDrift correction was updated
125Bool_t sddVDriftUpdOK = kFALSE; // flag that SDD VDrift object was updated
126//
127TString pathSDDRespOld = "";
128TString pathSDDVDriftOld = "";
129TString pathSDDCorrMapOld = "";
130//
131//--------------------------------------------------------------------
132void Process(const char* pathSDDResp=0,const char* pathSDDVDrift=0, const char* pathSDDCorrMap=0);
133void Process(TCollection* qa, const char* pathSDDResp=0,const char* pathSDDVDrift=0, const char* pathSDDCorrMap=0);
134Bool_t ProcSDDSideDXvsX();
135Bool_t ProcSDDSideDXvsZ();
136TProfile* GetQAHisto(Int_t qaType);
137TProfile* H2Profile(TH2* h2);
138TH1* H2ProfileAsTH1(TH2* h2);
139TH1* ProfileAsTH1(TProfile* prof, const char* addName);
140TH1* CleanProfile(TProfile* profR);
141TH1* Vdrift2Resid(TH1* vd);
142TH1* Resid2Vdrift(TH1* res);
143void RedoProfileErrors(TH1* profH1,TProfile* prof);
144void SetModuleID(Int_t mdID,Int_t side=-1);
145void PrepareModuleHistos(Int_t mdID, Int_t side, Bool_t vsZ);
146void CheckResultDX();
147void CalcDXCorrections();
148double GetVOld(double z);
149Int_t GetStoreID(int type, int imd=-1,int side=-1);
150void CleanPrev();
151void StoreCurrent();
152double ftVdZ(double *x, double *par);
153void PlotReport(const char* psname="repSDDQA.ps");
154TH1* GetPadBaseHisto(TPad* pad);
155Bool_t PlotHisto(TH1* h, Option_t* opt="", int mrkStyle=20,int mrkCol=kBlack, double mrkSize=1.);
156TLatex* AddPadLabel(const char*txt,float x=0.1,float y=0.9,int color=kBlack,float size=0.04);
157void GetHistoMinMaxInRange(TH1* h, double &mn,double &mx);
158double edgeLow(double* x, double *par);
159double edgeHigh(double* x, double *par);
160void CureEdges(TH1* prof);
161//
162Bool_t LoadSDDVDrift(TString& path, TObjArray *&arr);
163Bool_t LoadSDDResponse(TString& path, AliITSresponseSDD *&resp);
164Bool_t LoadSDDCorrMap(TString& path, TObjArray *&maps);
165AliCDBEntry* GetCDBEntry(const char* path);
166//
167void UpdateSDDResponse(AliITSresponseSDD *resp, Bool_t t0=kTRUE, Bool_t vdrift=kTRUE);
168void UpdateSDDVDrift(AliITSDriftSpeedArraySDD* vdArr, int imd, int side);
169TObjArray* UpdateSDDVDrift();
170TObjArray* CreateCorrMaps();
171AliITSresponseSDD* UpdateSDDResponse(Bool_t t0=kTRUE, Bool_t vdrift=kTRUE);
172AliITSCorrMap1DSDD* CreateCorrMap(TH1* mapHisto, int imd, int side, AliITSCorrMap1DSDD* updateMap=0);
173void PrepCDBObj(TObject *obj,const char* path,int firstrun=0,int lastrun=999999999,const char* comment="");
174
175//-------------------------------------------------------------------
176
177//_________________________________________________________________________
178void Process(TCollection* qa, const char* pathSDDResp, const char* pathSDDVDrift, const char* pathSDDCorrMap)
179{
180 // process all
181 qaHistos = qa;
182 Process(pathSDDResp,pathSDDVDrift,pathSDDCorrMap);
183}
184
185//_________________________________________________________________________
186void Process(const char* pathSDDResp, const char* pathSDDVDrift, const char* pathSDDCorrMap)
187{
188 // process all
189 procHistos.Clear();
190 pathSDDRespOld = pathSDDResp;
191 pathSDDVDriftOld = pathSDDVDrift;
192 pathSDDCorrMapOld = pathSDDCorrMap;
193 sddRespUpdT0OK = sddRespUpdVDOK = sddVDriftUpdOK = kFALSE;
194 //
195 if (pathSDDVDriftOld.IsNull()) printf("Attention: Old VDrift is missing!\n");
196 //
197 for (int imd=kSDDMin;imd<=kSDDMax;imd++) {
198 for (int ix=0;ix<2;ix++) {
199 CleanPrev(); // clean data from previous module
200 printf("Processing %d/%d\n",imd,ix);
201 PrepareModuleHistos(imd,ix, kTRUE);
202 ProcSDDSideDXvsX(); // process DX vs X histos to get corr.maps, deltaT0, deltaX0, deltaV_mean
203 StoreCurrent();
204 }
205 CalcDXCorrections(); // calculate delta's
206 //
207 for (int ix=0;ix<2;ix++) {
208 CleanPrev(); // clean data from previous module
209 PrepareModuleHistos(imd,ix, kFALSE);
210 ProcSDDSideDXvsZ(); // process deltaV vs anode profile
211 StoreCurrent();
212 }
213 }
214 //
215 corrMaps = CreateCorrMaps(); // create correction maps
216 PrepCDBObj(corrMaps,"ITS/Calib/MapsTimeSDD",firstRun,lastRun,cdbComment.Data());
217 //
218 if (!pathSDDVDriftOld.IsNull()) {
219 vDriftArr = UpdateSDDVDrift();
220 PrepCDBObj(vDriftArr,"ITS/Calib/DriftSpeedSDD",firstRun,lastRun,cdbComment.Data());
221 }
222 //
223 if (!pathSDDRespOld.IsNull()) {
224 sddResp = UpdateSDDResponse();
225 PrepCDBObj(sddResp,"ITS/Calib/RespSDD",firstRun,lastRun,cdbComment.Data());
226 }
227 //
228 PlotReport();
229}
230
231//_________________________________________________________________________
232Bool_t ProcSDDSideDXvsX()
233{
234 // process one side of the SDD module
235 currHDXvsXCorr = ProfileAsTH1(currHDXvsX,"corrCheck");
236 RedoProfileErrors(currHDXvsXCorr,currHDXvsX);
237 //
238 if ( (currHDXvsXclean=CleanProfile(currHDXvsX)) ) {
239 //
240 delete currHDXvsXCorr;
241 currHDXvsXCorr = (TH1*)currHDXvsXclean->Clone( Form("%s_corrCheck",currHDXvsX->GetName()) );
242 //
243 currHVDvxX = Resid2Vdrift(currHDXvsXclean);
244 currHCorrMapX = Vdrift2Resid(currHVDvxX);
245 //
246 currHDXvsXCorr->Add(currHCorrMapX,-1);
247 }
248 // check results
249 CheckResultDX();
250 //
251 return kTRUE;
252}
253
254//_________________________________________________________________________________
255Bool_t ProcSDDSideDXvsZ()
256{
257 // extract correction for Vdrift vs Z from DXvsZ profile and mean_drift_time vs Z
258 //
259 static TF1* fitP0 = new TF1("fitP0","pol0",-4000,4000);
260 double chi2s[maxVZOrd+1] = {0};
261 static TF1* fitVvsZs[maxVZOrd+1] = {0};
262 static Bool_t iniDone = kFALSE;
263 double zrange = sddSeg.Dz()/2 - 1.;
264 if (!iniDone) {
265 iniDone = kTRUE;
266 for (int iord=0;iord<=maxVZOrd;iord++) {
267 fitVvsZs[iord] = new TF1("fitVvsZ",ftVdZ, -zrange, zrange ,iord+2);
268 fitVvsZs[iord]->FixParameter(0, iord+0.1);
269 }
270 }
271 //
272 int nb = currHDXvsZ->GetNbinsX();
273 if (currHDXvsZ->GetEntries()<minBEntryZ*nb) return kFALSE;
274 //
275 currHDVvsZ = ProfileAsTH1(currHDXvsZ,"_vzcorr");
276 currHDVvsZ->Reset();
277 //
278 currHDXvsZCorrMT0 = ProfileAsTH1(currHDXvsZ,"_zcorrMapT0");
279 currHDXvsZCorrMT0->Reset();
280 //
281 int nbUse=0, ib0 = currHDVvsZ->FindBin(-zrange), ib1 = currHDVvsZ->FindBin( zrange);
282 double *statZB = new double[nb+1]; // entries per Z bin
283 memset(statZB,0,sizeof(double)*(nb+1));
284 double vmean=0, vmeanE=0, norm=0;
285 //
286 for (int ib=ib0;ib<=ib1;ib++) {
287 double ne = currHTDvsZ->GetBinEntries(ib);
288 if (ne<1) continue;
289 double dx = currHDXvsZ->GetBinContent(ib); // raw residual X vs Z
290 double dxe= currHDXvsZ->GetBinError(ib);
291 double t = currHTDvsZ->GetBinContent(ib); // mean assumed (TDrift - T0)
292 double te = currHTDvsZ->GetBinError(ib);
293 double vCorr = resVDCorr[currSDDSide]->GetBinContent(currMod+1);
294 dx -= vCorr*t; // subtract mean V correction
295 dxe = TMath::Sqrt(dxe*dxe + vCorr*vCorr*te*te);
296 currHDXvsZCorrMT0->SetBinContent(ib,dx);
297 currHDXvsZCorrMT0->SetBinError(ib,dxe);
298 }
299 currHDXvsZCorrMT0->Fit(fitP0,"q0","");
300 double pedestal = fitP0->GetParameter(0);
301 //
302 for (int ib=ib0;ib<=ib1;ib++) {
303 double ne = currHTDvsZ->GetBinEntries(ib);
304 if (ne<minBEntryZ) continue;
305 double dx = currHDXvsZCorrMT0->GetBinContent(ib); // raw residual X vs Z
306 double dxe= currHDXvsZCorrMT0->GetBinError(ib);
307 double t = currHTDvsZ->GetBinContent(ib); // mean assumed (TDrift - T0)
308 double te = currHTDvsZ->GetBinError(ib);
309 if (t<1) continue;
310 // corrections
311 // account for the corrections leading the mean DX to 0
312 dx -= pedestal;
313 double v = dx/t;
314 double ve = TMath::Sqrt(dxe*dxe + v*v*te*te)/t;
315 //
316 nbUse++;
317 vmeanE += 1./(ve*ve);
318 vmean += v/(ve*ve);
319 // printf("%d %f %f | %f %f | -> %f %f\n",ib,dx,dxe,t,te,v,ve);
320 int ibDest = currSDDSide==0 ? ib : nb+1-ib;
321 statZB[ibDest] = ne;
322 currHDVvsZ->SetBinContent(ibDest,v);
323 currHDVvsZ->SetBinError(ibDest,ve);
324 //
325 }
326 //
327 if (nbUse<maxVZOrd+1) {delete statZB; return kFALSE;}
328 if (vmeanE>0) { vmean /= vmeanE; vmeanE = 1./TMath::Sqrt(vmeanE); }
329 // fill empty and bad bins by mean value with large error
330 vmeanE *= nbUse;
331 /*
332 for (int ib=ib0;ib<=ib1;ib++) {
333 //if (statZB[ib]<minBEntryZ) currHDVvsZ->SetBinContent(ib,vmean);
334 if (currHDVvsZ->GetBinError(ib)>vmeanE || statZB[ib]<minBEntryZ) {
335 currHDVvsZ->SetBinError(ib,vmeanE*2);
336 currHDVvsZ->SetBinContent(ib,vmean);
337 }
338 }
339 */
340 // save original
341 currHDVvsZOrig = (TH1*) currHDVvsZ->Clone(Form("%s_orig",currHDVvsZ->GetName()));
342 // find leftmost good bin
343 int bL,bR;
344 double errL=0,errR=0, valL=0,valR=0;
345 for (bL=ib0;bL<=ib1;bL++) {
346 if (/*currHDVvsZ->GetBinError(bL)>vmeanE ||*/ statZB[bL]<minBEntryZ) continue;
347 errL = currHDVvsZ->GetBinError(bL);
348 valL = currHDVvsZ->GetBinContent(bL);
349 break;
350 }
351 for (bR=ib1;bR>=ib0;bR--) {
352 if (/*currHDVvsZ->GetBinError(bR)>vmeanE ||*/ statZB[bR]<minBEntryZ) continue;
353 errR = currHDVvsZ->GetBinError(bR);
354 valR = currHDVvsZ->GetBinContent(bR);
355 break;
356 }
357 // change bad bins
358 for (int ib=bL-1;ib>=ib0;ib--) {
359 double err = errL + (vmeanE-errL)*float(ib-ib0)/(bL-ib0+1);
360 double val = vmean + (valL-vmean)*float(ib-ib0)/(bL-ib0+1);
361 currHDVvsZ->SetBinError(ib,err);
362 currHDVvsZ->SetBinContent(ib,val);//currHDVvsZ->GetBinContent(ib+1));
363 }
364 for (int ib=bR+1;ib<=ib1;ib++) {
365 double err = errR + (vmeanE-errR)*float(ib1-ib)/(ib1-bR+1);
366 double val = vmean + (valR-vmean)*float(ib1-ib)/(ib1-bR+1);
367 currHDVvsZ->SetBinError(ib,err);
368 currHDVvsZ->SetBinContent(ib,val);//currHDVvsZ->GetBinContent(ib-1));
369 }
370 //
371 for (int iord=0;iord<=maxVZOrd;iord++) {
372 currHDVvsZ->Fit(fitVvsZs[iord],"q0");
373 chi2s[iord] = fitVvsZs[iord]->GetChisquare();
374 int ndf = fitVvsZs[iord]->GetNDF();
375 if (ndf>0) chi2s[iord] /= ndf;
376 }
377 //
378 // analyse chi2's
379 int bestOrd = 0;
380 double bestChi2 = 9999;
381 for (int iOrd=0;iOrd<=maxVZOrd;iOrd++) {
382 if (chi2s[iOrd]<bestChi2-0.5) {bestChi2 = chi2s[iOrd]; bestOrd = iOrd;}
383 }
384 TF1* fitVvsZ = fitVvsZs[bestOrd];
385 currHDVvsZ->Fit(fitVvsZ,"q");
386 //
387 // extract mean correction neglecting the constant term
388 vmean = 0;
389 double freePar = 0;//fitVvsZ->GetParameter(1);
390 for (int ib=1;ib<=nb;ib++) {
391 if (statZB[ib]<1) continue;
392 vmean += statZB[ib]*(fitVvsZ->Eval(currHDVvsZ->GetBinCenter(ib)) - freePar); // account only Z dependent part
393 norm += statZB[ib];
394 }
395 //
396 if (!resVDCorrZ[currSDDSide]) {
397 resVDCorrZ[currSDDSide] = new TH1F(Form("VDcorrvsZ%d",currSDDSide),Form("mean VDrift correction vs Z, side %d",currSDDSide),kNSDD,kSDDMin-0.5,kSDDMax+0.5);
398 resVDCorrZ[currSDDSide]->SetMarkerColor(2+currSDDSide);
399 resVDCorrZ[currSDDSide]->SetMarkerStyle(20+4*currSDDSide);
400 }
401 //
402 // if (currSDDSide) vmean = -vmean;
403 resVDCorrZ[currSDDSide]->SetBinContent(currMod+1, norm>0 ? vmean/norm : 0);
404 delete statZB;
405 return kTRUE;
406}
407
408//_________________________________________________________________________
409void PrepareModuleHistos(Int_t mdID, Int_t side, Bool_t vsX)
410{
411 // retrieve QA histos, if needed, convert to TH1
412 SetModuleID(mdID, side);
413 if (vsX) {
414 currHDXvsX = GetQAHisto(kDxX);
415 currHTDvsX = GetQAHisto(kTDX);
416 }
417 else {
418 currHTDvsZ = GetQAHisto(kTDZ);
419 currHDXvsZ = GetQAHisto(kDxZ);
420 }
421}
422
423//_________________________________________________________________________
424void SetModuleID(Int_t mdID,Int_t side)
425{
426 // set current module ID,side
427 currSDDId = mdID;
428 currSDDSide = side;
429 currMod = currSDDId - kSDDMin;
430}
431
432//_________________________________________________________________________
433TProfile* GetQAHisto(Int_t qaType)
434{
435 // retrieve needed QA histo
436 if (!qaHistos) {printf("QA Histos collection is not set\n"); exit(1);}
437 //
438 if (currSDDId<kSDDMin || currSDDId>kSDDMax || currSDDSide<0 || currSDDSide>2)
439 {printf("Illegal SDD module ID/Side: %d/%d\n",currSDDId,currSDDSide); exit(1);}
440 //
441 if (qaType<0 || qaType>=kNQATypes)
442 {printf("Illegal QA Histo Type %d\n",qaType); exit(1);}
443 //
444 const char* hname = Form("%s%d_%s",kQAHistoNames[qaType],currSDDId,kLRNames[currSDDSide]);
445 TH1* h = (TH1*) qaHistos->FindObject(hname);
446 if (!h) {printf("Did not find histo %s in QA output collection\n",hname); exit(1);}
447 //
448 if (h->InheritsFrom(TH2::Class())) h = H2Profile((TH2*)h); // make sure it is not TH2 histo
449 //
450 return (TProfile*)h;
451}
452
453//_________________________________________________________________________
454TH1* H2ProfileAsTH1(TH2* h2)
455{
456 // extract profile as TH1 histo
457 TProfile* prf = h2->ProfileX();
458 TH1* prof = ProfileAsTH1(prf, "_profH1");
459 delete prf;
460 return prof;
461}
462
463//_________________________________________________________________________
464TH1* ProfileAsTH1(TProfile* prf, const char* addName)
465{
466 // convert profile to TH1 histo
467 TString nm = prf->GetName(); nm += addName;
468 TAxis* xax = prf->GetXaxis();
469 TH1* prof = new TH1F(nm.Data(),nm.Data(),prf->GetNbinsX(),xax->GetXmin(),xax->GetXmax());
470 for (int i=1;i<=prof->GetNbinsX();i++) {
471 prof->SetBinContent(i, prf->GetBinContent(i));
472 prof->SetBinError(i, prf->GetBinError(i));
473 }
474 return prof;
475}
476
477//_________________________________________________________________________
478TProfile* H2Profile(TH2* h2)
479{
480 // extract profile
481 TString nm = h2->GetName(); nm += "_prof";
482 TProfile* prf = h2->ProfileX(nm.Data());
483 return prf;
484}
485
486//_________________________________________________________________________________
487TH1* CleanProfile(TProfile* profR)
488{
489 // clean profile copy
490 int ib0 = profR->FindBin(1), ib1 = profR->FindBin(sddSeg.Dx()-1);
491 int nbn = profR->GetNbinsX();
492 int nSkip = skipDXEdge/profR->GetBinWidth(1);
493 int nMean = wDXEdge/profR->GetBinWidth(1);
494 //
495 TH1* prof = ProfileAsTH1(profR,"_Clean");
496 RedoProfileErrors(prof,profR);
497 //
498 CureEdges(prof);
499 //
500 // for (int i=1;i<ib0;i++) {prof->SetBinContent(i,0); prof->SetBinError(i,0);}
501 // for (int i=ib1+1;i<=nbn;i++) {prof->SetBinContent(i,0); prof->SetBinError(i,0);}
502 // get mean occupancy skipping first and last nSkip occupied bins
503 //
504 int nbCntL=0, nbCntR=0,nskipL=0,nskipR=0;
505 double smL=0,smR=0;
506 for (int ib=ib0;ib<=ib1 && nbCntL<nMean;ib++) { // skipping left, get average of good bins
507 double vl = profR->GetBinEntries(ib);
508 if (vl<minBEntryX) continue;
509 if (nskipL<nSkip) {nskipL++; continue;}
510 smL += vl;
511 nbCntL++;
512 }
513 if (nbCntL<1) return 0;
514 smL /= nbCntL;
515 //
516 for (int ib=ib1+1;ib>ib0 && nbCntR<nMean;ib--) { // skipping right
517 double vl = profR->GetBinEntries(ib);
518 if (vl<minBEntryX) continue;
519 if (nskipR<nSkip) {nskipR++; continue;}
520 smR += vl;
521 nbCntR++;
522 }
523 //
524 if (nbCntR<1) return 0;
525 smR /= nbCntR;
526 //
527 prof->GetXaxis()->SetRange(ib0,ib1);
528 prof->Smooth(smoothLevel,"r");
529 prof->GetXaxis()->SetRange(1,nbn);
530 // kill bins with small statistics, from left and right
531 for (int ib=1;ib<ib1;ib++) if (profR->GetBinEntries(ib)<threshClean*smL) {prof->SetBinContent(ib,0);prof->SetBinError(ib,0);} else break;
532 for (int ib=ib1;ib>ib0;ib--) if (profR->GetBinEntries(ib)<threshClean*smR) {prof->SetBinContent(ib,0);prof->SetBinError(ib,0);} else break;
533 //
534 return prof;
535}
536
537
538//_________________________________________________________________________________
539TH1* Resid2Vdrift(TH1* res)
540{
541 // extract Vdrift profile from residuals profile
542 //
543 TSpline3 spl(res);
544 TString nm = res->GetName(); nm += "_vd";
545 TH1* vd = (TH1*) res->Clone(nm.Data());
546 vd->Reset();
547 int nb = vd->GetNbinsX();
548 //
549 int nbcnt=0, nbfl = wVDEdge/vd->GetBinWidth(1);
550 if (nbfl<nVDEdge) nbfl=nVDEdge;
551 //
552 double vAv=0,eAv=0; // average vdrift / v0
553 for (int i=1;i<=nb;i++) {
554 double deriv = useSplieForR2V ? spl.Derivative(vd->GetBinCenter(i)) : (res->GetBinContent(i)-res->GetBinContent(i-1))/res->GetBinWidth(i);
555 double v = TMath::Abs(deriv-1)>1e-6 ? 1./(1-deriv) : 1.;
556 vd->SetBinContent(i, v);
557 double err2 = 0; // set relative error
558 int nbd = 1;
559 for (int i1=i-1;i1<=i+1;i1++) {
560 if (i1<1 || i1>nb) continue;
561 double err = res->GetBinError(i);
562 if (err<1e-9) err = 1e4;
563 err2 += err*err;
564 nbd++;
565 }
566 if (nbd) err2/=nbd;
567 if (err2>0 && v>0 && v<3) {vAv += v/err2; eAv += 1./err2;}
568 err2 = TMath::Sqrt(err2)/res->GetBinWidth(1);
569 vd->SetBinError(i, err2 );
570 nbcnt++;
571 }
572 vAv = eAv>0 ? vAv/eAv : 1.0;
573 // printf("mean V = %f in %d bins\n",vAv,nbcnt);
574 // cure anomalous bins
575 // for (int i=1;i<=nb;i++) if (vd->GetBinContent(i)<0.1) vd->SetBinContent(i,vAv);
576 //
577 int ib = 0;
578 for (ib=1;ib<nb;ib++) { // detect up to which bin the left tail is unstable
579 double x = vd->GetBinCenter(ib);
580 if (x<0 || res->GetBinError(ib)<1e-9) {vd->SetBinContent(ib,1); vd->SetBinError(ib,0); continue;}
581 double vl = vd->GetBinContent(ib);
582 if (TMath::Abs(vl-vAv)<0.2) break;
583 vd->SetBinContent(ib,1);
584 vd->SetBinError(ib,0);
585 }
586 //
587 double vL=0,eL=0,vR=0,eR=0;
588 // get the mean of leftmost stable vdrift
589 int lastCounted = 0;
590 nbcnt = 0;
591 for (int i=ib;i<=nb && nbcnt<nbfl;i++) {
592 double berr = vd->GetBinError(i);
593 if (berr<1e-9 || berr>0.5 || TMath::Abs(vd->GetBinContent(i)-vAv)>0.2) continue;
594 vL += vd->GetBinContent(i)/berr/berr;
595 eL += 1./berr/berr;
596 nbcnt++;
597 lastCounted = i;
598 }
599 vL = eL>0 ? vL/eL : vAv;
600 //printf("VLeft: %f, in %d bins (%d %d)\n",vL,nbcnt,ib,lastCounted);
601 for (int i=1;i<=ib;i++) vd->SetBinContent(i, vL);
602 // for safety check if there are no outliers in first 3 "stable" bins
603 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);
604 vd->SetBinContent(vd->FindBin(1),1.); // no correction at t=0 !!
605 //
606 double lmax = sddSeg.Dx()-1;
607 for (ib=nb+1;ib--;) { // detect up to which bin the right tail is unstable
608 double x = vd->GetBinCenter(ib);
609 if (x>=lmax || res->GetBinError(ib)<1e-9) {vd->SetBinContent(ib,1);vd->SetBinError(ib,0); continue;}
610 if (TMath::Abs(vd->GetBinContent(ib)-vAv)<0.4) break;
611 vd->SetBinContent(ib,vAv);
612 vd->SetBinError(ib,0);
613 }
614 nbcnt= 0;
615 lastCounted = 0;
616 for (int i=ib;i>=1 && nbcnt<nbfl;i--) {
617 double berr = vd->GetBinError(i);
618 if (berr<1e-9 || berr>0.5 || TMath::Abs(vd->GetBinContent(i)-vAv)>0.4) continue;
619 vR += vd->GetBinContent(i)/berr/berr;
620 eR += 1./berr/berr;
621 nbcnt++;
622 lastCounted = i;
623 }
624 vR = eR>0 ? vR/eR : vAv;
625 // printf("VRight: %f, in %d bins (%d %d)\n",vR,nbcnt,lastCounted,ib);
626 for (int i=ib;i<=nb;i++) vd->SetBinContent(i, vR);
627 // for safety check if there are no outliers in first 3 "stable bins
628 for (int i=(lastCounted+ib)/2;i<ib;i++)
629 if ( vd->GetBinError(i)<1e-9 || (vd->GetBinError(i)>0.3 && TMath::Abs(vd->GetBinContent(i)-vR)>0.2)) vd->SetBinContent(i, vR);
630 // fit the empty bins on the right
631 // vd->Fit(fitP1,"+","",vd->GetBinCenter(ib-nbfl),vd->GetBinCenter(nb));
632 //for (int i=ib;i<=nb;i++) vd->SetBinContent(i, fitP1->Eval(vd->GetBinCenter(i)));
633 //
634 return vd;
635}
636
637
638//_________________________________________________________________________________
639TH1* Vdrift2Resid(TH1* vd)
640{
641 // convert Vdrift profile to residuals profile (final correction map)
642 //
643 TString nm = vd->GetName(); nm += "_vd2res";
644 TH1* res = (TH1*) vd->Clone(nm.Data());
645 res->Reset();
646 int ib0 = res->FindBin(1);
647 int ib1 = res->FindBin(sddSeg.Dx()-1);
648 double resv=0;
649 int lastB=0;
650 for (lastB=ib1+1;lastB--;) if (vd->GetBinError(lastB)>1e-9) break; // find last good bin
651 double lastX = vd->GetBinCenter(lastB);
652 // extend by 1mm
653 lastX += 1000;
654 if (lastX>sddSeg.Dx()) lastX = sddSeg.Dx();
655 lastB = vd->FindBin(lastX);
656 //
657 // 1st iteration : estimate correction at max Xdrift
658 for (int i=ib0;i<=lastB;i++) {
659 double dx = res->GetBinWidth(i);
660 double v = vd->GetBinContent(i);
661 resv += dx*(1.-1./v);
662 }
663 double vcorr = (resv)/lastX;
664 //
665 // 2nd iteration : create new residuals forcing them to be 0 at Xdrift=0 and maxXDrift
666 resv = res->GetBinWidth(ib0)*vcorr;
667 for (int i=ib0;i<=lastB;i++) {
668 double dx = res->GetBinWidth(i);
669 double v = vd->GetBinContent(i);
670 resv += dx*(1.-1./v - vcorr);
671 res->SetBinContent(i, resv);
672 }
673 //
674 return res;
675}
676
677//__________________________________________________________________________________
678void CheckResultDX()
679{
680 // check mean residuals before and after correction
681 const double kFOffs = 0.05; // skip edges
682 static TF1* fitP1 = new TF1("fitP1","pol1",-5000,40000);
683 static TF1* fitP0 = new TF1("fitP0","pol0",-5000,40000);
684 //
685 if (currHDXvsX->GetEntries()<minBEntryX*currHDXvsX->GetNbinsX()) {
686 if (currHCorrMapX) currHCorrMapX->Reset();
687 return;
688 }
689 //
690 // vdrift correction
691 int b0 = currHDXvsXCorr->FindBin(1),b1 = currHDXvsXCorr->FindBin(sddSeg.Dx()-1);
692 int nb = 0;
693 for (int ib=b0;ib<b1;ib++) if (currHTDvsX->GetBinEntries(ib)>=minBEntryX) nb++;
694 currHDXvsX->Fit(fitP0,"q","");
695 double offsRaw = fitP0->GetParameter(0);
696 currHDXvsXCorr->Fit(fitP0,"q0","");
697 double offsAMap = fitP0->GetParameter(0);
698
699 //
700 currGrDxx = new TGraphErrors(nb); // residual vs Xdrift
701 currGrTDx = new TGraphErrors(nb); // residual vs tdrift
702 currGrTXCorr = new TGraphErrors(nb); // Xdrift vs tdrift
703 double tmin = 1e6, tmax = -1e6, xmin = 1e6, xmax = -1e6;
704 int ip = 0;
705 for (int i=b0;i<b1;i++) {
706 if (currHTDvsX->GetBinEntries(i)<minBEntryX) continue;
707 double t = currHTDvsX->GetBinContent(i);
708 double x = currHDXvsXCorr->GetBinCenter(i);
709 if (tmin>t) tmin = t;
710 if (tmax<t) tmax = t;
711 if (xmin>x) xmin = x;
712 if (xmax<x) xmax = x;
713 currGrDxx->SetPoint(ip,x,currHDXvsXCorr->GetBinContent(i));
714 currGrDxx->SetPointError(ip,currHDXvsXCorr->GetBinWidth(i),currHDXvsXCorr->GetBinError(i));
715 //
716 currGrTDx->SetPoint(ip,t, currHDXvsXCorr->GetBinContent(i));
717 currGrTDx->SetPointError(ip, currHTDvsX->GetBinError(i), currHDXvsXCorr->GetBinError(i));
718 //
719 currGrTXCorr->SetPoint(ip,t, currHTDvsX->GetBinCenter(i));
720 currGrTXCorr->SetPointError(ip,currHTDvsX->GetBinError(i), currHTDvsX->GetBinWidth(i));
721 //
722 ip++;
723 }
724 double del = tmax-tmin;
725 tmin += kFOffs*del;
726 tmax -= kFOffs*del;
727 del = xmax - xmin;
728 xmin += kFOffs*del;
729 xmax -= kFOffs*del;
730 //
731 currGrDxx->Fit(fitP1,"q","",xmin, xmax);
732 double offs = fitP1->GetParameter(0); // offset of correction line at Xdrift=0
733 //
734 // printf("Fitting VD correction in the range %.1f:%.1f\n",tmin,tmax);
735 currGrTDx->Fit(fitP1,"q","",tmin,tmax);
736 double vcor = fitP1->GetParameter(1);
737 //
738 currGrTXCorr->Fit(fitP1,"q","",tmin,tmax);
739 double vav = fitP1->GetParameter(1);
740 //
741 // store results
742 if (!resOffsDXraw[currSDDSide]) {
743 resOffsDXraw[currSDDSide] = new TH1F(Form("OffsRaw%d",currSDDSide),Form("DX Raw Offset, side %d",currSDDSide),kNSDD,kSDDMin-0.5,kSDDMax+0.5);
744 resOffsDXraw[currSDDSide]->SetMarkerColor(2+currSDDSide);
745 resOffsDXraw[currSDDSide]->SetMarkerStyle(20+4*currSDDSide);
746 }
747 //
748 if (!resOffsDXAMap[currSDDSide]) {
749 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);
750 resOffsDXAMap[currSDDSide]->SetMarkerColor(3+currSDDSide);
751 resOffsDXAMap[currSDDSide]->SetMarkerStyle(20+4*currSDDSide);
752 }
753 //
754 if (!resOffsDX[currSDDSide]) {
755 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);
756 resOffsDX[currSDDSide]->SetMarkerColor(2+currSDDSide);
757 resOffsDX[currSDDSide]->SetMarkerStyle(20+4*currSDDSide);
758 }
759 //
760 if (!resVDCorr[currSDDSide]) {
761 resVDCorr[currSDDSide] = new TH1F(Form("VDcorr%d",currSDDSide),Form("VDrift correction, side %d",currSDDSide),kNSDD,kSDDMin-0.5,kSDDMax+0.5);
762 resVDCorr[currSDDSide]->SetMarkerColor(2+currSDDSide);
763 resVDCorr[currSDDSide]->SetMarkerStyle(20+4*currSDDSide);
764 }
765 //
766 if (!resVDMean[currSDDSide]) {
767 resVDMean[currSDDSide] = new TH1F(Form("VDmean%d",currSDDSide),Form("VDrift mean, side %d",currSDDSide),kNSDD,kSDDMin-0.5,kSDDMax+0.5);
768 resVDMean[currSDDSide]->SetMarkerColor(2+currSDDSide);
769 resVDMean[currSDDSide]->SetMarkerStyle(20+4*currSDDSide);
770 }
771 //
772 resOffsDXraw[currSDDSide]->SetBinContent(currMod+1, offsRaw);
773 resOffsDXAMap[currSDDSide]->SetBinContent(currMod+1, offsAMap);
774 resOffsDX[currSDDSide]->SetBinContent(currMod+1, offs);
775 resVDCorr[currSDDSide]->SetBinContent(currMod+1, vcor);
776 resVDMean[currSDDSide]->SetBinContent(currMod+1, vav);
777 //
778}
779
780//__________________________________________________________________________
781void CalcDXCorrections()
782{
783 // estimate time0 and alignment correction for the whole module
784 if (!resT0Corr) {
785 resT0Corr = new TH1F("T0Corr","T0 Correction",kNSDD,kSDDMin-0.5,kSDDMax+0.5);
786 resT0Corr->SetMarkerColor(2);
787 resT0Corr->SetMarkerStyle(20);
788 }
789 //
790 if (!resXLocCorr) {
791 resXLocCorr = new TH1F("XLocCorr","XLoc Correction",kNSDD,kSDDMin-0.5,kSDDMax+0.5);
792 resXLocCorr->SetMarkerColor(2);
793 resXLocCorr->SetMarkerStyle(20);
794 }
795 //
796 if (!resVDMean[0] || !resVDMean[1]) return;
797 if (!resOffsDX[0] || !resOffsDX[1]) return;
798 double vL = resVDMean[0]->GetBinContent(currMod+1); // true mean VL
799 double vR = resVDMean[1]->GetBinContent(currMod+1); // true mean VR
800 double offsL = resOffsDX[0]->GetBinContent(currMod+1);
801 double offsR = resOffsDX[1]->GetBinContent(currMod+1);
802 //
803 double vsum=0,t0Corr=0,xlCorr=0;
804 if (vL>1 && vR>1) { // both sides available
805 vsum = vL + vR;
806 t0Corr = -(offsL+offsR)/vsum;
807 xlCorr = -(offsL*vR - offsR*vL)/vsum;
808 }
809 /*
810 else if (vL>1) t0Corr = -offsL/vL; // only one side is available
811 else if (vR>1) t0Corr = -offsR/vR;
812 */
813 else if (vL>1) xlCorr = -offsL; // only one side is available
814 else if (vR>1) xlCorr = offsR;
815 //
816 // printf("SDD%d VL:%f VR:%f offsL:%+f offsR:%+f dT:%+f dX:%+f\n",currSDDId, vL,vR, offsL,offsR, t0Corr,xlCorr);
817 resT0Corr->SetBinContent(currMod+1, t0Corr); // T0 correction
818 resXLocCorr->SetBinContent(currMod+1, xlCorr); // X alignment correction
819 //
820 double addMap[2]={0,0};
821 Bool_t redoMaps = kFALSE;
822 //
823 if (forceT0CorrTo0) { // T0 correction was forced to be 0, attribute this to map
824 addMap[0] -= vL*t0Corr;
825 addMap[1] -= vR*t0Corr;
826 redoMaps = kTRUE;
827 }
828 if (forceRPhiCorrTo0) { // alignment correction was forced to be 0, attribute this to map
829 addMap[0] -= xlCorr;
830 addMap[1] -= -xlCorr;
831 redoMaps = kTRUE;
832 }
833 //
834 if (redoMaps) {
835 for (int ix=0;ix<2;ix++) {
836 TH1* map = (TH1*)procHistos.At( GetStoreID(kSCorrMapX, currSDDId, ix) );
837 TH1* mapc = (TH1*)procHistos.At( GetStoreID(kSXvsXCorr, currSDDId, ix) );
838 if (!map || !mapc) continue;
839 int ib0 = map->FindBin(1);
840 int ib1 = map->FindBin(sddSeg.Dx()-1);
841 for (int ib=ib0+1;ib<ib1;ib++) {
842 map->AddBinContent(ib, addMap[ix]);
843 mapc->AddBinContent(ib, -addMap[ix]);
844 }
845 }
846 }
847
848 //
849}
850
851//______________________________________________________________
852Int_t GetStoreID(int type, int imd,int side)
853{
854 // entry of the histo/graph of type in the procHistos array
855 //
856 if (imd<0) imd = currSDDId;
857 if (side<0) side = currSDDSide;
858 if (type<0||type>=kNStore || imd<kSDDMin || imd>kSDDMax || side<0 || side>1) {
859 printf("Wrong object requested: type: %d, Mod:%d/%d\n",type,imd,side);
860 exit(1);
861 }
862 return (2*(imd-kSDDMin)+side)*kNStore + type;
863}
864
865//______________________________________________________________
866void CleanPrev()
867{
868 // clean "current" objects from last event
869 currHDXvsX = 0;
870 currHDXvsZ = 0;
871 currHTDvsX = 0;
872 currHTDvsZ = 0;
873 currHDXvsXclean = 0;
874 currHVDvxX = 0;
875 currHCorrMapX = 0;
876 currHDXvsXCorr = 0;
877 currHDVvsZ = 0;
878 currGrDxx = 0;
879 currGrTDx = 0;
880 currGrTXCorr = 0;
881 //
882}
883
884//______________________________________________________________
885void StoreCurrent()
886{
887 // store "current" objects in procHistos
888 if (currHDXvsX) procHistos.AddAtAndExpand(currHDXvsX, GetStoreID(kSXvsX));
889 if (currHDXvsZ) procHistos.AddAtAndExpand(currHDXvsZ, GetStoreID(kSXvsZ));
890 if (currHTDvsX) procHistos.AddAtAndExpand(currHTDvsX, GetStoreID(kSTDvsX));
891 if (currHTDvsZ) procHistos.AddAtAndExpand(currHTDvsZ, GetStoreID(kSTDvsZ));
892 if (currHDXvsXclean) procHistos.AddAtAndExpand(currHDXvsXclean, GetStoreID(kSDXvsXclean));
893 if (currHVDvxX) procHistos.AddAtAndExpand(currHVDvxX, GetStoreID(kSVDvxX));
894 if (currHCorrMapX) procHistos.AddAtAndExpand(currHCorrMapX, GetStoreID(kSCorrMapX));
895 if (currHDXvsXCorr) procHistos.AddAtAndExpand(currHDXvsXCorr, GetStoreID(kSXvsXCorr));
896 if (currHDXvsZCorrMT0) procHistos.AddAtAndExpand(currHDXvsZCorrMT0, GetStoreID(kSXvsZCorr));
897 if (currHDVvsZ) procHistos.AddAtAndExpand(currHDVvsZ, GetStoreID(kSDVvsZ));
898 if (currHDVvsZOrig) procHistos.AddAtAndExpand(currHDVvsZOrig, GetStoreID(kSDVvsZOrig));
899 if (currGrDxx) procHistos.AddAtAndExpand(currGrDxx, GetStoreID(kSGrDxx));
900 if (currGrTDx) procHistos.AddAtAndExpand(currGrTDx, GetStoreID(kSGrTDx));
901 if (currGrTXCorr) procHistos.AddAtAndExpand(currGrTXCorr, GetStoreID(kSGrTXCorr));
902 //
903}
904
905//_________________________________________________________________________________
906TObjArray* CreateCorrMaps()
907{
908 // create correction maps for all modules
909 printf("Creating correction maps (update %s)\n",pathSDDCorrMapOld.Data());
910 TObjArray *dest = new TObjArray(2*kNSDD);
911 TObjArray* update = 0;
912 if (!pathSDDCorrMapOld.IsNull() && !LoadSDDCorrMap(pathSDDCorrMapOld,update)) {
913 printf("The update of correction map was requested but the source %s is not found\n",pathSDDCorrMapOld.Data());
914 exit(1);
915 }
916 //
917 dest->Clear();
918 AliITSCorrMap1DSDD *updMap = 0;
919 for (int imd=kSDDMin;imd<=kSDDMax;imd++) {
920 for (int side=0;side<2;side++) {
921 TH1* mph = (TH1*)procHistos.At( GetStoreID(kSCorrMapX,imd,side) );
922 //if (!mph) printf("Correction map is missing for module %d/%d\n",imd,side);
923 if (update) updMap = (AliITSCorrMap1DSDD*)update->At(2*(imd-kSDDMin) + side);
924 AliITSCorrMap1DSDD* mp = CreateCorrMap(mph,imd,side, updMap);
925 dest->AddAtAndExpand(mp, 2*(imd-kSDDMin) + side);
926 }
927 }
928 //
929 return dest;
930}
931
932//_________________________________________________________________________________
933AliITSCorrMap1DSDD* CreateCorrMap(TH1* mapHisto, int imd, int side, AliITSCorrMap1DSDD* updateMap)
934{
935 // create or update correction map from histo
936 int nbCorr = 1, nbOld = 0;
937 int b0=0,b1=0;
938 if (mapHisto) {
939 b0 = mapHisto->FindBin(1);
940 b1 = mapHisto->FindBin(sddSeg.Dx()-1);
941 }
942 nbCorr = b1-b0+1;
943 AliITSCorrMap1DSDD* mpCorr = 0;
944 //
977bf0d4 945 // check if the updateMap is meaningful
946 if (updateMap && updateMap->GetNBinsDrift()>2 && nbCorr>1) {
ba2089c4 947 if (mapHisto) {
948 TSpline3 spl(mapHisto);
949 nbOld = updateMap->GetNBinsDrift();
950 double dx = sddSeg.Dx()/nbOld;
951 for (int ip=0;ip<nbOld;ip++) {
952 double x = dx*(0.5+ip);
953 updateMap->SetCellContent(0,ip,updateMap->GetCellContent(0,ip)-spl.Eval(x));
954 }
955 }
956 mpCorr = updateMap;
957 }
958 else {
959 mpCorr = new AliITSCorrMap1DSDD(Form("DriftTimeMap_%d_%d",imd,side),nbCorr);
960 if (side==0) mpCorr->SetInversionBit(); // !!! left side should return correction*-1
961 if (mapHisto) for (int ib=b0;ib<=b1;ib++) mpCorr->SetCellContent(0,ib-b0,-mapHisto->GetBinContent(ib));
962 }
963 //
964 return mpCorr;
965}
966
967//_________________________________________________________________________________
968TObjArray* UpdateSDDVDrift()
969{
970 // retrieve SDD VDrift object and update it
971 if (!vdarrayOld && !LoadSDDVDrift(pathSDDVDriftOld,vdarrayOld)) return 0;
972 TObjArray *vdarrayUp = new TObjArray(2*kNSDD);
973 //
974 for (int imd=kSDDMin;imd<=kSDDMax;imd++) {
975 for (int side=0;side<2;side++) {
976 int iad = 2*(imd-kSDDMin)+side;
977 AliITSDriftSpeedArraySDD* drarr = (AliITSDriftSpeedArraySDD*) vdarrayOld->At( iad );
978 AliITSDriftSpeedArraySDD* drarrUp = new AliITSDriftSpeedArraySDD();
979 AliITSDriftSpeedSDD* vOr = drarr->GetDriftSpeedObject(0);
980 AliITSDriftSpeedSDD* vUp = new AliITSDriftSpeedSDD(*vOr);
981 drarrUp->AddDriftSpeed(vUp);
982 vdarrayUp->AddAt(drarrUp, iad);
983 UpdateSDDVDrift(drarrUp, imd, side);
984 }
985 }
986 //
987 sddVDriftUpdOK = kTRUE;
988 return vdarrayUp;
989}
990
991
992//_________________________________________________________________________________
993void UpdateSDDVDrift(AliITSDriftSpeedArraySDD* vdArr, int imd, int side)
994{
995 // update vdrift vs anode in the object
996 AliITSDriftSpeedSDD* ds;
997 if (!vdArr || !(ds=vdArr->GetDriftSpeedObject(0))) {printf("No VDrift object for module %d/%d\n",imd,side); exit(1);}
998 TH1* vdh = (TH1*)procHistos.At( GetStoreID(kSDVvsZ,imd,side) );
999 if (!vdh) {
1000 //printf("VDrift vs Z correction is not processed for module %d/%d\n",imd,side);
1001 return;
1002 }
1003 TF1* fp = vdh->GetFunction("fitVvsZ");
1004 if (!fp) {printf("VDrift vs Z correction fit is missing SDD%d/%d\n",imd,side); return;}
1005 //
1006 int ord = (int)fp->GetParameter(0); // 1st param is the order of poly
1007 int ordOld = ds->GetDegreeofPoly();
1008 if (ord>ordOld) ds->SetDegreeofPoly(ord);
1009 for (int ip=0;ip<ord+1;ip++) { // don't store offset (par[1])
1010 double par = ds->GetDriftSpeedParameter(ip) - fp->GetParameter(ip+1);
1011 ds->SetDriftSpeedParameter(ip, par);
1012 }
1013 //
1014}
1015
1016//_________________________________________________________________________________
1017AliITSresponseSDD* UpdateSDDResponse(Bool_t t0, Bool_t vdrift)
1018{
1019 // retrieve RespSDD object and update it
1020 AliITSresponseSDD* resp = 0;
1021 if (!LoadSDDResponse(pathSDDRespOld, resp)) return 0;
1022 UpdateSDDResponse(resp, t0, vdrift);
1023 sddRespUpdT0OK = t0;
1024 sddRespUpdVDOK = vdrift;
1025 //
1026 return resp;
1027}
1028
1029//_________________________________________________________________________________
1030void UpdateSDDResponse(AliITSresponseSDD *resp, Bool_t t0, Bool_t vdrift)
1031{
1032 // update the map with extracted values
1033 printf("Updating RespSDD object: T0:%s VDrift:%s\n",t0?"ON":"OFF",vdrift?"ON":"OFF");
1034 //
1035 if (t0 && !resT0Corr)
1036 {printf("T0 update is requested but corrections were not processed"); exit(1);}
1037 if (vdrift && !(resVDCorr[0] && resVDCorr[1]))
1038 {printf("VDrift update is requested but corrections were not processed"); exit(1);}
1039 //
1040 for (int imd=kSDDMin;imd<=kSDDMax;imd++) {
1041 if (t0 && !forceT0CorrTo0) resp->SetModuleTimeZero(imd, resp->GetTimeZero(imd) - resT0Corr->GetBinContent(imd-kSDDMin+1));
1042 if (vdrift) {
1043 for (int ix=0;ix<2;ix++) {
1044 double vdZ = sddVDriftUpdOK&&resVDCorrZ[ix] ? resVDCorrZ[ix]->GetBinContent(imd-kSDDMin+1) : 0; // contribution from DXvsZ correction
1045 double vdX = resVDCorr[ix]->GetBinContent(imd-kSDDMin+1); // contribution from DXvsX correction
1046 resp->SetDeltaVDrift(imd, resp->GetDeltaVDrift(imd,ix) - (vdX-vdZ), ix);
1047 }
1048 }
1049 }
1050 //
1051}
1052
1053//___________________________________________________________________
1054double GetVOld(double z)
1055{
1056 // return VDrift assumed in reconstruction
1057 if (!vdarrayOld && !LoadSDDVDrift(pathSDDVDriftOld,vdarrayOld)) return 0;
1058 AliITSDriftSpeedArraySDD* drarr = (AliITSDriftSpeedArraySDD*) vdarrayOld->At( 2*currMod + currSDDSide);
1059 float anode = sddSeg.GetAnodeFromLocal( currSDDSide==0 ? 1.:-1. ,z*1e-4);
1060 double v = drarr->GetDriftSpeed(0, anode);
1061 return v;
1062}
1063
1064//___________________________________________________________________
1065double ftVdZ(double *x, double *par)
1066{
1067 // function to fit the vdrift dependence on Z
1068 //
1069 // convert Z to anode
1070 double z = x[0];
1071 double ian = (z/sddSeg.Dz() + 0.5);
1072 if (ian<0) ian = 0.;
1073 else if (ian>1) ian = 1.;
1074 ian *= sddSeg.NpzHalf();
1075 //
1076 int ord = int(par[0]);
1077 double v = par[ord+1];
1078 for (int i=ord;i--;) v = par[i+1]+ian*v;
1079 return v;
1080}
1081
1082//________________________________________________________________________________________________________
1083Bool_t LoadSDDVDrift(TString& path, TObjArray *&arr)
1084{
1085 // load VDrift object
1086 if (path.IsNull()) return kFALSE;
1087 printf("Loading SDD VDrift from %s\n",path.Data());
1088 //
1089 AliCDBEntry *entry = 0;
1090 delete arr;
1091 arr = 0;
1092 while(1) {
1093 if (path.BeginsWith("path: ")) { // must load from OCDB
1094 entry = GetCDBEntry(path.Data());
1095 if (!entry) break;
1096 arr = (TObjArray*) entry->GetObject();
1097 entry->SetObject(NULL);
1098 entry->SetOwner(kTRUE);
1099 break;
1100 }
1101 //
1102 if (gSystem->AccessPathName(path.Data())) break;
1103 TFile* precf = TFile::Open(path.Data());
1104 if (precf->FindKey("TObjArray")) arr = (TObjArray*)precf->Get("TObjArray");
1105 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
1106 arr = (TObjArray*) entry->GetObject();
1107 if (arr && arr->InheritsFrom(TObjArray::Class())) entry->SetObject(NULL);
1108 else arr = 0;
1109 entry->SetObject(NULL);
1110 entry->SetOwner(kTRUE);
1111 delete entry;
1112 }
1113 //
1114 precf->Close();
1115 delete precf;
1116 break;
1117 }
1118 //
1119 if (!arr) {printf("Failed to load SDD vdrift from %s\n",path.Data()); return kFALSE;}
1120 arr->SetOwner(kTRUE);
1121 return kTRUE;
1122}
1123
1124//________________________________________________________________________________________________________
1125Bool_t LoadSDDResponse(TString& path, AliITSresponseSDD *&resp)
1126{
1127 // load SDD response
1128 if (path.IsNull()) return kFALSE;
1129 printf("Loading SDD response from %s\n",path.Data());
1130 //
1131 AliCDBEntry *entry = 0;
1132 delete resp;
1133 resp = 0;
1134 //
1135 while(1) {
1136 if (path.BeginsWith("path: ")) { // must load from OCDB
1137 entry = GetCDBEntry(path.Data());
1138 if (!entry) break;
1139 resp = (AliITSresponseSDD*) entry->GetObject();
1140 entry->SetObject(NULL);
1141 entry->SetOwner(kTRUE);
1142 break;
1143 }
1144 //
1145 if (gSystem->AccessPathName(path.Data())) break;
1146 TFile* precf = TFile::Open(path.Data());
1147 if (precf->FindKey("AliITSresponseSDD")) resp = (AliITSresponseSDD*)precf->Get("AliITSresponseSDD");
1148 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
1149 resp = (AliITSresponseSDD*) entry->GetObject();
1150 if (resp && resp->InheritsFrom(AliITSresponseSDD::Class())) entry->SetObject(NULL);
1151 else resp = 0;
1152 entry->SetObject(NULL);
1153 entry->SetOwner(kTRUE);
1154 delete entry;
1155 }
1156 //
1157 precf->Close();
1158 delete precf;
1159 break;
1160 }
1161 //
1162 if (!resp) {printf("Error: Failed to load SDD response from %s\n",path.Data()); return kFALSE;}
1163 return kTRUE;
1164}
1165
1166//________________________________________________________________________________________________________
1167Bool_t LoadSDDCorrMap(TString& path, TObjArray *&maps)
1168{
1169 // Load SDD correction map
1170 //
1171 if (path.IsNull()) return kFALSE;
1172 printf("Loading SDD Correction Maps from %s\n",path.Data());
1173 //
1174 AliCDBEntry *entry = 0;
1175 delete maps;
1176 maps = 0;
1177 while(1) {
1178 if (path.BeginsWith("path: ")) { // must load from OCDB
1179 entry = GetCDBEntry(path.Data());
1180 if (!entry) break;
1181 maps = (TObjArray*) entry->GetObject();
1182 entry->SetObject(NULL);
1183 entry->SetOwner(kTRUE);
1184 break;
1185 }
1186 //
1187 if (gSystem->AccessPathName(path.Data())) break;
1188 TFile* precf = TFile::Open(path.Data());
1189 if (precf->FindKey("TObjArray")) maps = (TObjArray*)precf->Get("TObjArray");
1190 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
1191 maps = (TObjArray*) entry->GetObject();
1192 if (maps && maps->InheritsFrom(TObjArray::Class())) entry->SetObject(NULL);
1193 else maps = 0;
1194 entry->SetObject(NULL);
1195 entry->SetOwner(kTRUE);
1196 delete entry;
1197 }
1198 //
1199 precf->Close();
1200 delete precf;
1201 break;
1202 }
1203 //
1204 if (!maps) {printf("Failed to load SDD Correction Map from %s\n",path.Data()); return kFALSE;}
1205
1206 return kTRUE;
1207}
1208
1209//_______________________________________________________________________________________
1210AliCDBEntry* GetCDBEntry(const char* path)
1211{
1212 // return object from the OCDB
1213 AliCDBEntry *entry = 0;
1214 printf("Loading object %s\n",path);
1215 AliCDBManager* man = AliCDBManager::Instance();
1216 AliCDBId* cdbId = AliCDBId::MakeFromString(path);
1217 if (!cdbId) {
1218 printf("Failed to create cdbId\n");
1219 return 0;
1220 }
1221 //
1222 AliCDBStorage* stor = man->GetDefaultStorage();
1223 if (!stor && !man->GetRaw()) man->SetDefaultStorage("raw://");
1224 if (man->GetRaw()) man->SetRun(cdbId->GetFirstRun());
1225 if (stor) {
1226 TString tp = stor->GetType();
1227 if (tp.Contains("alien",TString::kIgnoreCase) && !gGrid) TGrid::Connect("alien:");
1228 }
1229 entry = man->Get(cdbId->GetPath(),cdbId->GetFirstRun(),cdbId->GetVersion(),cdbId->GetSubVersion());
1230 // entry = man->Get( *cdbId );
1231 man->ClearCache();
1232 //
1233 delete cdbId;
1234 return entry;
1235 //
1236}
1237//
1238
1239//_______________________________________________________________________________________
1240Bool_t PlotHisto(TH1* h, Option_t* opt, int mrkStyle,int mrkCol, double mrkSize)
1241{
1242 const double kOffsH = 0.15;
1243 if (!h) return kFALSE;
1244 TString opts = opt; opts.ToLower();
1245 if (opts.Contains("p")) {
1246 h->SetMarkerStyle(mrkStyle);
1247 h->SetMarkerColor(mrkCol);
1248 h->SetMarkerSize(mrkSize);
1249 }
1250 h->SetLineColor(mrkCol);
1251 h->Draw(opt);
1252 //
1253 h->SetMinimum(); h->SetMaximum();
1254 double hmn=h->GetMinimum(),hmx=h->GetMaximum(); // new histo min/max
1255 //
1256 TH1* hbase = GetPadBaseHisto((TPad*)gPad); if (!hbase) return 0;
1257 double smn = hbase->GetMinimum(),smx = hbase->GetMaximum(); // base set min/max?
1258 hbase->SetMinimum(); hbase->SetMaximum();
1259 double omn = hbase->GetMinimum(),omx = hbase->GetMaximum(); // base real min max
1260 if (smn<omn && smx>omx) { // min/max for bas histo was set by hand: extract original min/max
1261 omx = (smn*kOffsH+smx*(1+kOffsH))/(1+2*kOffsH);
1262 omn = (smn-kOffsH*omx)/(1+kOffsH);
1263 }
1264 if (hmn<omn) omn = hmn;
1265 if (hmx>omx) omx = hmx;
1266 double del = omx-omn;
1267 hbase->SetMinimum( omn - kOffsH*del );
1268 hbase->SetMaximum( omx + kOffsH*del );
1269 gPad->Update();
1270 return kTRUE;
1271}
1272
1273//_______________________________________________________________________________________
1274void GetHistoMinMaxInRange(TH1* h, double &mn,double &mx)
1275{
1276 // compute min/max of histo in the user range
1277 mn = 1e50;
1278 mx =-1e50;
1279 int b0 = h->GetXaxis()->GetFirst(), b1 = h->GetXaxis()->GetLast();
1280 for (int i=b0;i<=b1;i++) {
1281 double e = h->GetBinError(i);
1282 if (TMath::Abs(e)<1e-9) continue;
1283 double v = h->GetBinContent(i);
1284 if (mn>v-e) mn = v-e;
1285 if (mx<v+e) mx = v+e;
1286 }
1287}
1288
1289//_______________________________________________________________________________________
1290void PlotReport(const char* psname)
1291{
1292 // report results
1293 sddCanv = new TCanvas("sddCanv","sddCanv",700,900);
1294 //
1295 gStyle->SetOptStat(0);
1296 gStyle->SetOptTitle(0);
1297 //
1298 TString psnm1 = psname;
1299 if (psnm1.IsNull()) psnm1 = "sddQAreport.ps";
1300 TString psnm0 = psnm1 + "[";
1301 TString psnm2 = psnm1 + "]";
1302 sddCanv->Print(psnm0.Data());
1303 //
1304 // mean corrections per module/side
1305 sddCanv->Clear();
1306 sddCanv->Divide(2,3);
1307 int cntPad = 0;
1308 //
1309 for (int ix=0;ix<2;ix++) { // mean residuals before/after correction
1310 sddCanv->cd(++cntPad);
1311 PlotHisto(resOffsDXraw[ix],"p" ,7,kBlack,1);
1312 PlotHisto(resOffsDX[ix] ,"p same",7,kRed ,1);
1313 AddPadLabel(Form("<#DeltaX> %s : Raw",ix?"Right":"Left"), 0.1,0.93,kBlack,0.05);
1314 AddPadLabel("After Map", 0.5,0.93,kRed,0.05);
1315 }
1316 //
1317 for (int ix=0;ix<2;ix++) { // mean residuals before/after correction
1318 sddCanv->cd(++cntPad);
1319 PlotHisto(resVDCorr[ix] ,"p" ,7,kBlack,1);
1320 PlotHisto(resVDCorrZ[ix],"p same",7,kRed ,1);
1321 AddPadLabel(Form("<#DeltaV> %s : from #DeltaX vs X",ix?"Right":"Left"), 0.1,0.93,kBlack,0.05);
1322 AddPadLabel("from #DeltaX vs Z", 0.6,0.93,kRed,0.05);
1323 }
1324 //
1325 sddCanv->cd(++cntPad);
1326 PlotHisto(resT0Corr,"p",7,kBlue,1);
1327 AddPadLabel("T0 Correction", 0.3,0.93,kRed,0.05);
1328 if (forceT0CorrTo0) AddPadLabel("Forced to 0 by transferring to maps", 0.15,0.88,kRed,0.05);
1329 //
1330 sddCanv->cd(++cntPad);
1331 PlotHisto(resXLocCorr,"p",7,kBlack,1);
1332 AddPadLabel("#DeltaX Correction", 0.3,0.93,kRed,0.05);
1333 if (forceRPhiCorrTo0) AddPadLabel("Forced to 0 by transferring to maps", 0.15,0.88,kRed,0.05);
1334 //
1335 sddCanv->cd();
1336 sddCanv->Print(psnm1.Data());
1337 //
1338 //-------------------------------------------------------------------
1339 TH1* hsdd = 0;
1340 //
1341 cntPad = 999;
1342 int nModPerPage = 3;
1343 int nRowPerMod = 2;
1344 Bool_t saved = kFALSE;
1345 int ib0=1,ib1=999;
1346 //
1347 for (int imd=kSDDMin;imd<=kSDDMax;imd++) {
1348 if (cntPad>=2*nModPerPage*nRowPerMod) {
1349 sddCanv->cd();
1350 if (imd!=kSDDMin) sddCanv->Print(psnm1.Data());
1351 sddCanv->Clear();
1352 sddCanv->Divide(2,nModPerPage*nRowPerMod);
1353 cntPad = 0;
1354 saved = kTRUE;
1355 }
1356 for (int ix=0;ix<2;ix++) {
1357 sddCanv->cd(++cntPad);
1358 hsdd = (TH1*)procHistos.At(GetStoreID(kSXvsX,imd,ix)); // raw residuals
1359 if (hsdd) {
1360 ib0 = hsdd->FindBin(1);
1361 ib1 = hsdd->FindBin(sddSeg.Dx()-1);
1362 hsdd->GetXaxis()->SetRange(ib0,ib1);
1363 }
1364 PlotHisto(hsdd," ",7,kBlack,1);
1365 hsdd = (TH1*)procHistos.At(GetStoreID(kSCorrMapX,imd,ix)); // map
1366 if (hsdd) hsdd->GetXaxis()->SetRange(ib0,ib1);
1367 PlotHisto(hsdd,"same",7,kRed,1);
1368 hsdd = (TH1*)procHistos.At(GetStoreID(kSXvsXCorr,imd,ix)); // map
1369 if (hsdd) hsdd->GetXaxis()->SetRange(ib0,ib1);
1370 PlotHisto(hsdd,"same",7,kBlue,1);
1371 //
1372 AddPadLabel(Form("<#DeltaX> %d %s: Raw",imd,ix?"Right":"Left"), 0.1,0.93,kBlack,0.07);
1373 AddPadLabel("Map", 0.4,0.93,kRed,0.07);
1374 AddPadLabel("Corr.Map", 0.5,0.93,kBlue,0.07);
1375 //
1376 AddPadLabel(Form("#DeltaV:%+.4f | #DeltaT0:%+5.0f | #DeltaX:%+4.0f",
1377 resVDCorr[ix] ? resVDCorr[ix]->GetBinContent(imd-kSDDMin+1):0,
1378 resT0Corr ? resT0Corr->GetBinContent(imd-kSDDMin+1):0,
1379 resXLocCorr ? resXLocCorr->GetBinContent(imd-kSDDMin+1):0),
1380 0.5, 0.15, kRed, 0.07);
1381 //
1382 saved = kFALSE;
1383 }
1384 //
1385 for (int ix=0;ix<2;ix++) {
1386 sddCanv->cd(++cntPad);
1387 hsdd = (TH1*)procHistos.At(GetStoreID(kSDVvsZ,imd,ix)); // correction Vd vs Z
1388 PlotHisto(hsdd," ",7,kBlack,1);
1389 hsdd = (TH1*)procHistos.At(GetStoreID(kSDVvsZOrig,imd,ix)); // correction Vd vs Z
1390 PlotHisto(hsdd,"same",24,kBlue,1);
1391 AddPadLabel(Form("<#DeltaV> vs Z %d %s | Stat:%.2e",imd,ix?"Right":"Left",
1392 ((TH1*)procHistos.At(GetStoreID(kSXvsX,imd,ix)))->GetEntries()), 0.1,0.93,kBlack,0.07);
1393 //
1394 AddPadLabel(Form("<#DeltaV>:%+.4f",resVDCorrZ[ix] ? resVDCorrZ[ix]->GetBinContent(imd-kSDDMin+1):0), 0.5, 0.15, kRed, 0.07);
1395 //
1396 saved = kFALSE;
1397 }
1398 //
1399 }
1400 //
1401 sddCanv->cd();
1402 if (!saved) sddCanv->Print(psnm1.Data());
1403 sddCanv->Print(psnm2.Data());
1404}
1405
1406//__________________________________
1407TH1* GetPadBaseHisto(TPad* pad)
1408{
1409 if (!pad) pad = (TPad*)gPad;
1410 if (!pad) return 0;
1411 TList* lst = pad->GetListOfPrimitives();
1412 int size = lst->GetSize();
1413 TH1* hst=0;
1414 for (int i=0;i<size;i++) {
1415 TObject* obj = lst->At(i);
1416 if (!obj) continue;
1417 if (obj->InheritsFrom("TH1")) {hst = (TH1*)obj; break;}
1418 }
1419 return hst;
1420}
1421
1422//__________________________________
1423TLatex* AddPadLabel(const char*txt,float x,float y,int color,float size)
1424{
1425 TLatex* lt = new TLatex(x,y,txt);
1426 lt->SetNDC();
1427 lt->SetTextColor(color);
1428 lt->SetTextSize(size);
1429 lt->Draw();
1430 return lt;
1431}
1432
1433//__________________________________
1434void SetCDBObjData(int firstrun,int lastrun,const char* comment)
1435{
1436 // change range and comment of the objects to store
1437 firstRun = firstrun;
1438 lastRun = lastrun;
1439 cdbComment = comment;
1440}
1441
1442//__________________________________
1443void PrepCDBObj(TObject *obj,const char* path,int firstrun,int lastrun,const char* comment)
1444{
1445 if (firstrun<0) firstrun = 0;
1446 //
1447 AliCDBManager* man = AliCDBManager::Instance();
1448 man->UnsetDefaultStorage();
1449 man->SetDefaultStorage("local://");
1450 AliCDBMetaData* md = new AliCDBMetaData();
1451 md->SetResponsible("Ruben Shahoyan");
1452 md->SetComment(comment);
1453 AliCDBId id(path,firstrun,lastrun<=0 ? (AliCDBRunRange::Infinity()) : lastrun);
1454 //AliCDBStorage* st = man->GetStorage("local//.");
1455 man->Put(obj,id,md);
1456 //
1457}
1458
1459//__________________________________________________________
1460double edgeLow(double* x, double *par)
1461{
1462 // Low TDrift edge:
1463 // residuals assuming linear dependence of "natural" residual vs Xtrue and smeared
1464 // by the finite track resolution
1465 double x0 = par[0];
1466 double sigma = par[1];
1467 double offs = par[2];
1468 double slop = par[3];
1469 //
1470 if (sigma<1) return 0;
1471 if (x0<-sigma) return 0;
1472 //
1473 double xex = x[0];
1474 xex -= x0;
1475 //
1476 double arg = xex/sigma;
1477 arg *= arg/2;
1478 double res = arg<50 ? slop*sigma*TMath::Exp(-arg)/TMath::Sqrt(2*TMath::Pi()) : 0;
1479 double erftrm = 1.+TMath::Erf(xex/sigma/TMath::Sqrt(2));
1480 //printf("%+e %+e %+e\n",x[0],x0,erftrm);
1481 if (xex<0 && TMath::Abs(erftrm)<1e-10) res = -xex*slop;
1482 else res /= erftrm/2.;
1483 res += (offs + (slop-1.)*xex);
1484 return res;
1485 //
1486}
1487
1488//__________________________________________________________
1489double edgeHigh(double* x, double *par)
1490{
1491 // High TDrift edge
1492 // residual assuming linear dependence of "natural" residual vs Xtrue and smeared
1493 // by the finite track resolution
1494 double x0 = par[0];
1495 double sigma = par[1];
1496 double offs = par[2];
1497 double slop = par[3];
1498 double tailCorr = par[4];
1499 //
1500 if (sigma<1) return 0;
1501 if (x0<-sigma) return 0;
1502 //
1503 double xex = (x0 - x[0])*tailCorr;
1504 //
1505 double arg = xex/sigma;
1506 arg *= arg/2;
1507 double res = arg<50 ? slop*sigma*TMath::Exp(-arg)/TMath::Sqrt(2*TMath::Pi()) : 0;
1508 double erftrm = 1.+TMath::Erf(xex/sigma/TMath::Sqrt(2));
1509 if (xex<0 && TMath::Abs(erftrm)<1e-10) res = xex*slop;
1510 else res /= -erftrm/2.;
1511 res += (offs + (slop-1.)*xex);
1512 return res;
1513 //
1514}
1515
1516//_________________________________________________________________________
1517void RedoProfileErrors(TH1* profH, TProfile* prof)
1518{
1519 // cure errors of low.stat bins
1520 int nbCnt = 0, nbn = prof->GetNbinsX();
1521 double meanStat = 0, meanSpread = 0, wghStat = 0;
1522 for (int i=1;i<=nbn;i++) {
1523 double stat = prof->GetBinEntries(i);
1524 if (stat>0) {meanStat+=stat; nbCnt++;}
1525 }
1526 if (nbCnt>0) meanStat/= nbCnt; // mean occupancy
1527 //
1528 for (int i=1;i<=nbn;i++) {
1529 double stat = prof->GetBinEntries(i);
1530 if (stat<meanStat/2) continue;
1531 meanSpread += prof->GetBinError(i)*TMath::Sqrt(stat)*stat;
1532 wghStat += stat;
1533 }
1534 if (wghStat) meanSpread /= wghStat; // mean spread
1535 //
1536 for (int i=1;i<=nbn;i++) { // assign error acording to occupancy
1537 double stat = prof->GetBinEntries(i);
1538 if (stat>meanStat/2 || stat<1) continue;
1539 profH->SetBinError(i, meanSpread/TMath::Sqrt(stat));
1540 }
1541}
1542
1543//_________________________________________________________________________
1544void CureEdges(TH1* prof)
1545{
1546 // cure edges of the profile histo
1547 const double kMaxChi2 = 20.;
1548 const double kSlpDf = 0.05;
1549 static TF1* fitEdgeLow = new TF1("fitEdgeLow" ,edgeLow , -5000, sddSeg.Dx()+5000,4);
1550 static TF1* fitEdgeHigh = new TF1("fitEdgeHigh",edgeHigh, -5000, sddSeg.Dx()+5000,5);
1551 //
1552 int ndf,ib0,ib1,nbn = prof->GetNbinsX();
1553 double sigma,offs,slp,chi2,x0;
1554 //
1555 // LowT edge
1556 // find 1st non-empty bin
1557 for (ib0=1;ib0<=nbn;ib0++) if (prof->GetBinError(ib0)>1e-9) break;
1558 x0 = prof->GetBinCenter(ib0);
1559 ib1 = prof->FindBin(wDXEdge);
1560 if (ib1-ib0<minBinsEdge) ib1 = ib0+minBinsEdge;
1561 //
1562 fitEdgeLow->SetParameters(100,100,0,1);
1563 fitEdgeLow->SetParLimits(0,0, sddSeg.Dx());
1564 fitEdgeLow->SetParLimits(1,edgeSmearMinErr, edgeSmearMaxErr);
1565 fitEdgeLow->SetParLimits(3,1.-kSlpDf, 1.+kSlpDf);
1566 //
1567 prof->Fit(fitEdgeLow,"q","",prof->GetBinLowEdge(ib0)+1, prof->GetBinCenter(ib1+1)-1);
1568 chi2 = fitEdgeLow->GetChisquare();
1569 ndf = fitEdgeLow->GetNDF();
1570 if (ndf>0) chi2 /= ndf;
1571 //
1572 x0 = fitEdgeLow->GetParameter(0);
1573 sigma = fitEdgeLow->GetParameter(1);
1574 offs = fitEdgeLow->GetParameter(2);
1575 slp = fitEdgeLow->GetParameter(3);
1576 if ( chi2<kMaxChi2) {
1577 x0 += 3*sigma;
1578 ib1 = prof->FindBin(x0);
1579 for (int i=ib0;i<=ib1;i++) {
1580 if (prof->GetBinError(i)<1e-9) continue;
1581 double xb = prof->GetBinCenter(i);
1582 double polval = offs+(slp-1.)*(xb-x0);
1583 if (xb>0) prof->AddBinContent(i, polval - fitEdgeLow->Eval( xb ) );
1584 else prof->SetBinContent(i,polval);
1585 }
1586 }
1587 //
1588 // find last non-empty bin
1589 for (ib1=nbn;ib1>=1;ib1--) if (prof->GetBinError(ib1)>1e-9) break;
1590 x0 = prof->GetBinCenter(ib1);
1591 ib0 = prof->FindBin(sddSeg.Dx() - wDXEdge);
1592 if (ib1-ib0<minBinsEdge) ib0 = ib1-minBinsEdge;
1593 //
1594 fitEdgeHigh->SetParameters(prof->GetBinCenter(ib0)+wDXEdge-100,100,0,1,1.);
1595 fitEdgeHigh->SetParLimits(0,0, sddSeg.Dx()+150);
1596 fitEdgeHigh->SetParLimits(1,edgeSmearMinErr, edgeSmearMaxErr);
1597 fitEdgeHigh->SetParLimits(3,1.-kSlpDf, 1.+kSlpDf);
1598 fitEdgeHigh->SetParLimits(4,0.3, 3.);
1599 prof->Fit(fitEdgeHigh,"q+","",prof->GetBinLowEdge(ib0)+1, prof->GetBinCenter(ib1+1)+1);
1600 //
1601 chi2 = fitEdgeHigh->GetChisquare();
1602 ndf = fitEdgeHigh->GetNDF();
1603 if (ndf>0) chi2 /= ndf;
1604 //
1605 x0 = fitEdgeHigh->GetParameter(0);
1606 sigma = fitEdgeHigh->GetParameter(1);
1607 offs = fitEdgeHigh->GetParameter(2);
1608 slp = fitEdgeHigh->GetParameter(3);
1609 if ( chi2<kMaxChi2 ) {
1610 x0 -= 3*sigma;
1611 ib0 = prof->FindBin(x0);
1612 for (int i=ib0;i<=ib1;i++) {
1613 if (prof->GetBinError(i)<1e-9) continue;
1614 double xb = prof->GetBinCenter(i);
1615 double polval = offs+(slp-1.)*(xb-x0);
1616 if (xb<sddSeg.Dx()) prof->AddBinContent(i, polval - fitEdgeHigh->Eval( xb ) );
1617 else prof->SetBinContent(i,polval);
1618 }
1619 }
1620 //
1621}