ITS/macrosSDD/CalibrateSDD.C

   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 <TMinuit.h>
  16 #include "AliCDBManager.h"
  17 #include "AliCDBEntry.h"
  18 #include "AliCDBStorage.h"
  19 #include "AliITSsegmentationSDD.h"
  20 #include "AliITSCorrMapSDD.h"
  21 #include "AliITSCorrMap1DSDD.h"
  22 #include "AliITSDriftSpeedArraySDD.h"
  23 #include "AliITSDriftSpeedSDD.h"
  24 #include "AliITSresponseSDD.h"
  25 #endif
  26
  27 Bool_t verbose = kFALSE;
  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 ---------------
  51 int firstRun = -1;
  52 int lastRun  = -1;
  53 TString cdbComment = "";
  54 //--------------------------------------------------------------------
  55
  56 const Int_t kSDDMin=240,kSDDMax=499,kNSDD=kSDDMax-kSDDMin+1;
  57 enum {kSXvsX,kSXvsZ,kSTDvsX,kSTDvsZ,kSDXvsXclean,kSVDvsX,kSCorrMapX,kSXvsXCorr,kSXvsZCorr,kSDVvsZ,kSDVvsZOrig,kSGrDxx,kSGrTDx,kSGrTXCorr, kNStore};
  58 const Int_t maxVZOrd    = 4;            // max order of polinomial for VD vs Z correction (normal sensors)
  59 const Double_t chi2BadVZ = 2.9;         // treat as bad module if chi2/ndf is worse
  60 //
  61 enum {kDxX,kDxZ,kTDX,kTDZ,kNQATypes};
  62 // Prefixes for residuals histos prepared by PrepSDDResidMap.C from MP2 output
  63 // For QA histograms the name should be changed to "hpSDDResXvsXD","hpSDDResXvsZ","hpSDDDrTimevsXD","hpSDDDrTimevsZ" respectively
  64 const Char_t* kQAHistoNames[] = {    //
  65   "hpSDDResXvsXD",                   // DX vs Xdrift
  66   "hpSDDResXvsZ",                    // DX vs local Z
  67   "hpSDDDrTimevsXD",                 // TDrift-T0 vs Xdrift
  68   "hpSDDDrTimevsZ"                   // TDrift-T0 vs local Z
  69 };
  70 const Char_t* kLRNames[] = {"0","1"};  // identifiers for left/right sides in histo names
  71 /*
  72 const Char_t* kQAHistoNames[] = {    //
  73   "Xdrift_Mod_",                     // DX vs Xdrift
  74   "ZAnode_Mod_",                     // DX vs local Z
  75   "TD_vs_Xdrift_Mod_",               // TDrift-T0 vs Xdrift
  76   "TD_vs_ZAnode_Mod_"                // TDrift-T0 vs local Z
  77 };
  78 const Char_t* kLRNames[] = {"left","right"};  // identifiers for left/right sides in histo names
  79 */
  80 //
  81 const double kMaxChi2SimpleMap = 1.0; // if polinomial fit gives chi2 below this value, make map out of it
  82 double minBEntryX = 10;        // min entries in the bin to account (vs Xdrift)
  83 double minBEntryZ = 20;        // min entries in the bin to account (vs Z)
  84 double skipDXEdge = 1000.;     // microns to skip from each side of XDrift when smoothing
  85 double edgeSmearMinErr = 60.;        // minimal error for extrapolation
  86 double edgeSmearMaxErr = 600.;       // maximal error for extrapolation
  87 double wDXEdge    = 5000.;     // width in microns for XDrift edges to average
  88 double wVDEdge    = 2500;      // width in microns for VDrift averaging region on the edges
  89 double threshClean= 0.3;       // clean edge bins with occupancy smaller than in the neighbourhood
  90 int    nVDEdge    = 20;        // min width in n bins for VDrift averaging region on the edges
  91 int    smoothLevel= 5;         // smoothing level
  92 int    minBinsEdge = 5;        // min number of X>0 good bins to fit the edge
  93 Bool_t useSplieForR2V = kFALSE;      // if true, extract Vd profile using Derivateve of spline (may be inadequate for sharp turns)
  94 //
  95
  96 Bool_t userLeftRightFromTASK = kTRUE;   // VERY IMPORTANT: histos produced by AliAnalysisTaskITSAlignQA have inverted left/right side definitions
  97 Bool_t forceT0CorrTo0   = kFALSE;//kTRUE;
  98 Bool_t forceRPhiCorrTo0 = kTRUE;
  99 Bool_t userDummyCorrMap = kFALSE;//kTRUE;
 100 Bool_t userDummyt0Corr  = kFALSE;//kTRUE;
 101 Bool_t userDummyxlCorr  = kFALSE;//kTRUE;
 102 Int_t  userRebinX = 1;
 103 Int_t  userRebinZ = 2;
 104 //
 105 //
 106 AliITSsegmentationSDD sddSeg;
 107 //
 108 //----------------- working variables --------------------------------
 109 Int_t currSDDId=-1,currSDDSide=-1,currMod=-1;
 110 TProfile*  currHDXvsX = 0;                // raw DX vs X histo from QA
 111 TProfile*  currHDXvsZ = 0;                // raw DX vs Z histo from QA
 112 TProfile*  currHTDvsX = 0;                // raw DriftTime-T0 vs X histo from QA
 113 TProfile*  currHTDvsZ = 0;                // raw DriftTime-T0 vs Z histo from QA
 114 TH1*       currHDXvsXclean = 0;           // smoothed DX vs X histo
 115 TH1*       currHVDvsX = 0;                // extracted VDrift vs X profile
 116 TH1*       currHCorrMapX = 0;             // extracted correction map
 117 TH1*       currHDXvsXCorr = 0;            // DX vs X histo after application of the map
 118 TH1*       currHDXvsZCorrMT0 = 0;         // DX vs X histo after application of the map and t0 correction
 119 TH1*       currHDVvsZ = 0;                // correction for the VDrift vs Z
 120 TH1*       currHDVvsZOrig = 0;            // correction for the VDrift vs Z (before error processing)
 121 TGraphErrors* currGrDxx = 0;              // DX final vs XDrift
 122 TGraphErrors* currGrTDx = 0;              // DX final vs TDrift
 123 TGraphErrors* currGrTXCorr = 0;           // TDrift vs XDritf
 124 TCollection* qaHistos = 0;                // input QA collection
 125 TObjArray procHistos;                     // processed histos buffer
 126 AliITSresponseSDD* sddResp=0;             // SDD response, updated
 127 TObjArray *vdarrayOld = 0;                // olf VDrifts array
 128 TObjArray* vDriftArr=0;                   // VDrifts array, updated
 129 TObjArray* corrMaps=0;                    // holder for correction maps
 130 //
 131 TH1*       resOffsDXraw[2]={0};           // Dx vs X offset at Xdrift=0 raw
 132 TH1*       resOffsDXAMap[2]={0};          // Dx vs X offset (mean) after correction by map
 133 TH1*       resOffsDX[2]={0};              // Dx vs X offset at Xdrift=0 after correction by map (for each side)
 134 TH1*       resVDCorr[2]={0};              // VDrift correction  (for each side)
 135 TH1*       resVDMean[2]={0};              // average VDrift (for each side)
 136 TH1*       resVDCorrZ[2]={0};             // mean VDrift vs Z corrections (for each side)
 137 TH1*       resT0Corr = 0;                 // correction to modules T0
 138 TH1*       resXLocCorr = 0;               // correction to module location
 139 TCanvas*   sddCanv = 0;                   // report canvas
 140 //
 141 Bool_t     sddRespUpdT0OK = kFALSE;       // flag that SDDresponse object T0 was updated
 142 Bool_t     sddRespUpdVDOK = kFALSE;       // flag that SDDresponse object VDrift correction was updated
 143 Bool_t     sddVDriftUpdOK = kFALSE;       // flag that SDD VDrift object was updated
 144 //
 145 TString pathSDDRespOld    = "";
 146 TString pathSDDVDriftOld  = "";
 147 TString pathSDDCorrMapOld = "";
 148 //
 149 //--------------------------------------------------------------------
 150 void      Process(const char* pathSDDResp=0,const char* pathSDDVDrift=0, const char* pathSDDCorrMap=0);
 151 void      Process(TCollection* qa, const char* pathSDDResp=0,const char* pathSDDVDrift=0, const char* pathSDDCorrMap=0);
 152 Bool_t    ProcSDDSideDXvsX();
 153 Bool_t    ProcSDDSideDXvsZ();
 154 TProfile* GetQAHisto(Int_t qaType);
 155 TProfile* H2Profile(TH2* h2);
 156 TH1*      H2ProfileAsTH1(TH2* h2);
 157 TH1*      GetProfHEntries(TProfile* prof);
 158 TH1*      ProfileAsTH1(TProfile* prof, const char* addName);
 159 TH1*      CleanProfile(TProfile* profR);
 160 TH1*      Vdrift2Resid(TH1* vd);
 161 TH1*      Resid2Vdrift(TH1* res);
 162 void      RedoProfileErrors(TH1* profH1,TProfile* prof);
 163 void      SetModuleID(Int_t mdID,Int_t side=-1);
 164 void      PrepareModuleHistos(Int_t mdID, Int_t side, Bool_t vsZ);
 165 void      CheckResultDX();
 166 void      CalcDXCorrections();
 167 double    GetVOld(double z);
 168 Int_t     GetStoreID(int type, int imd=-1,int side=-1);
 169 void      CleanPrev();
 170 void      StoreCurrent();
 171 double    ftVdZ(double *x, double *par);
 172 void      PlotReport(const char* psname="repSDDQA.ps");
 173 TH1*      GetPadBaseHisto(TPad* pad);
 174 Bool_t    PlotHisto(TH1* h, Option_t* opt="", int mrkStyle=20,int mrkCol=kBlack, double mrkSize=1.);
 175 TLatex*   AddPadLabel(const char*txt,float x=0.1,float y=0.9,int color=kBlack,float size=0.04);
 176 void      GetHistoMinMaxInRange(TH1* h, double &mn,double &mx);
 177 double    edgeLow(double* x, double *par);
 178 double    edgeHigh(double* x, double *par);
 179 void      CureEdges(TH1* prof);
 180 void      SafeRebin(TProfile* prof, Int_t factor, Bool_t xprof);
 181 TH1*      FitDXEdges(TProfile* prof);
 182 Bool_t    TestMapFunction(TH1* smap, TF1* fun, double lft, double rgt);
 183 double    ftPolComb(double* x, double *par);
 184 TH1*      SimpleMap(TH1* prof);
 185 //
 186 Bool_t    LoadSDDVDrift(TString& path, TObjArray *&arr);
 187 Bool_t    LoadSDDResponse(TString& path, AliITSresponseSDD *&resp);
 188 Bool_t    LoadSDDCorrMap(TString& path, TObjArray *&maps);
 189 AliCDBEntry* GetCDBEntry(const char* path);
 190 //
 191 void                UpdateSDDResponse(AliITSresponseSDD *resp, Bool_t t0=kTRUE, Bool_t vdrift=kTRUE);
 192 void                UpdateSDDVDrift(AliITSDriftSpeedArraySDD* vdArr, int imd, int side);
 193 TObjArray*          UpdateSDDVDrift();
 194 TObjArray*          CreateCorrMaps();
 195 AliITSresponseSDD*  UpdateSDDResponse(Bool_t t0=kTRUE, Bool_t vdrift=kTRUE);
 196 AliITSCorrMap1DSDD* CreateCorrMap(TH1* mapHisto, int imd, int side, AliITSCorrMap1DSDD* updateMap=0);
 197 void      PrepCDBObj(TObject *obj,const char* path,int firstrun=0,int lastrun=999999999,const char* comment="");
 198
 199 //-------------------------------------------------------------------
 200
 201 //_________________________________________________________________________
 202 void 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 //_________________________________________________________________________
 210 void 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 //_________________________________________________________________________
 256 Bool_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     //
 267     // try simple solution
 268     if (!(currHCorrMapX=SimpleMap(currHDXvsXclean))) {
 269       currHVDvsX      = Resid2Vdrift(currHDXvsXclean);
 270       currHCorrMapX   = Vdrift2Resid(currHVDvsX);
 271     }
 272     //
 273     currHDXvsXCorr->Add(currHCorrMapX,-1);
 274   }
 275   // check results
 276   CheckResultDX();
 277   //
 278   return kTRUE;
 279 }
 280
 281 //_________________________________________________________________________________
 282 Bool_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 //_________________________________________________________________________
 436 void 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 //_________________________________________________________________________
 451 void 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 //_________________________________________________________________________
 460 TProfile* 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   //
 471   int trueSide = userLeftRightFromTASK ? (1-currSDDSide) : currSDDSide;
 472   const char* hname = Form("%s%d_%s",kQAHistoNames[qaType],currSDDId,kLRNames[trueSide]);
 473   TH1* h = (TH1*) qaHistos->FindObject(hname);
 474   if (!h) {printf("Did not find histo %s in QA output collection\n",hname); exit(1);}
 475   //
 476   if (h->InheritsFrom(TH2::Class())) h = H2Profile((TH2*)h); // make sure it is not TH2 histo
 477   //
 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   //
 481   return (TProfile*)h;
 482 }
 483
 484 //_________________________________________________________________________
 485 TH1* 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 //_________________________________________________________________________
 495 TH1* 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();
 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   }
 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 //_________________________________________________________________________
 516 TProfile* 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 //_________________________________________________________________________________
 525 TH1* CleanProfile(TProfile* profR)
 526 {
 527   // clean profile copy
 528   TH1* prof = FitDXEdges(profR); // cure edges
 529
 530   int ib0 = profR->FindBin(1), ib1 = profR->FindBin(sddSeg.Dx()-1);
 531   int nbn = profR->GetNbinsX();
 532   int nSkip = skipDXEdge/profR->GetBinWidth(nbn/2);
 533   int nMean = wDXEdge/profR->GetBinWidth(nbn/2);
 534   //
 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 //_________________________________________________________________________________
 572 TH1* 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   //
 582   int nbcnt=0, nbfl = wVDEdge/vd->GetBinWidth(nb/2);
 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;}
 601     err2 = TMath::Sqrt(err2)/res->GetBinWidth(nb/2);
 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 //_________________________________________________________________________________
 672 TH1* 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();
 679   if (userDummyCorrMap) return res;
 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 //__________________________________________________________________________________
 712 void 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++;
 728   currHDXvsX->Fit(fitP0,"q0N","");
 729   double offsRaw = fitP0->GetParameter(0);
 730   currHDXvsXCorr->Fit(fitP0,"q0N","");
 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   //
 765   fitP1->SetParameters(0,0);
 766   currGrDxx->Fit(fitP1,"q","",xmin, xmax);
 767   double offs = fitP1->GetParameter(0);           // offset of correction line at Xdrift=0
 768   //
 769   //  printf("Fitting VD correction in the range %.1f:%.1f\n",tmin,tmax);
 770   fitP1->SetParameters(0,0);
 771   currGrTDx->Fit(fitP1,"q","",tmin,tmax);
 772   double vcor = fitP1->GetParameter(1);
 773   //
 774   fitP1->SetParameters(0,0);
 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 //__________________________________________________________________________
 818 void 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   //
 853   if (userDummyt0Corr) t0Corr = 0;
 854   if (userDummyxlCorr) xlCorr = 0;
 855   //  printf("SDD%d VL:%f VR:%f offsL:%+f offsR:%+f  dT:%+f dX:%+f\n",currSDDId, vL,vR, offsL,offsR, t0Corr,xlCorr);
 856   resT0Corr->SetBinContent(currMod+1, t0Corr);        // T0 correction
 857   resXLocCorr->SetBinContent(currMod+1, xlCorr);      // X alignment correction
 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 //______________________________________________________________
 891 Int_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 //______________________________________________________________
 905 void CleanPrev()
 906 {
 907   // clean "current" objects from last event
 908   currHDXvsX = 0;
 909   currHDXvsZ = 0;
 910   currHTDvsX = 0;
 911   currHTDvsZ = 0;
 912   currHDXvsXclean = 0;
 913   currHVDvsX = 0;
 914   currHCorrMapX = 0;
 915   currHDXvsXCorr = 0;
 916   currHDVvsZ = 0;
 917   currGrDxx = 0;
 918   currGrTDx = 0;
 919   currGrTXCorr = 0;
 920   //
 921 }
 922
 923 //______________________________________________________________
 924 void 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));
 932   if (currHVDvsX)       procHistos.AddAtAndExpand(currHVDvsX,      GetStoreID(kSVDvsX));
 933   if (currHCorrMapX)    procHistos.AddAtAndExpand(currHCorrMapX,   GetStoreID(kSCorrMapX));
 934   if (currHDXvsXCorr)   procHistos.AddAtAndExpand(currHDXvsXCorr,  GetStoreID(kSXvsXCorr));
 935   if (currHDXvsZCorrMT0) procHistos.AddAtAndExpand(currHDXvsZCorrMT0,  GetStoreID(kSXvsZCorr));
 936   if (currHDVvsZ)       procHistos.AddAtAndExpand(currHDVvsZ,      GetStoreID(kSDVvsZ));
 937   if (currHDVvsZOrig)   procHistos.AddAtAndExpand(currHDVvsZOrig,  GetStoreID(kSDVvsZOrig));
 938   if (currGrDxx)        procHistos.AddAtAndExpand(currGrDxx,       GetStoreID(kSGrDxx));
 939   if (currGrTDx)        procHistos.AddAtAndExpand(currGrTDx,       GetStoreID(kSGrTDx));
 940   if (currGrTXCorr)     procHistos.AddAtAndExpand(currGrTXCorr,    GetStoreID(kSGrTXCorr));
 941   //
 942 }
 943
 944 //_________________________________________________________________________________
 945 TObjArray* 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 //_________________________________________________________________________________
 972 AliITSCorrMap1DSDD* 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   //
 984   // check if the updateMap is meaningful
 985   if (updateMap && updateMap->GetNBinsDrift()>2 && nbCorr>1) {
 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 //_________________________________________________________________________________
1007 TObjArray* 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 //_________________________________________________________________________________
1032 void 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 //_________________________________________________________________________________
1056 AliITSresponseSDD* 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 //_________________________________________________________________________________
1069 void 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 //___________________________________________________________________
1093 double 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 //___________________________________________________________________
1104 double 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 //________________________________________________________________________________________________________
1122 Bool_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 //________________________________________________________________________________________________________
1164 Bool_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 //________________________________________________________________________________________________________
1206 Bool_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 //_______________________________________________________________________________________
1249 AliCDBEntry* 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 //_______________________________________________________________________________________
1279 Bool_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 //_______________________________________________________________________________________
1313 void 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 //_______________________________________________________________________________________
1329 void 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);
1350     PlotHisto(resOffsDXraw[ix],"p"     ,7,kBlack,0.5);
1351     PlotHisto(resOffsDX[ix]   ,"p same",7,kRed  ,1);
1352     AddPadLabel(Form("<#DeltaX> %s : Raw",ix?"Right":"Left"), 0.1,0.93,kBlack,0.05);
1353     AddPadLabel("After Map", 0.5,0.93,kRed,0.05);
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);
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       //
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       }
1411       PlotHisto(hsdd,"p same",20,kBlack,0.4);
1412       PlotHisto(hsddcl,"p sames",24,kGreen+2,0.5);
1413       //
1414       //
1415       hsdd = (TH1*)procHistos.At(GetStoreID(kSCorrMapX,imd,ix)); // map
1416       if (hsdd) hsdd->GetXaxis()->SetRange(ib0,ib1);
1417       PlotHisto(hsdd,"same",7,kRed,0.5);
1418       hsdd = (TH1*)procHistos.At(GetStoreID(kSXvsXCorr,imd,ix));
1419       if (hsdd) hsdd->GetXaxis()->SetRange(ib0,ib1);
1420       PlotHisto(hsdd,"histo same",7,kBlue,0.5);
1421       //
1422       AddPadLabel(Form("<#DeltaX> %d %s: Raw",imd,ix?"Right":"Left"), 0.1,0.93,kBlack,0.07);
1423       AddPadLabel("Clean", 0.35,0.93,kGreen+2,0.07);
1424       AddPadLabel("Map", 0.42,0.93,kRed,0.07);
1425       AddPadLabel("+Map", 0.5,0.93,kBlue,0.07);
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 //__________________________________
1458 TH1* 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 //__________________________________
1474 TLatex* 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 //__________________________________
1485 void 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 //__________________________________
1494 void 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 //__________________________________________________________
1511 double 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 //__________________________________________________________
1540 double 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 //_________________________________________________________________________
1568 void 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 //_________________________________________________________________________
1595 void 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 }
1673
1674 //_________________________________________________________________________
1675 void 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 //____________________________________________________
1747 Double_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 //____________________________________________________
1841 Double_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 //____________________________________________________________
1901 TH1* 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 //___________________________________________________________
1913 TH1* 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
2066 double 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 //______________________________________________________________
2085 TH1* 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 //______________________________________________________________________
2166 Bool_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 }