TRD/AliTRDCalibraFit.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////////////////////////
  19 //
  20 // AliTRDCalibraFit
  21 //
  22 // This class is for the TRD calibration of the relative gain factor, the drift velocity,
  23 // the time 0 and the pad response function. It fits the histos.
  24 // The 2D histograms or vectors (first converted in 1D histos) will be fitted
  25 // if they have enough entries, otherwise the (default) value of the choosen database
  26 // will be put. For the relative gain calibration the resulted factors will be globally
  27 // normalized to the gain factors of the choosen database. It unables to precise
  28 // previous calibration procedure.
  29 // The function SetDebug enables the user to see:
  30 // _fDebug = 0: nothing, only the values are written in the tree if wanted
  31 // _fDebug = 1: a comparaison of the coefficients found and the default values
  32 //              in the choosen database.
  33 //              fCoef , histogram of the coefs as function of the calibration group number
  34 //              fDelta , histogram of the relative difference of the coef with the default
  35 //                        value in the database as function of the calibration group number
  36 //              fError , dirstribution of this relative difference
  37 // _fDebug = 2: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
  38 // _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the
  39 //              pad row and col number
  40 // _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with
  41 //              also the comparaison histograms of the 1 for this detector
  42 //
  43 //
  44 // Author:
  45 //   R. Bailhache (R.Bailhache@gsi.de)
  46 //
  47 //////////////////////////////////////////////////////////////////////////////////////
  48
  49 #include <TLine.h>
  50 #include <TH1I.h>
  51 #include <TStyle.h>
  52 #include <TProfile2D.h>
  53 #include <TCanvas.h>
  54 #include <TGraphErrors.h>
  55 #include <TObjArray.h>
  56 #include <TH1.h>
  57 #include <TH1F.h>
  58 #include <TF1.h>
  59 #include <TAxis.h>
  60 #include <TMath.h>
  61 #include <TDirectory.h>
  62 #include <TROOT.h>
  63 #include <TTreeStream.h>
  64 #include <TLinearFitter.h>
  65 #include <TVectorD.h>
  66 #include <TArrayF.h>
  67
  68 #include "AliLog.h"
  69 #include "AliMathBase.h"
  70
  71 #include "AliTRDCalibraFit.h"
  72 #include "AliTRDCalibraMode.h"
  73 #include "AliTRDCalibraVector.h"
  74 #include "AliTRDCalibraVdriftLinearFit.h"
  75 #include "AliTRDcalibDB.h"
  76 #include "AliTRDgeometry.h"
  77 #include "AliTRDpadPlane.h"
  78 #include "AliTRDgeometry.h"
  79 #include "./Cal/AliTRDCalROC.h"
  80 #include "./Cal/AliTRDCalPad.h"
  81 #include "./Cal/AliTRDCalDet.h"
  82
  83
  84 ClassImp(AliTRDCalibraFit)
  85
  86 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
  87 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
  88
  89 //_____________singleton implementation_________________________________________________
  90 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
  91 {
  92   //
  93   // Singleton implementation
  94   //
  95
  96   if (fgTerminated != kFALSE) {
  97     return 0;
  98   }
  99
 100   if (fgInstance == 0) {
 101     fgInstance = new AliTRDCalibraFit();
 102   }
 103
 104   return fgInstance;
 105
 106 }
 107
 108 //______________________________________________________________________________________
 109 void AliTRDCalibraFit::Terminate()
 110 {
 111   //
 112   // Singleton implementation
 113   // Deletes the instance of this class
 114   //
 115
 116   fgTerminated = kTRUE;
 117
 118   if (fgInstance != 0) {
 119     delete fgInstance;
 120     fgInstance = 0;
 121   }
 122
 123 }
 124
 125 //______________________________________________________________________________________
 126 AliTRDCalibraFit::AliTRDCalibraFit()
 127   :TObject()
 128   ,fGeo(0)
 129   ,fNumberOfBinsExpected(0)
 130   ,fMethod(0)
 131   ,fBeginFitCharge(3.5)
 132   ,fFitPHPeriode(1)
 133   ,fTakeTheMaxPH(kTRUE)
 134   ,fT0Shift0(0.124797)
 135   ,fT0Shift1(0.267451)
 136   ,fRangeFitPRF(1.0)
 137   ,fAccCDB(kFALSE)
 138   ,fMinEntries(800)
 139   ,fRebin(1)
 140   ,fNumberFit(0)
 141   ,fNumberFitSuccess(0)
 142   ,fNumberEnt(0)
 143   ,fStatisticMean(0.0)
 144   ,fDebugStreamer(0x0)
 145   ,fDebugLevel(0)
 146   ,fFitVoir(0)
 147   ,fMagneticField(0.5)
 148   ,fCalibraMode(new AliTRDCalibraMode())
 149   ,fCurrentCoefE(0.0)
 150   ,fCurrentCoefE2(0.0)
 151   ,fDect1(0)
 152   ,fDect2(0)
 153   ,fScaleFitFactor(0.0)
 154   ,fEntriesCurrent(0)
 155   ,fCountDet(0)
 156   ,fCount(0)
 157   ,fCalDet(0x0)
 158   ,fCalROC(0x0)
 159   ,fCalDet2(0x0)
 160   ,fCalROC2(0x0)
 161   ,fCurrentCoefDetector(0x0)
 162   ,fCurrentCoefDetector2(0x0)
 163   ,fVectorFit(0)
 164   ,fVectorFit2(0)
 165 {
 166   //
 167   // Default constructor
 168   //
 169
 170   fGeo         = new AliTRDgeometry();
 171
 172   // Current variables initialised
 173   for (Int_t k = 0; k < 2; k++) {
 174     fCurrentCoef[k]  = 0.0;
 175     fCurrentCoef2[k] = 0.0;
 176   }
 177   for (Int_t i = 0; i < 3; i++) {
 178     fPhd[i]          = 0.0;
 179     fDet[i]          = 0;
 180   }
 181
 182 }
 183 //______________________________________________________________________________________
 184 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
 185 :TObject(c)
 186 ,fGeo(0)
 187 ,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
 188 ,fMethod(c.fMethod)
 189 ,fBeginFitCharge(c.fBeginFitCharge)
 190 ,fFitPHPeriode(c.fFitPHPeriode)
 191 ,fTakeTheMaxPH(c.fTakeTheMaxPH)
 192 ,fT0Shift0(c.fT0Shift0)
 193 ,fT0Shift1(c.fT0Shift1)
 194 ,fRangeFitPRF(c.fRangeFitPRF)
 195 ,fAccCDB(c.fAccCDB)
 196 ,fMinEntries(c.fMinEntries)
 197 ,fRebin(c.fRebin)
 198 ,fNumberFit(c.fNumberFit)
 199 ,fNumberFitSuccess(c.fNumberFitSuccess)
 200 ,fNumberEnt(c.fNumberEnt)
 201 ,fStatisticMean(c.fStatisticMean)
 202 ,fDebugStreamer(0x0)
 203 ,fDebugLevel(c.fDebugLevel)
 204 ,fFitVoir(c.fFitVoir)
 205 ,fMagneticField(c.fMagneticField)
 206 ,fCalibraMode(0x0)
 207 ,fCurrentCoefE(c.fCurrentCoefE)
 208 ,fCurrentCoefE2(c.fCurrentCoefE2)
 209 ,fDect1(c.fDect1)
 210 ,fDect2(c.fDect2)
 211 ,fScaleFitFactor(c.fScaleFitFactor)
 212 ,fEntriesCurrent(c.fEntriesCurrent)
 213 ,fCountDet(c.fCountDet)
 214 ,fCount(c.fCount)
 215 ,fCalDet(0x0)
 216 ,fCalROC(0x0)
 217 ,fCalDet2(0x0)
 218 ,fCalROC2(0x0)
 219 ,fCurrentCoefDetector(0x0)
 220 ,fCurrentCoefDetector2(0x0)
 221 ,fVectorFit(0)
 222 ,fVectorFit2(0)
 223 {
 224   //
 225   // Copy constructor
 226   //
 227
 228   if(c.fCalibraMode)   fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
 229
 230   //Current variables initialised
 231   for (Int_t k = 0; k < 2; k++) {
 232     fCurrentCoef[k]  = 0.0;
 233     fCurrentCoef2[k] = 0.0;
 234   }
 235   for (Int_t i = 0; i < 3; i++) {
 236     fPhd[i]          = 0.0;
 237     fDet[i]          = 0;
 238   }
 239   if(c.fCalDet) fCalDet   = new AliTRDCalDet(*c.fCalDet);
 240   if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
 241
 242   if(c.fCalROC) fCalROC   = new AliTRDCalROC(*c.fCalROC);
 243   if(c.fCalROC2) fCalROC  = new AliTRDCalROC(*c.fCalROC2);
 244
 245   fVectorFit.SetName(c.fVectorFit.GetName());
 246   for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
 247     AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
 248     Int_t detector         = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
 249     Int_t ntotal = 1;
 250     if (GetChamber(detector) == 2) {
 251       ntotal = 1728;
 252     }
 253     else {
 254       ntotal = 2304;
 255     }
 256     Float_t *coef = new Float_t[ntotal];
 257     for (Int_t i = 0; i < ntotal; i++) {
 258       coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
 259     }
 260     fitInfo->SetCoef(coef);
 261     fitInfo->SetDetector(detector);
 262     fVectorFit.Add((TObject *) fitInfo);
 263   }
 264   fVectorFit.SetName(c.fVectorFit.GetName());
 265   for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
 266     AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
 267     Int_t detector         = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
 268     Int_t ntotal = 1;
 269     if (GetChamber(detector) == 2) {
 270       ntotal = 1728;
 271     }
 272     else {
 273       ntotal = 2304;
 274     }
 275     Float_t *coef = new Float_t[ntotal];
 276     for (Int_t i = 0; i < ntotal; i++) {
 277       coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
 278     }
 279     fitInfo->SetCoef(coef);
 280     fitInfo->SetDetector(detector);
 281     fVectorFit2.Add((TObject *) fitInfo);
 282   }
 283   if (fGeo) {
 284     delete fGeo;
 285   }
 286   fGeo = new AliTRDgeometry();
 287
 288 }
 289 //____________________________________________________________________________________
 290 AliTRDCalibraFit::~AliTRDCalibraFit()
 291 {
 292   //
 293   // AliTRDCalibraFit destructor
 294   //
 295   if ( fDebugStreamer ) delete fDebugStreamer;
 296   if ( fCalDet )  delete fCalDet;
 297   if ( fCalDet2 ) delete fCalDet2;
 298   if ( fCalROC )  delete fCalROC;
 299   if ( fCalROC2 ) delete fCalROC2;
 300   fVectorFit.Delete();
 301   fVectorFit2.Delete();
 302   if (fGeo) {
 303     delete fGeo;
 304   }
 305
 306 }
 307 //_____________________________________________________________________________
 308 void AliTRDCalibraFit::Destroy()
 309 {
 310   //
 311   // Delete instance
 312   //
 313
 314   if (fgInstance) {
 315     delete fgInstance;
 316     fgInstance = 0x0;
 317   }
 318
 319 }
 320 //__________________________________________________________________________________
 321 void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end)
 322 {
 323   //
 324   // From the drift velocity and t0
 325   // return the position of the peak and maximum negative slope
 326   //
 327
 328   const Float_t kDrWidth = AliTRDgeometry::DrThick();    // drift region
 329   Double_t widbins = 0.1;                                // 0.1 mus
 330
 331   //peak and maxnegslope in mus
 332   Double_t begind = t0*widbins + fT0Shift0;
 333   Double_t peakd  = t0*widbins + fT0Shift1;
 334   Double_t maxnegslope = (kDrWidth + vdrift*peakd)/vdrift;
 335
 336   // peak and maxnegslope in timebin
 337   begin = TMath::Nint(begind*widbins);
 338   peak  = TMath::Nint(peakd*widbins);
 339   end   = TMath::Nint(maxnegslope*widbins);
 340
 341 }
 342 //____________Functions fit Online CH2d________________________________________
 343 Bool_t AliTRDCalibraFit::AnalyseCH(TH2I *ch)
 344 {
 345   //
 346   // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
 347   // calibration group normalized the resulted coefficients (to 1 normally)
 348   //
 349
 350   // Set the calibration mode
 351   const char *name = ch->GetTitle();
 352   SetModeCalibration(name,0);
 353
 354   // Number of Ybins (detectors or groups of pads)
 355   Int_t    nbins   = ch->GetNbinsX();// charge
 356   Int_t    nybins  = ch->GetNbinsY();// groups number
 357   if (!InitFit(nybins,0)) {
 358     return kFALSE;
 359   }
 360   if (!InitFitCH()) {
 361     return kFALSE;
 362   }
 363   fStatisticMean        = 0.0;
 364   fNumberFit            = 0;
 365   fNumberFitSuccess     = 0;
 366   fNumberEnt            = 0;
 367   // Init fCountDet and fCount
 368   InitfCountDetAndfCount(0);
 369   // Beginning of the loop betwwen dect1 and dect2
 370   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 371     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
 372     UpdatefCountDetAndfCount(idect,0);
 373     ReconstructFitRowMinRowMax(idect, 0);
 374     // Take the histo
 375     TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
 376     projch->SetDirectory(0);
 377     // Number of entries for this calibration group
 378     Double_t nentries = 0.0;
 379     Double_t mean = 0.0;
 380     for (Int_t k = 0; k < nbins; k++) {
 381       Int_t binnb = (nbins+2)*(idect+1)+(k+1);
 382       nentries += ch->GetBinContent(binnb);
 383       mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
 384       projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
 385     }
 386     projch->SetEntries(nentries);
 387     //printf("The number of entries for the group %d is %f\n",idect,nentries);
 388     if (nentries > 0) {
 389       fNumberEnt++;
 390       mean /= nentries;
 391     }
 392     // Rebin and statistic stuff
 393     if (fRebin > 1) {
 394       projch = ReBin((TH1I *) projch);
 395     }
 396     // This detector has not enough statistics or was off
 397     if (nentries <= fMinEntries) {
 398       NotEnoughStatisticCH(idect);
 399       if (fDebugLevel != 1) {
 400         delete projch;
 401       }
 402       continue;
 403     }
 404     // Statistics of the group fitted
 405     fStatisticMean += nentries;
 406     fNumberFit++;
 407     //Method choosen
 408     switch(fMethod)
 409       {
 410       case 0: FitMeanW((TH1 *) projch, nentries); break;
 411       case 1: FitMean((TH1 *) projch, nentries, mean); break;
 412       case 2: FitCH((TH1 *) projch, mean); break;
 413       case 3: FitBisCH((TH1 *) projch, mean); break;
 414       default: return kFALSE;
 415       }
 416     // Fill Infos Fit
 417     FillInfosFitCH(idect);
 418     // Memory!!!
 419     if (fDebugLevel != 1) {
 420       delete projch;
 421     }
 422   } // Boucle object
 423   // Normierungcharge
 424   if (fDebugLevel != 1) {
 425     NormierungCharge();
 426   }
 427   // Mean Statistic
 428   if (fNumberFit > 0) {
 429     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
 430     fStatisticMean = fStatisticMean / fNumberFit;
 431   }
 432   else {
 433     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 434   }
 435   delete fDebugStreamer;
 436   fDebugStreamer = 0x0;
 437
 438   return kTRUE;
 439 }
 440 //____________Functions fit Online CH2d________________________________________
 441 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
 442 {
 443   //
 444   // Reconstruct a 1D histo from the vectorCH for each calibration group,
 445   // fit the histo, normalized the resulted coefficients (to 1 normally)
 446   //
 447
 448   // Set the calibraMode
 449   const char *name = calvect->GetNameCH();
 450   SetModeCalibration(name,0);
 451
 452   // Number of Xbins (detectors or groups of pads)
 453   if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
 454     return kFALSE;
 455   }
 456   if (!InitFitCH()) {
 457     return kFALSE;
 458   }
 459   fStatisticMean        = 0.0;
 460   fNumberFit            = 0;
 461   fNumberFitSuccess     = 0;
 462   fNumberEnt            = 0;
 463   // Init fCountDet and fCount
 464   InitfCountDetAndfCount(0);
 465   // Beginning of the loop between dect1 and dect2
 466   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 467     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
 468     UpdatefCountDetAndfCount(idect,0);
 469     ReconstructFitRowMinRowMax(idect,0);
 470     // Take the histo
 471     Double_t nentries = 0.0;
 472     Double_t mean = 0.0;
 473     TH1F *projch = 0x0;
 474     Bool_t something = kTRUE;
 475     if(!calvect->GetCHEntries(fCountDet)) something = kFALSE;
 476     if(something){
 477       TString name("CH");
 478       name += idect;
 479       projch  = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) name);
 480       projch->SetDirectory(0);
 481       for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
 482         nentries += projch->GetBinContent(k+1);
 483         mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
 484       }
 485       if (nentries > 0) {
 486         fNumberEnt++;
 487         mean /= nentries;
 488       }
 489       //printf("The number of entries for the group %d is %f\n",idect,nentries);
 490       // Rebin
 491       if (fRebin >  1) {
 492         projch = ReBin((TH1F *) projch);
 493       }
 494     }
 495     // This detector has not enough statistics or was not found in VectorCH
 496     if (nentries <= fMinEntries) {
 497       NotEnoughStatisticCH(idect);
 498       if (fDebugLevel != 1) {
 499         if(projch) delete projch;
 500       }
 501       continue;
 502     }
 503     // Statistic of the histos fitted
 504     fStatisticMean += nentries;
 505     fNumberFit++;
 506     //Method choosen
 507     switch(fMethod)
 508       {
 509       case 0: FitMeanW((TH1 *) projch, nentries); break;
 510       case 1: FitMean((TH1 *) projch, nentries, mean); break;
 511       case 2: FitCH((TH1 *) projch, mean); break;
 512       case 3: FitBisCH((TH1 *) projch, mean); break;
 513       default: return kFALSE;
 514       }
 515     // Fill Infos Fit
 516     FillInfosFitCH(idect);
 517     // Memory!!!
 518     if (fDebugLevel != 1) {
 519       delete projch;
 520     }
 521   } // Boucle object
 522   // Normierungcharge
 523   if (fDebugLevel != 1) {
 524     NormierungCharge();
 525   }
 526   // Mean Statistics
 527   if (fNumberFit > 0) {
 528     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
 529     fStatisticMean = fStatisticMean / fNumberFit;
 530   }
 531   else {
 532     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 533   }
 534   delete fDebugStreamer;
 535   fDebugStreamer = 0x0;
 536   return kTRUE;
 537 }
 538 //________________functions fit Online PH2d____________________________________
 539 Bool_t AliTRDCalibraFit::AnalysePH(TProfile2D *ph)
 540 {
 541   //
 542   // Take the 1D profiles (average pulse height), projections of the 2D PH
 543   // on the Xaxis, for each calibration group
 544   // Reconstruct a drift velocity
 545   // A first calibration of T0 is also made  using the same method
 546   //
 547
 548   // Set the calibration mode
 549   const char *name = ph->GetTitle();
 550   SetModeCalibration(name,1);
 551
 552   // Number of Xbins (detectors or groups of pads)
 553   Int_t    nbins   = ph->GetNbinsX();// time
 554   Int_t    nybins  = ph->GetNbinsY();// calibration group
 555   if (!InitFit(nybins,1)) {
 556     return kFALSE;
 557   }
 558   if (!InitFitPH()) {
 559     return kFALSE;
 560   }
 561   fStatisticMean        = 0.0;
 562   fNumberFit            = 0;
 563   fNumberFitSuccess     = 0;
 564   fNumberEnt            = 0;
 565   // Init fCountDet and fCount
 566   InitfCountDetAndfCount(1);
 567   // Beginning of the loop
 568   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 569     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
 570     UpdatefCountDetAndfCount(idect,1);
 571     ReconstructFitRowMinRowMax(idect,1);
 572     // Take the histo
 573     TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
 574     projph->SetDirectory(0);
 575     // Number of entries for this calibration group
 576     Double_t nentries = 0;
 577     for (Int_t k = 0; k < nbins; k++) {
 578       Int_t binnb = (nbins+2)*(idect+1)+(k+1);
 579       nentries += ph->GetBinEntries(binnb);
 580     }
 581     if (nentries > 0) {
 582       fNumberEnt++;
 583     }
 584     //printf("The number of entries for the group %d is %f\n",idect,nentries);
 585     // This detector has not enough statistics or was off
 586     if (nentries  <= fMinEntries) {
 587       //printf("Not enough statistic!\n");
 588       NotEnoughStatisticPH(idect);
 589       if (fDebugLevel != 1) {
 590         delete projph;
 591       }
 592       continue;
 593     }
 594     // Statistics of the histos fitted
 595     fNumberFit++;
 596     fStatisticMean += nentries;
 597     // Calcul of "real" coef
 598     CalculVdriftCoefMean();
 599     CalculT0CoefMean();
 600     //Method choosen
 601     switch(fMethod)
 602       {
 603       case 0: FitLagrangePoly((TH1 *) projph); break;
 604       case 1: FitPente((TH1 *) projph); break;
 605       case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
 606       default: return kFALSE;
 607       }
 608     // Fill the tree if end of a detector or only the pointer to the branch!!!
 609     FillInfosFitPH(idect);
 610     // Memory!!!
 611     if (fDebugLevel != 1) {
 612       delete projph;
 613     }
 614   } // Boucle object
 615   // Mean Statistic
 616   if (fNumberFit > 0) {
 617     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
 618     fStatisticMean = fStatisticMean / fNumberFit;
 619   }
 620   else {
 621     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 622   }
 623   delete fDebugStreamer;
 624   fDebugStreamer = 0x0;
 625   return kTRUE;
 626 }
 627 //____________Functions fit Online PH2d________________________________________
 628 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
 629 {
 630   //
 631   // Reconstruct the average pulse height from the vectorPH for each
 632   // calibration group
 633   // Reconstruct a drift velocity
 634   // A first calibration of T0 is also made  using the same method (slope method)
 635   //
 636
 637   // Set the calibration mode
 638   const char *name = calvect->GetNamePH();
 639   SetModeCalibration(name,1);
 640
 641   // Number of Xbins (detectors or groups of pads)
 642   if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
 643     return kFALSE;
 644   }
 645   if (!InitFitPH()) {
 646     return kFALSE;
 647   }
 648   fStatisticMean        = 0.0;
 649   fNumberFit            = 0;
 650   fNumberFitSuccess     = 0;
 651   fNumberEnt            = 0;
 652   // Init fCountDet and fCount
 653   InitfCountDetAndfCount(1);
 654   // Beginning of the loop
 655   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 656     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
 657     UpdatefCountDetAndfCount(idect,1);
 658     ReconstructFitRowMinRowMax(idect,1);
 659     // Take the histo
 660     TH1F *projph = 0x0;
 661     fEntriesCurrent = 0;
 662     Bool_t something = kTRUE;
 663     if(!calvect->GetPHEntries(fCountDet)) something = kFALSE;
 664     if(something){
 665       TString name("PH");
 666       name += idect;
 667       projph  = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) name)));
 668       projph->SetDirectory(0);
 669     }
 670     //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
 671     // This detector has not enough statistics or was off
 672     if (fEntriesCurrent <=  fMinEntries) {
 673       //printf("Not enough stat!\n");
 674       NotEnoughStatisticPH(idect);
 675       if (fDebugLevel != 1) {
 676         if(projph) delete projph;
 677       }
 678       continue;
 679     }
 680     // Statistic of the histos fitted
 681     fNumberFit++;
 682     fStatisticMean += fEntriesCurrent;
 683     // Calcul of "real" coef
 684     CalculVdriftCoefMean();
 685     CalculT0CoefMean();
 686     //Method choosen
 687     switch(fMethod)
 688       {
 689       case 0: FitLagrangePoly((TH1 *) projph); break;
 690       case 1: FitPente((TH1 *) projph); break;
 691       case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
 692       default: return kFALSE;
 693       }
 694     // Fill the tree if end of a detector or only the pointer to the branch!!!
 695     FillInfosFitPH(idect);
 696     // Memory!!!
 697     if (fDebugLevel != 1) {
 698       delete projph;
 699     }
 700   } // Boucle object
 701
 702   // Mean Statistic
 703   if (fNumberFit > 0) {
 704     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
 705     fStatisticMean = fStatisticMean / fNumberFit;
 706   }
 707   else {
 708     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 709   }
 710   delete fDebugStreamer;
 711   fDebugStreamer = 0x0;
 712   return kTRUE;
 713 }
 714 //____________Functions fit Online PRF2d_______________________________________
 715 Bool_t AliTRDCalibraFit::AnalysePRF(TProfile2D *prf)
 716 {
 717   //
 718   // Take the 1D profiles (pad response function), projections of the 2D PRF
 719   // on the Xaxis, for each calibration group
 720   // Fit with a gaussian to reconstruct the sigma of the pad response function
 721   //
 722
 723   // Set the calibration mode
 724   const char *name = prf->GetTitle();
 725   SetModeCalibration(name,2);
 726
 727   // Number of Ybins (detectors or groups of pads)
 728   Int_t    nybins  = prf->GetNbinsY();// calibration groups
 729   Int_t    nbins   = prf->GetNbinsX();// bins
 730   Int_t    nbg     = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
 731   if((nbg > 0) || (nbg == -1)) return kFALSE;
 732   if (!InitFit(nybins,2)) {
 733     return kFALSE;
 734   }
 735   if (!InitFitPRF()) {
 736     return kFALSE;
 737   }
 738   fStatisticMean        = 0.0;
 739   fNumberFit            = 0;
 740   fNumberFitSuccess     = 0;
 741   fNumberEnt            = 0;
 742   // Init fCountDet and fCount
 743   InitfCountDetAndfCount(2);
 744   // Beginning of the loop
 745   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 746     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
 747     UpdatefCountDetAndfCount(idect,2);
 748     ReconstructFitRowMinRowMax(idect,2);
 749     // Take the histo
 750     TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
 751     projprf->SetDirectory(0);
 752     // Number of entries for this calibration group
 753     Double_t nentries = 0;
 754     for (Int_t k = 0; k < nbins; k++) {
 755       Int_t binnb = (nbins+2)*(idect+1)+(k+1);
 756       nentries += prf->GetBinEntries(binnb);
 757     }
 758     if(nentries > 0) fNumberEnt++;
 759     // This detector has not enough statistics or was off
 760     if (nentries <= fMinEntries) {
 761       NotEnoughStatisticPRF(idect);
 762       if (fDebugLevel != 1) {
 763         delete projprf;
 764       }
 765       continue;
 766     }
 767     // Statistics of the histos fitted
 768     fNumberFit++;
 769     fStatisticMean += nentries;
 770     // Calcul of "real" coef
 771     CalculPRFCoefMean();
 772     //Method choosen
 773     switch(fMethod)
 774       {
 775       case 0: FitPRF((TH1 *) projprf); break;
 776       case 1: RmsPRF((TH1 *) projprf); break;
 777       default: return kFALSE;
 778       }
 779     // Fill the tree if end of a detector or only the pointer to the branch!!!
 780     FillInfosFitPRF(idect);
 781     // Memory!!!
 782     if (fDebugLevel != 1) {
 783       delete projprf;
 784     }
 785   } // Boucle object
 786   // Mean Statistic
 787   if (fNumberFit > 0) {
 788     AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
 789     AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
 790     AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
 791                 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
 792     fStatisticMean = fStatisticMean / fNumberFit;
 793   }
 794   else {
 795     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 796   }
 797   delete fDebugStreamer;
 798   fDebugStreamer = 0x0;
 799   return kTRUE;
 800 }
 801 //____________Functions fit Online PRF2d_______________________________________
 802 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(TProfile2D *prf)
 803 {
 804   //
 805   // Take the 1D profiles (pad response function), projections of the 2D PRF
 806   // on the Xaxis, for each calibration group
 807   // Fit with a gaussian to reconstruct the sigma of the pad response function
 808   //
 809
 810   // Set the calibration mode
 811   const char *name = prf->GetTitle();
 812   SetModeCalibration(name,2);
 813
 814   // Number of Ybins (detectors or groups of pads)
 815   TAxis   *xprf    = prf->GetXaxis();
 816   TAxis   *yprf    = prf->GetYaxis();
 817   Int_t    nybins  = yprf->GetNbins();// calibration groups
 818   Int_t    nbins   = xprf->GetNbins();// bins
 819   Float_t  lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
 820   Float_t  upedge  = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
 821   Int_t    nbg     = GetNumberOfGroupsPRF((const char *)name);
 822   if(nbg == -1) return kFALSE;
 823   if(nbg > 0) fMethod = 1;
 824   else fMethod = 0;
 825   if (!InitFit(nybins,2)) {
 826     return kFALSE;
 827   }
 828   if (!InitFitPRF()) {
 829     return kFALSE;
 830   }
 831   fStatisticMean        = 0.0;
 832   fNumberFit            = 0;
 833   fNumberFitSuccess     = 0;
 834   fNumberEnt            = 0;
 835   // Init fCountDet and fCount
 836   InitfCountDetAndfCount(2);
 837   // Beginning of the loop
 838   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 839     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
 840     UpdatefCountDetAndfCount(idect,2);
 841     ReconstructFitRowMinRowMax(idect,2);
 842     // Build the array of entries and sum
 843     TArrayD arraye  = TArrayD(nbins);
 844     TArrayD arraym  = TArrayD(nbins);
 845     TArrayD arrayme = TArrayD(nbins);
 846     Double_t nentries = 0;
 847     //printf("nbins %d\n",nbins);
 848     for (Int_t k = 0; k < nbins; k++) {
 849       Int_t binnb = (nbins+2)*(idect+1)+(k+1);
 850       Double_t entries = (Double_t)prf->GetBinEntries(binnb);
 851       Double_t mean    = (Double_t)prf->GetBinContent(binnb);
 852       Double_t error   = (Double_t)prf->GetBinError(binnb);
 853       //printf("for %d we have %f\n",k,entries);
 854       nentries += entries;
 855       arraye.AddAt(entries,k);
 856       arraym.AddAt(mean,k);
 857       arrayme.AddAt(error,k);
 858     }
 859     if(nentries > 0) fNumberEnt++;
 860     //printf("The number of entries for the group %d is %f\n",idect,nentries);
 861     // This detector has not enough statistics or was off
 862     if (nentries <= fMinEntries) {
 863       NotEnoughStatisticPRF(idect);
 864       continue;
 865     }
 866     // Statistics of the histos fitted
 867     fNumberFit++;
 868     fStatisticMean += nentries;
 869     // Calcul of "real" coef
 870     CalculPRFCoefMean();
 871     //Method choosen
 872     switch(fMethod)
 873       {
 874       case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
 875       case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
 876       default: return kFALSE;
 877       }
 878     // Fill the tree if end of a detector or only the pointer to the branch!!!
 879     FillInfosFitPRF(idect);
 880   } // Boucle object
 881   // Mean Statistic
 882   if (fNumberFit > 0) {
 883     AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
 884     AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
 885     AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
 886                 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
 887     fStatisticMean = fStatisticMean / fNumberFit;
 888   }
 889   else {
 890     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 891   }
 892   delete fDebugStreamer;
 893   fDebugStreamer = 0x0;
 894   return kTRUE;
 895 }
 896 //____________Functions fit Online PRF2d_______________________________________
 897 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
 898 {
 899   //
 900   // Reconstruct the 1D histo (pad response function) from the vectorPRD for
 901   // each calibration group
 902   // Fit with a gaussian to reconstruct the sigma of the pad response function
 903   //
 904
 905   // Set the calibra mode
 906   const char *name = calvect->GetNamePRF();
 907   SetModeCalibration(name,2);
 908   //printf("test0 %s\n",name);
 909
 910   // Number of Xbins (detectors or groups of pads)
 911   if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
 912     //printf("test1\n");
 913     return kFALSE;
 914   }
 915   if (!InitFitPRF()) {
 916     ///printf("test2\n");
 917     return kFALSE;
 918   }
 919   fStatisticMean        = 0.0;
 920   fNumberFit            = 0;
 921   fNumberFitSuccess     = 0;
 922   fNumberEnt            = 0;
 923   // Init fCountDet and fCount
 924   InitfCountDetAndfCount(2);
 925   // Beginning of the loop
 926   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 927     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
 928     UpdatefCountDetAndfCount(idect,2);
 929     ReconstructFitRowMinRowMax(idect,2);
 930     // Take the histo
 931     TH1F *projprf = 0x0;
 932     fEntriesCurrent = 0;
 933     Bool_t something = kTRUE;
 934     if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE;
 935     if(something){
 936       TString name("PRF");
 937       name += idect;
 938       projprf  = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) name)));
 939       projprf->SetDirectory(0);
 940     }
 941     // This detector has not enough statistics or was off
 942     if (fEntriesCurrent <= fMinEntries) {
 943       NotEnoughStatisticPRF(idect);
 944       if (fDebugLevel != 1) {
 945         if(projprf) delete projprf;
 946       }
 947       continue;
 948     }
 949     // Statistic of the histos fitted
 950     fNumberFit++;
 951     fStatisticMean += fEntriesCurrent;
 952     // Calcul of "real" coef
 953     CalculPRFCoefMean();
 954     //Method choosen
 955     switch(fMethod)
 956       {
 957       case 1: FitPRF((TH1 *) projprf); break;
 958       case 2: RmsPRF((TH1 *) projprf); break;
 959       default: return kFALSE;
 960       }
 961     // Fill the tree if end of a detector or only the pointer to the branch!!!
 962     FillInfosFitPRF(idect);
 963     // Memory!!!
 964     if (fDebugLevel != 1) {
 965       delete projprf;
 966     }
 967   } // Boucle object
 968   // Mean Statistics
 969   if (fNumberFit > 0) {
 970     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
 971   }
 972   else {
 973     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 974   }
 975   delete fDebugStreamer;
 976   fDebugStreamer = 0x0;
 977   return kTRUE;
 978 }
 979 //____________Functions fit Online PRF2d_______________________________________
 980 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
 981 {
 982   //
 983   // Reconstruct the 1D histo (pad response function) from the vectorPRD for
 984   // each calibration group
 985   // Fit with a gaussian to reconstruct the sigma of the pad response function
 986   //
 987
 988   // Set the calibra mode
 989   const char *name = calvect->GetNamePRF();
 990   SetModeCalibration(name,2);
 991   //printf("test0 %s\n",name);
 992   Int_t    nbg     = GetNumberOfGroupsPRF((const char *)name);
 993   printf("test1 %d\n",nbg);
 994   if(nbg == -1) return kFALSE;
 995   if(nbg > 0) fMethod = 1;
 996   else fMethod = 0;
 997   // Number of Xbins (detectors or groups of pads)
 998   if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
 999     //printf("test2\n");
1000     return kFALSE;
1001   }
1002   if (!InitFitPRF()) {
1003     //printf("test3\n");
1004     return kFALSE;
1005   }
1006   fStatisticMean        = 0.0;
1007   fNumberFit            = 0;
1008   fNumberFitSuccess     = 0;
1009   fNumberEnt            = 0;
1010   // Variables
1011   Int_t nbins           = 0;
1012   Double_t *arrayx       = 0;
1013   Double_t *arraye       = 0;
1014   Double_t *arraym       = 0;
1015   Double_t *arrayme      = 0;
1016   Float_t lowedge       = 0.0;
1017   Float_t upedge        = 0.0;
1018   // Init fCountDet and fCount
1019   InitfCountDetAndfCount(2);
1020   // Beginning of the loop
1021   for (Int_t idect = fDect1; idect < fDect2; idect++) {
1022     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1023     UpdatefCountDetAndfCount(idect,2);
1024     ReconstructFitRowMinRowMax(idect,2);
1025     // Take the histo
1026     TGraphErrors *projprftree = 0x0;
1027     fEntriesCurrent  = 0;
1028     Bool_t something = kTRUE;
1029     if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE;
1030     if(something){
1031       TString name("PRF");
1032       name += idect;
1033       projprftree  = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) name);
1034       nbins   = projprftree->GetN();
1035       arrayx  = (Double_t *)projprftree->GetX();
1036       arraye  = (Double_t *)projprftree->GetEX();
1037       arraym  = (Double_t *)projprftree->GetY();
1038       arrayme = (Double_t *)projprftree->GetEY();
1039       Float_t step = arrayx[1]-arrayx[0];
1040       lowedge = arrayx[0] - step/2.0;
1041       upedge  = arrayx[(nbins-1)] + step/2.0;
1042       //printf("nbins est %d\n",nbins);
1043       for(Int_t k = 0; k < nbins; k++){
1044         fEntriesCurrent += (Int_t)arraye[k];
1045         //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1046         if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1047       }
1048       if(fEntriesCurrent > 0) fNumberEnt++;
1049     }
1050     //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1051     // This detector has not enough statistics or was off
1052     if (fEntriesCurrent <= fMinEntries) {
1053       NotEnoughStatisticPRF(idect);
1054       if(projprftree) delete projprftree;
1055       continue;
1056     }
1057     // Statistic of the histos fitted
1058     fNumberFit++;
1059     fStatisticMean += fEntriesCurrent;
1060     // Calcul of "real" coef
1061     CalculPRFCoefMean();
1062     //Method choosen
1063     switch(fMethod)
1064       {
1065       case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1066       case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1067       default: return kFALSE;
1068       }
1069     // Fill the tree if end of a detector or only the pointer to the branch!!!
1070     FillInfosFitPRF(idect);
1071     // Memory!!!
1072     if (fDebugLevel != 1) {
1073       delete projprftree;
1074     }
1075   } // Boucle object
1076   // Mean Statistics
1077   if (fNumberFit > 0) {
1078     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1079   }
1080   else {
1081     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1082   }
1083   delete fDebugStreamer;
1084   fDebugStreamer = 0x0;
1085   return kTRUE;
1086 }
1087 //____________Functions fit Online CH2d________________________________________
1088 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1089 {
1090   //
1091   // The linear method
1092   //
1093
1094   fStatisticMean        = 0.0;
1095   fNumberFit            = 0;
1096   fNumberFitSuccess     = 0;
1097   fNumberEnt            = 0;
1098   if(!InitFitLinearFitter()) return kFALSE;
1099
1100
1101   for(Int_t idet = 0; idet < 540; idet++){
1102
1103
1104     //printf("detector number %d\n",idet);
1105
1106     // Take the result
1107     TVectorD param(2);
1108     TVectorD error(3);
1109     fEntriesCurrent = 0;
1110     fCountDet       = idet;
1111     Bool_t here     = calivdli->GetParam(idet,&param);
1112     Bool_t heree    = calivdli->GetError(idet,&error);
1113     //printf("here %d and heree %d\n",here, heree);
1114     if(heree) {
1115       fEntriesCurrent = (Int_t) error[2];
1116       fNumberEnt++;
1117     }
1118     //printf("Number of entries %d\n",fEntriesCurrent);
1119     // Nothing found or not enough statistic
1120     if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1121       NotEnoughStatisticLinearFitter();
1122       continue;
1123     }
1124     //param.Print();
1125     //error.Print();
1126     //Statistics
1127     fNumberFit++;
1128     fStatisticMean += fEntriesCurrent;
1129
1130     // Check the fit
1131     if((-(param[1])) <= 0.0) {
1132       NotEnoughStatisticLinearFitter();
1133       continue;
1134     }
1135
1136     // CalculDatabaseVdriftandTan
1137     CalculVdriftLorentzCoef();
1138
1139     // Statistics
1140     fNumberFitSuccess ++;
1141
1142     // Put the fCurrentCoef
1143     fCurrentCoef[0]  = -param[1];
1144     // here the database must be the one of the reconstruction for the lorentz angle....
1145     fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1146     fCurrentCoefE    = error[1];
1147     fCurrentCoefE2   = error[0];
1148     if((fCurrentCoef2[0] != 0.0) && (param[0] != 0.0)){
1149       fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1150     }
1151
1152     // Fill
1153     FillInfosFitLinearFitter();
1154
1155
1156   }
1157   // Mean Statistics
1158   if (fNumberFit > 0) {
1159     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1160   }
1161   else {
1162     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1163   }
1164   delete fDebugStreamer;
1165   fDebugStreamer = 0x0;
1166   return kTRUE;
1167
1168 }
1169 //____________Functions for seeing if the pad is really okey___________________
1170 //_____________________________________________________________________________
1171 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle)
1172 {
1173   //
1174   // Get numberofgroupsprf
1175   //
1176
1177   // Some patterns
1178   const Char_t *pattern0 = "Ngp0";
1179   const Char_t *pattern1 = "Ngp1";
1180   const Char_t *pattern2 = "Ngp2";
1181   const Char_t *pattern3 = "Ngp3";
1182   const Char_t *pattern4 = "Ngp4";
1183   const Char_t *pattern5 = "Ngp5";
1184   const Char_t *pattern6 = "Ngp6";
1185
1186   // Nrphi mode
1187   if (strstr(nametitle,pattern0)) {
1188     return 0;
1189   }
1190   if (strstr(nametitle,pattern1)) {
1191     return 1;
1192   }
1193   if (strstr(nametitle,pattern2)) {
1194     return 2;
1195   }
1196   if (strstr(nametitle,pattern3)) {
1197     return 3;
1198   }
1199   if (strstr(nametitle,pattern4)) {
1200     return 4;
1201   }
1202   if (strstr(nametitle,pattern5)) {
1203     return 5;
1204   }
1205   if (strstr(nametitle,pattern6)){
1206     return 6;
1207   }
1208   else return -1;
1209
1210
1211 }
1212 //_____________________________________________________________________________
1213 Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i)
1214 {
1215   //
1216   // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1217   // corresponding to the given name
1218   //
1219
1220   if(!SetNzFromTObject(name,i)) return kFALSE;
1221   if(!SetNrphiFromTObject(name,i)) return kFALSE;
1222
1223   return kTRUE;
1224
1225 }
1226 //_____________________________________________________________________________
1227 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i)
1228 {
1229   //
1230   // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1231   // corresponding to the given TObject
1232   //
1233
1234   // Some patterns
1235   const Char_t *patternrphi0 = "Nrphi0";
1236   const Char_t *patternrphi1 = "Nrphi1";
1237   const Char_t *patternrphi2 = "Nrphi2";
1238   const Char_t *patternrphi3 = "Nrphi3";
1239   const Char_t *patternrphi4 = "Nrphi4";
1240   const Char_t *patternrphi5 = "Nrphi5";
1241   const Char_t *patternrphi6 = "Nrphi6";
1242
1243   // Nrphi mode
1244   if (strstr(name,patternrphi0)) {
1245     fCalibraMode->SetNrphi(i ,0);
1246     return kTRUE;
1247   }
1248   if (strstr(name,patternrphi1)) {
1249     fCalibraMode->SetNrphi(i, 1);
1250     return kTRUE;
1251   }
1252   if (strstr(name,patternrphi2)) {
1253     fCalibraMode->SetNrphi(i, 2);
1254     return kTRUE;
1255   }
1256   if (strstr(name,patternrphi3)) {
1257     fCalibraMode->SetNrphi(i, 3);
1258     return kTRUE;
1259   }
1260   if (strstr(name,patternrphi4)) {
1261     fCalibraMode->SetNrphi(i, 4);
1262     return kTRUE;
1263   }
1264   if (strstr(name,patternrphi5)) {
1265     fCalibraMode->SetNrphi(i, 5);
1266     return kTRUE;
1267   }
1268   if (strstr(name,patternrphi6)) {
1269     fCalibraMode->SetNrphi(i, 6);
1270     return kTRUE;
1271   }
1272
1273   fCalibraMode->SetNrphi(i ,0);
1274   return kFALSE;
1275
1276 }
1277 //_____________________________________________________________________________
1278 Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i)
1279 {
1280   //
1281   // Set fNz[i] of the AliTRDCalibraFit::Instance()
1282   // corresponding to the given TObject
1283   //
1284
1285   // Some patterns
1286   const Char_t *patternz0    = "Nz0";
1287   const Char_t *patternz1    = "Nz1";
1288   const Char_t *patternz2    = "Nz2";
1289   const Char_t *patternz3    = "Nz3";
1290   const Char_t *patternz4    = "Nz4";
1291
1292   if (strstr(name,patternz0)) {
1293     fCalibraMode->SetNz(i, 0);
1294     return kTRUE;
1295   }
1296   if (strstr(name,patternz1)) {
1297     fCalibraMode->SetNz(i ,1);
1298     return kTRUE;
1299   }
1300   if (strstr(name,patternz2)) {
1301     fCalibraMode->SetNz(i ,2);
1302     return kTRUE;
1303   }
1304   if (strstr(name,patternz3)) {
1305     fCalibraMode->SetNz(i ,3);
1306     return kTRUE;
1307   }
1308   if (strstr(name,patternz4)) {
1309     fCalibraMode->SetNz(i ,4);
1310     return kTRUE;
1311   }
1312
1313   fCalibraMode->SetNz(i ,0);
1314   return kFALSE;
1315 }
1316 //_____________________________________________________________________________
1317 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(TObjArray *vectorFit, Bool_t perdetector)
1318 {
1319   //
1320   // It creates the AliTRDCalDet object from the AliTRDFitInfo
1321   // It takes the mean value of the coefficients per detector
1322   // This object has to be written in the database
1323   //
1324
1325   // Create the DetObject
1326   AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1327
1328   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1329   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1330   Int_t detector = -1;
1331   Float_t value  = 0.0;
1332
1333   for (Int_t k = 0; k < loop; k++) {
1334     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1335     Float_t mean  = 0.0;
1336     if(perdetector){
1337       mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
1338     }
1339     else {
1340       Int_t   count = 0;
1341       Int_t rowMax    = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1342       Int_t colMax    = fGeo->GetColMax(GetPlane(detector));
1343       for (Int_t row = 0; row < rowMax; row++) {
1344         for (Int_t col = 0; col < colMax; col++) {
1345           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1346           mean += TMath::Abs(value);
1347           count++;
1348         } // Col
1349       } // Row
1350       if(count > 0) mean = mean/count;
1351     }
1352     object->SetValue(detector,mean);
1353   }
1354
1355   return object;
1356 }
1357 //_____________________________________________________________________________
1358 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, Bool_t perdetector)
1359 {
1360   //
1361   // It creates the AliTRDCalDet object from the AliTRDFitInfo
1362   // It takes the mean value of the coefficients per detector
1363   // This object has to be written in the database
1364   //
1365
1366   // Create the DetObject
1367   AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1368
1369
1370   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1371   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1372   Int_t detector = -1;
1373   Float_t value  = 0.0;
1374
1375   for (Int_t k = 0; k < loop; k++) {
1376     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1377     Float_t mean  = 0.0;
1378     if(perdetector){
1379       value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1380       if(value > 0) value = value*scaleFitFactor;
1381       mean = TMath::Abs(value);
1382     }
1383     else{
1384       Int_t   count = 0;
1385       Int_t rowMax    = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1386       Int_t colMax    = fGeo->GetColMax(GetPlane(detector));
1387       for (Int_t row = 0; row < rowMax; row++) {
1388         for (Int_t col = 0; col < colMax; col++) {
1389           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1390           if(value > 0) value = value*scaleFitFactor;
1391           mean += TMath::Abs(value);
1392           count++;
1393         } // Col
1394       } // Row
1395       if(count > 0) mean = mean/count;
1396     }
1397     object->SetValue(detector,mean);
1398   }
1399
1400   return object;
1401 }
1402 //_____________________________________________________________________________
1403 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(TObjArray *vectorFit, Bool_t perdetector)
1404 {
1405   //
1406   // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1407   // It takes the min value of the coefficients per detector
1408   // This object has to be written in the database
1409   //
1410
1411   // Create the DetObject
1412   AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1413
1414   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1415   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1416   Int_t detector = -1;
1417   Float_t value  = 0.0;
1418
1419   for (Int_t k = 0; k < loop; k++) {
1420     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1421     Float_t min  = 100.0;
1422     if(perdetector){
1423       min = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1424     }
1425     else{
1426       Int_t rowMax    = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1427       Int_t colMax    = fGeo->GetColMax(GetPlane(detector));
1428       for (Int_t row = 0; row < rowMax; row++) {
1429         for (Int_t col = 0; col < colMax; col++) {
1430           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1431           if(min > value) min = value;
1432         } // Col
1433       } // Row
1434     }
1435     object->SetValue(detector,min);
1436   }
1437
1438   return object;
1439
1440 }
1441 //_____________________________________________________________________________
1442 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(TObjArray *vectorFit)
1443 {
1444   //
1445   // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1446   // It takes the min value of the coefficients per detector
1447   // This object has to be written in the database
1448   //
1449
1450   // Create the DetObject
1451   AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
1452
1453
1454   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1455   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1456   Int_t detector = -1;
1457   Float_t value  = 0.0;
1458
1459   for (Int_t k = 0; k < loop; k++) {
1460     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1461     /*
1462       Int_t rowMax    = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1463       Int_t colMax    = fGeo->GetColMax(GetPlane(detector));
1464       Float_t min  = 100.0;
1465       for (Int_t row = 0; row < rowMax; row++) {
1466       for (Int_t col = 0; col < colMax; col++) {
1467       value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1468       mean += -TMath::Abs(value);
1469       count++;
1470       } // Col
1471       } // Row
1472       if(count > 0) mean = mean/count;
1473     */
1474     value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1475     object->SetValue(detector,-TMath::Abs(value));
1476   }
1477
1478   return object;
1479
1480 }
1481 //_____________________________________________________________________________
1482 TObject *AliTRDCalibraFit::CreatePadObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, AliTRDCalDet *detobject)
1483 {
1484   //
1485   // It Creates the AliTRDCalPad object from AliTRDFitInfo
1486   // You need first to create the object for the detectors,
1487   // where the mean value is put.
1488   // This object has to be written in the database
1489   //
1490
1491   // Create the DetObject
1492   AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
1493
1494   if(!vectorFit){
1495     for(Int_t k = 0; k < 540; k++){
1496       AliTRDCalROC *calROC = object->GetCalROC(k);
1497       Int_t nchannels = calROC->GetNchannels();
1498       for(Int_t ch = 0; ch < nchannels; ch++){
1499         calROC->SetValue(ch,1.0);
1500       }
1501     }
1502   }
1503   else{
1504
1505     Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1506     if(loop != 540) AliInfo("The Vector Fit is not complete!");
1507     Int_t detector = -1;
1508     Float_t value  = 0.0;
1509
1510     for (Int_t k = 0; k < loop; k++) {
1511       detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1512       AliTRDCalROC *calROC = object->GetCalROC(detector);
1513       Float_t mean         = detobject->GetValue(detector);
1514       if(mean == 0) continue;
1515       Int_t rowMax    = calROC->GetNrows();
1516       Int_t colMax    = calROC->GetNcols();
1517       for (Int_t row = 0; row < rowMax; row++) {
1518         for (Int_t col = 0; col < colMax; col++) {
1519           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1520           if(value > 0) value = value*scaleFitFactor;
1521           calROC->SetValue(col,row,TMath::Abs(value)/mean);
1522         } // Col
1523       } // Row
1524     }
1525   }
1526
1527   return object;
1528 }
1529 //_____________________________________________________________________________
1530 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(TObjArray *vectorFit, AliTRDCalDet *detobject)
1531 {
1532   //
1533   // It Creates the AliTRDCalPad object from AliTRDFitInfo
1534   // You need first to create the object for the detectors,
1535   // where the mean value is put.
1536   // This object has to be written in the database
1537   //
1538
1539   // Create the DetObject
1540   AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
1541
1542   if(!vectorFit){
1543     for(Int_t k = 0; k < 540; k++){
1544       AliTRDCalROC *calROC = object->GetCalROC(k);
1545       Int_t nchannels = calROC->GetNchannels();
1546       for(Int_t ch = 0; ch < nchannels; ch++){
1547         calROC->SetValue(ch,1.0);
1548       }
1549     }
1550   }
1551   else {
1552
1553     Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1554     if(loop != 540) AliInfo("The Vector Fit is not complete!");
1555     Int_t detector = -1;
1556     Float_t value  = 0.0;
1557
1558     for (Int_t k = 0; k < loop; k++) {
1559       detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1560       AliTRDCalROC *calROC = object->GetCalROC(detector);
1561       Float_t mean         = detobject->GetValue(detector);
1562       if(mean == 0) continue;
1563       Int_t rowMax    = calROC->GetNrows();
1564       Int_t colMax    = calROC->GetNcols();
1565       for (Int_t row = 0; row < rowMax; row++) {
1566         for (Int_t col = 0; col < colMax; col++) {
1567           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1568           calROC->SetValue(col,row,TMath::Abs(value)/mean);
1569         } // Col
1570       } // Row
1571     }
1572   }
1573   return object;
1574
1575 }
1576 //_____________________________________________________________________________
1577 TObject *AliTRDCalibraFit::CreatePadObjectT0(TObjArray *vectorFit, AliTRDCalDet *detobject)
1578 {
1579   //
1580   // It Creates the AliTRDCalPad object from AliTRDFitInfo2
1581   // You need first to create the object for the detectors,
1582   // where the mean value is put.
1583   // This object has to be written in the database
1584   //
1585
1586   // Create the DetObject
1587   AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
1588
1589   if(!vectorFit){
1590     for(Int_t k = 0; k < 540; k++){
1591       AliTRDCalROC *calROC = object->GetCalROC(k);
1592       Int_t nchannels = calROC->GetNchannels();
1593       for(Int_t ch = 0; ch < nchannels; ch++){
1594         calROC->SetValue(ch,0.0);
1595       }
1596     }
1597   }
1598   else {
1599
1600     Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1601     if(loop != 540) AliInfo("The Vector Fit is not complete!");
1602     Int_t detector = -1;
1603     Float_t value  = 0.0;
1604
1605     for (Int_t k = 0; k < loop; k++) {
1606       detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1607       AliTRDCalROC *calROC = object->GetCalROC(detector);
1608       Float_t min          = detobject->GetValue(detector);
1609       Int_t rowMax    = calROC->GetNrows();
1610       Int_t colMax    = calROC->GetNcols();
1611       for (Int_t row = 0; row < rowMax; row++) {
1612         for (Int_t col = 0; col < colMax; col++) {
1613           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1614           calROC->SetValue(col,row,value-min);
1615         } // Col
1616       } // Row
1617     }
1618   }
1619   return object;
1620
1621 }
1622 //_____________________________________________________________________________
1623 TObject *AliTRDCalibraFit::CreatePadObjectPRF(TObjArray *vectorFit)
1624 {
1625   //
1626   // It Creates the AliTRDCalPad object from AliTRDFitInfo
1627   // This object has to be written in the database
1628   //
1629
1630   // Create the DetObject
1631   AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
1632
1633   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1634   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1635   Int_t detector = -1;
1636   Float_t value  = 0.0;
1637
1638   for (Int_t k = 0; k < loop; k++) {
1639     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1640     AliTRDCalROC *calROC = object->GetCalROC(detector);
1641     Int_t rowMax    = calROC->GetNrows();
1642     Int_t colMax    = calROC->GetNcols();
1643     for (Int_t row = 0; row < rowMax; row++) {
1644       for (Int_t col = 0; col < colMax; col++) {
1645         value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1646         calROC->SetValue(col,row,TMath::Abs(value));
1647       } // Col
1648     } // Row
1649   }
1650
1651   return object;
1652
1653 }
1654 //_____________________________________________________________________________
1655 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(TObjArray *vectorFit, const char *name, Double_t &mean)
1656 {
1657   //
1658   // It Creates the AliTRDCalDet object from AliTRDFitInfo
1659   // 0 successful fit 1 not successful fit
1660   // mean is the mean value over the successful fit
1661   // do not use it for t0: no meaning
1662   //
1663
1664   // Create the CalObject
1665   AliTRDCalDet *object = new AliTRDCalDet(name,name);
1666   mean = 0.0;
1667   Int_t count = 0;
1668
1669   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1670   if(loop != 540) {
1671     AliInfo("The Vector Fit is not complete! We initialise all outliers");
1672     for(Int_t k = 0; k < 540; k++){
1673       object->SetValue(k,1.0);
1674     }
1675   }
1676   Int_t detector = -1;
1677   Float_t value  = 0.0;
1678
1679   for (Int_t k = 0; k < loop; k++) {
1680     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1681     value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1682     if(value <= 0) object->SetValue(detector,1.0);
1683     else {
1684       object->SetValue(detector,0.0);
1685       mean += value;
1686       count++;
1687     }
1688   }
1689   if(count > 0) mean /= count;
1690   return object;
1691 }
1692 //_____________________________________________________________________________
1693 TObject *AliTRDCalibraFit::MakeOutliersStatPad(TObjArray *vectorFit, const char *name, Double_t &mean)
1694 {
1695   //
1696   // It Creates the AliTRDCalPad object from AliTRDFitInfo
1697   // 0 not successful fit 1 successful fit
1698   // mean mean value over the successful fit
1699   //
1700
1701   // Create the CalObject
1702   AliTRDCalPad *object = new AliTRDCalPad(name,name);
1703   mean = 0.0;
1704   Int_t count = 0;
1705
1706   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1707   if(loop != 540) {
1708     AliInfo("The Vector Fit is not complete! We initialise all outliers");
1709     for(Int_t k = 0; k < 540; k++){
1710       AliTRDCalROC *calROC = object->GetCalROC(k);
1711       Int_t nchannels = calROC->GetNchannels();
1712       for(Int_t ch = 0; ch < nchannels; ch++){
1713         calROC->SetValue(ch,1.0);
1714       }
1715     }
1716   }
1717   Int_t detector = -1;
1718   Float_t value  = 0.0;
1719
1720   for (Int_t k = 0; k < loop; k++) {
1721     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1722     AliTRDCalROC *calROC = object->GetCalROC(detector);
1723     Int_t nchannels    = calROC->GetNchannels();
1724     for (Int_t ch = 0; ch < nchannels; ch++) {
1725       value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
1726       if(value <= 0) calROC->SetValue(ch,1.0);
1727       else {
1728         calROC->SetValue(ch,0.0);
1729         mean += value;
1730         count++;
1731       }
1732     } // channels
1733   }
1734   if(count > 0) mean /= count;
1735   return object;
1736 }
1737 //_____________________________________________________________________________
1738 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
1739 {
1740   //
1741   // Set FitPH if 1 then each detector will be fitted
1742   //
1743
1744   if (periodeFitPH > 0) {
1745     fFitPHPeriode   = periodeFitPH;
1746   }
1747   else {
1748     AliInfo("periodeFitPH must be higher than 0!");
1749   }
1750
1751 }
1752 //_____________________________________________________________________________
1753 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
1754 {
1755   //
1756   // The fit of the deposited charge distribution begins at
1757   // histo->Mean()/beginFitCharge
1758   // You can here set beginFitCharge
1759   //
1760
1761   if (beginFitCharge > 0) {
1762     fBeginFitCharge = beginFitCharge;
1763   }
1764   else {
1765     AliInfo("beginFitCharge must be strict positif!");
1766   }
1767
1768 }
1769
1770 //_____________________________________________________________________________
1771 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
1772 {
1773   //
1774   // The t0 calculated with the maximum positif slope is shift from t0Shift0
1775   // You can here set t0Shift0
1776   //
1777
1778   if (t0Shift > 0) {
1779     fT0Shift0 = t0Shift;
1780   }
1781   else {
1782     AliInfo("t0Shift0 must be strict positif!");
1783   }
1784
1785 }
1786
1787 //_____________________________________________________________________________
1788 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
1789 {
1790   //
1791   // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
1792   // You can here set t0Shift1
1793   //
1794
1795   if (t0Shift > 0) {
1796     fT0Shift1 = t0Shift;
1797   }
1798   else {
1799     AliInfo("t0Shift must be strict positif!");
1800   }
1801
1802 }
1803
1804 //_____________________________________________________________________________
1805 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
1806 {
1807   //
1808   // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
1809   // You can here set rangeFitPRF
1810   //
1811
1812   if ((rangeFitPRF >    0) &&
1813       (rangeFitPRF <= 1.5)) {
1814     fRangeFitPRF = rangeFitPRF;
1815   }
1816   else {
1817     AliInfo("rangeFitPRF must be between 0 and 1.0");
1818   }
1819
1820 }
1821
1822 //_____________________________________________________________________________
1823 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
1824 {
1825   //
1826   // Minimum entries for fitting
1827   //
1828
1829   if (minEntries >    0) {
1830     fMinEntries = minEntries;
1831   }
1832   else {
1833     AliInfo("fMinEntries must be >= 0.");
1834   }
1835
1836 }
1837
1838 //_____________________________________________________________________________
1839 void AliTRDCalibraFit::SetRebin(Short_t rebin)
1840 {
1841   //
1842   // Rebin with rebin time less bins the Ch histo
1843   // You can set here rebin that should divide the number of bins of CH histo
1844   //
1845
1846   if (rebin > 0) {
1847     fRebin = rebin;
1848     AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
1849   }
1850   else {
1851     AliInfo("You have to choose a positiv value!");
1852   }
1853
1854 }
1855 //_____________________________________________________________________________
1856 Bool_t AliTRDCalibraFit::FillVectorFit()
1857 {
1858   //
1859   // For the Fit functions fill the vector Fit
1860   //
1861
1862   AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
1863
1864   Int_t ntotal = 1;
1865   if (GetChamber(fCountDet) == 2) {
1866     ntotal = 1728;
1867   }
1868   else {
1869     ntotal = 2304;
1870   }
1871
1872   //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
1873   Float_t *coef = new Float_t[ntotal];
1874   for (Int_t i = 0; i < ntotal; i++) {
1875     coef[i] = fCurrentCoefDetector[i];
1876   }
1877
1878   Int_t detector = fCountDet;
1879   // Set
1880   fitInfo->SetCoef(coef);
1881   fitInfo->SetDetector(detector);
1882   fVectorFit.Add((TObject *) fitInfo);
1883
1884   return kTRUE;
1885
1886 }
1887 //_____________________________________________________________________________
1888 Bool_t AliTRDCalibraFit::FillVectorFit2()
1889 {
1890   //
1891   // For the Fit functions fill the vector Fit
1892   //
1893
1894   AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
1895
1896   Int_t ntotal = 1;
1897   if (GetChamber(fCountDet) == 2) {
1898     ntotal = 1728;
1899   }
1900   else {
1901     ntotal = 2304;
1902   }
1903
1904   //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
1905   Float_t *coef = new Float_t[ntotal];
1906   for (Int_t i = 0; i < ntotal; i++) {
1907     coef[i] = fCurrentCoefDetector2[i];
1908   }
1909
1910   Int_t detector = fCountDet;
1911   // Set
1912   fitInfo->SetCoef(coef);
1913   fitInfo->SetDetector(detector);
1914   fVectorFit2.Add((TObject *) fitInfo);
1915
1916   return kTRUE;
1917
1918 }
1919 //____________Functions for initialising the AliTRDCalibraFit in the code_________
1920 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
1921 {
1922   //
1923   // Init the number of expected bins and fDect1[i] fDect2[i]
1924   //
1925
1926   gStyle->SetPalette(1);
1927   gStyle->SetOptStat(1111);
1928   gStyle->SetPadBorderMode(0);
1929   gStyle->SetCanvasColor(10);
1930   gStyle->SetPadLeftMargin(0.13);
1931   gStyle->SetPadRightMargin(0.01);
1932
1933   // Mode groups of pads: the total number of bins!
1934   CalculNumberOfBinsExpected(i);
1935
1936   // Quick verification that we have the good pad calibration mode!
1937   if (fNumberOfBinsExpected != nbins) {
1938     AliInfo("It doesn't correspond to the mode of pad group calibration!");
1939     return kFALSE;
1940   }
1941
1942   // Security for fDebug 3 and 4
1943   if ((fDebugLevel >= 3) &&
1944       ((fDet[0] >  5) ||
1945        (fDet[1] >  4) ||
1946        (fDet[2] > 17))) {
1947     AliInfo("This detector doesn't exit!");
1948     return kFALSE;
1949   }
1950
1951   // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
1952   CalculDect1Dect2(i);
1953
1954
1955   return kTRUE;
1956 }
1957 //____________Functions for initialising the AliTRDCalibraFit in the code_________
1958 Bool_t AliTRDCalibraFit::InitFitCH()
1959 {
1960   //
1961   // Init the fVectorFitCH for normalisation
1962   // Init the histo for debugging
1963   //
1964
1965   gDirectory = gROOT;
1966
1967   fScaleFitFactor = 0.0;
1968   fCurrentCoefDetector   = new Float_t[2304];
1969   for (Int_t k = 0; k < 2304; k++) {
1970     fCurrentCoefDetector[k] = 0.0;
1971   }
1972   fVectorFit.SetName("gainfactorscoefficients");
1973
1974   // fDebug == 0 nothing
1975   // fDebug == 1 and fFitVoir no histo
1976   if (fDebugLevel == 1) {
1977     if(!CheckFitVoir()) return kFALSE;
1978   }
1979   //Get the CalDet object
1980   if(fAccCDB){
1981     AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
1982     if (!cal) {
1983       AliInfo("Could not get calibDB");
1984       return kFALSE;
1985     }
1986     if(fCalDet) delete fCalDet;
1987     fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
1988   }
1989   else{
1990     Float_t devalue = 1.0;
1991     if(fCalDet) delete fCalDet;
1992     fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1993     for(Int_t k = 0; k < 540; k++){
1994       fCalDet->SetValue(k,devalue);
1995     }
1996   }
1997   return kTRUE;
1998
1999 }
2000 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2001 Bool_t AliTRDCalibraFit::InitFitPH()
2002 {
2003   //
2004   // Init the arrays of results
2005   // Init the histos for debugging
2006   //
2007
2008   gDirectory = gROOT;
2009   fVectorFit.SetName("driftvelocitycoefficients");
2010   fVectorFit2.SetName("t0coefficients");
2011
2012   fCurrentCoefDetector   = new Float_t[2304];
2013   for (Int_t k = 0; k < 2304; k++) {
2014     fCurrentCoefDetector[k] = 0.0;
2015   }
2016
2017   fCurrentCoefDetector2   = new Float_t[2304];
2018   for (Int_t k = 0; k < 2304; k++) {
2019     fCurrentCoefDetector2[k] = 0.0;
2020   }
2021
2022   //fDebug == 0 nothing
2023   // fDebug == 1 and fFitVoir no histo
2024   if (fDebugLevel == 1) {
2025     if(!CheckFitVoir()) return kFALSE;
2026   }
2027   //Get the CalDet object
2028   if(fAccCDB){
2029     AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
2030     if (!cal) {
2031       AliInfo("Could not get calibDB");
2032       return kFALSE;
2033     }
2034     if(fCalDet) delete fCalDet;
2035     if(fCalDet2) delete fCalDet2;
2036     fCalDet  = new AliTRDCalDet(*(cal->GetVdriftDet()));
2037     fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2038   }
2039   else{
2040     Float_t devalue  = 1.5;
2041     Float_t devalue2 = 0.0;
2042     if(fCalDet) delete fCalDet;
2043     if(fCalDet2) delete fCalDet2;
2044     fCalDet  = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2045     fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2046     for(Int_t k = 0; k < 540; k++){
2047       fCalDet->SetValue(k,devalue);
2048       fCalDet2->SetValue(k,devalue2);
2049     }
2050   }
2051   return kTRUE;
2052 }
2053 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2054 Bool_t AliTRDCalibraFit::InitFitPRF()
2055 {
2056   //
2057   // Init the calibration mode (Nz, Nrphi), the histograms for
2058   // debugging the fit methods if fDebug > 0,
2059   //
2060
2061   gDirectory = gROOT;
2062   fVectorFit.SetName("prfwidthcoefficients");
2063
2064   fCurrentCoefDetector   = new Float_t[2304];
2065   for (Int_t k = 0; k < 2304; k++) {
2066     fCurrentCoefDetector[k] = 0.0;
2067   }
2068
2069   // fDebug == 0 nothing
2070   // fDebug == 1 and fFitVoir no histo
2071   if (fDebugLevel == 1) {
2072     if(!CheckFitVoir()) return kFALSE;
2073   }
2074   return kTRUE;
2075 }
2076 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2077 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2078 {
2079   //
2080   // Init the fCalDet, fVectorFit fCurrentCoefDetector
2081   //
2082
2083   gDirectory = gROOT;
2084
2085   fCurrentCoefDetector   = new Float_t[2304];
2086   fCurrentCoefDetector2  = new Float_t[2304];
2087   for (Int_t k = 0; k < 2304; k++) {
2088     fCurrentCoefDetector[k]  = 0.0;
2089     fCurrentCoefDetector2[k] = 0.0;
2090   }
2091
2092   //printf("test0\n");
2093
2094   AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
2095   if (!cal) {
2096     AliInfo("Could not get calibDB");
2097     return kFALSE;
2098   }
2099
2100   //Get the CalDet object
2101   if(fAccCDB){
2102     if(fCalDet) delete fCalDet;
2103     if(fCalDet2) delete fCalDet2;
2104     fCalDet  = new AliTRDCalDet(*(cal->GetVdriftDet()));
2105     //printf("test1\n");
2106     fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2107     //printf("test2\n");
2108     for(Int_t k = 0; k < 540; k++){
2109       fCalDet2->SetValue(k,cal->GetOmegaTau(fCalDet->GetValue(k),-fMagneticField));
2110     }
2111     //printf("test3\n");
2112   }
2113   else{
2114     Float_t devalue  = 1.5;
2115     Float_t devalue2 = cal->GetOmegaTau(1.5,-fMagneticField);
2116     if(fCalDet) delete fCalDet;
2117     if(fCalDet2) delete fCalDet2;
2118     //printf("test1\n");
2119     fCalDet  = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2120     fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2121     //printf("test2\n");
2122     for(Int_t k = 0; k < 540; k++){
2123       fCalDet->SetValue(k,devalue);
2124       fCalDet2->SetValue(k,devalue2);
2125     }
2126     //printf("test3\n");
2127   }
2128   return kTRUE;
2129 }
2130
2131 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2132 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2133 {
2134   //
2135   // Init the current detector where we are fCountDet and the
2136   // next fCount for the functions Fit...
2137   //
2138
2139   // Loop on the Xbins of ch!!
2140   fCountDet = -1; // Current detector
2141   fCount    =  0; // To find the next detector
2142
2143   // If fDebug >= 3
2144   if (fDebugLevel >= 3) {
2145     // Set countdet to the detector
2146     fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2147     // Set counter to write at the end of the detector
2148     fCount = fDect2;
2149     // Get the right calib objects
2150     SetCalROC(i);
2151   }
2152   if(fDebugLevel == 1) {
2153     fCountDet = 0;
2154     fCalibraMode->CalculXBins(fCountDet,i);
2155     while(fCalibraMode->GetXbins(i) <=fFitVoir){
2156       fCountDet++;
2157       fCalibraMode->CalculXBins(fCountDet,i);
2158     }
2159     fCount    = fCalibraMode->GetXbins(i);
2160     fCountDet--;
2161     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2162     fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet),i);
2163     fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet)
2164                                        ,(Int_t) GetChamber(fCountDet)
2165                                        ,(Int_t) GetSector(fCountDet),i);
2166   }
2167 }
2168 //_______________________________________________________________________________
2169 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2170 {
2171   //
2172   // Calculate the number of bins expected (calibration groups)
2173   //
2174
2175   fNumberOfBinsExpected = 0;
2176   fCalibraMode->ModePadCalibration(2,i);
2177   fCalibraMode->ModePadFragmentation(0,2,0,i);
2178   fCalibraMode->SetDetChamb2(i);
2179   if (fDebugLevel > 1) {
2180     AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2181   }
2182   fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2183   fCalibraMode->ModePadCalibration(0,i);
2184   fCalibraMode->ModePadFragmentation(0,0,0,i);
2185   fCalibraMode->SetDetChamb0(i);
2186   if (fDebugLevel > 1) {
2187     AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2188   }
2189   fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2190
2191 }
2192 //_______________________________________________________________________________
2193 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2194 {
2195   //
2196   // Calculate the range of fits
2197   //
2198
2199   fDect1 = -1;
2200   fDect2 = -1;
2201   if (fDebugLevel == 1) {
2202     fDect1 = fFitVoir;
2203     fDect2 = fDect1 +1;
2204   }
2205   if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2206     fDect1 = 0;
2207     fDect2 = fNumberOfBinsExpected;
2208   }
2209   if (fDebugLevel >= 3) {
2210     fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2211     fCalibraMode->CalculXBins(fCountDet,i);
2212     fDect1 = fCalibraMode->GetXbins(i);
2213     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2214     fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet),i);
2215     fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet)
2216                                        ,(Int_t) GetChamber(fCountDet)
2217                                        ,(Int_t) GetSector(fCountDet),i);
2218     // Set for the next detector
2219     fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2220   }
2221 }
2222 //_______________________________________________________________________________
2223 Bool_t AliTRDCalibraFit::CheckFitVoir()
2224 {
2225   //
2226   // Check if fFitVoir is in the range
2227   //
2228
2229   if (fFitVoir < fNumberOfBinsExpected) {
2230     AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2231   }
2232   else {
2233     AliInfo("fFitVoir is out of range of the histo!");
2234     return kFALSE;
2235   }
2236   return kTRUE;
2237 }
2238 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2239 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2240 {
2241   //
2242   // See if we are in a new detector and update the
2243   // variables fNfragZ and fNfragRphi if yes
2244   // Will never happen for only one detector (3 and 4)
2245   // Doesn't matter for 2
2246   //
2247   if (fCount == idect) {
2248      // On en est au detector
2249      fCountDet += 1;
2250      // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2251      fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet),i);
2252      fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet)
2253                                 ,(Int_t) GetChamber(fCountDet)
2254                           ,(Int_t) GetSector(fCountDet),i);
2255      // Set for the next detector
2256      fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2257      // calib objects
2258      SetCalROC(i);
2259     }
2260 }
2261 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2262 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2263 {
2264   //
2265   // Reconstruct the min pad row, max pad row, min pad col and
2266   // max pad col of the calibration group for the Fit functions
2267   //
2268   if (fDebugLevel !=  1) {
2269     fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2270   }
2271 }
2272 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2273 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2274 {
2275   //
2276   // For the case where there are not enough entries in the histograms
2277   // of the calibration group, the value present in the choosen database
2278   // will be put. A negativ sign enables to know that a fit was not possible.
2279   //
2280
2281   if (fDebugLevel == 1) {
2282     AliInfo("The element has not enough statistic to be fitted");
2283   }
2284
2285   else {
2286
2287     AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2288                  ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2289
2290     // Calcul the coef from the database choosen
2291     CalculChargeCoefMean(kFALSE);
2292
2293     //chamber 2, not chamber 2
2294     Int_t factor = 0;
2295     if(GetChamber(fCountDet) == 2) factor = 12;
2296     else factor = 16;
2297
2298     // Fill the fCurrentCoefDetector with negative value to say: not fitted
2299     for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2300       for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2301         fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2302       }
2303     }
2304
2305     //Put default value negative
2306     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2307     fCurrentCoefE   = 0.0;
2308
2309     FillFillCH(idect);
2310   }
2311
2312   return kTRUE;
2313 }
2314
2315
2316 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2317 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect)
2318 {
2319   //
2320   // For the case where there are not enough entries in the histograms
2321   // of the calibration group, the value present in the choosen database
2322   // will be put. A negativ sign enables to know that a fit was not possible.
2323   //
2324   if (fDebugLevel == 1) {
2325     AliInfo("The element has not enough statistic to be fitted");
2326   }
2327   else {
2328
2329     AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2330                  ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
2331
2332     CalculVdriftCoefMean();
2333     CalculT0CoefMean();
2334
2335     //chamber 2 and not chamber 2
2336     Int_t factor = 0;
2337     if(GetChamber(fCountDet) == 2) factor = 12;
2338     else factor = 16;
2339
2340
2341     // Fill the fCurrentCoefDetector 2
2342     for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2343       for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2344         fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2345         fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1];
2346       }
2347     }
2348
2349     // Put the default value
2350     fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
2351     fCurrentCoefE    = 0.0;
2352     fCurrentCoef2[0] = fCurrentCoef2[1];
2353     fCurrentCoefE2   = 0.0;
2354
2355     FillFillPH(idect);
2356
2357   }
2358
2359   return kTRUE;
2360
2361 }
2362
2363
2364 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2365 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
2366 {
2367   //
2368   // For the case where there are not enough entries in the histograms
2369   // of the calibration group, the value present in the choosen database
2370   // will be put. A negativ sign enables to know that a fit was not possible.
2371   //
2372
2373   if (fDebugLevel == 1) {
2374     AliInfo("The element has not enough statistic to be fitted");
2375   }
2376   else {
2377
2378     AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2379                  ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
2380
2381     CalculPRFCoefMean();
2382
2383     // chamber 2 and not chamber 2
2384     Int_t factor = 0;
2385     if(GetChamber(fCountDet) == 2) factor = 12;
2386     else factor = 16;
2387
2388
2389     // Fill the fCurrentCoefDetector
2390     for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2391       for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2392         fCurrentCoefDetector[(Int_t)(j*factor+k)] = -fCurrentCoef[1];
2393       }
2394     }
2395
2396     // Put the default value
2397     fCurrentCoef[0] = -fCurrentCoef[1];
2398     fCurrentCoefE   = 0.0;
2399
2400     FillFillPRF(idect);
2401   }
2402
2403   return kTRUE;
2404
2405 }
2406 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2407 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
2408 {
2409   //
2410   // For the case where there are not enough entries in the histograms
2411   // of the calibration group, the value present in the choosen database
2412   // will be put. A negativ sign enables to know that a fit was not possible.
2413   //
2414
2415   // Calcul the coef from the database choosen
2416   CalculVdriftLorentzCoef();
2417
2418   Int_t factor = 0;
2419   if(GetChamber(fCountDet) == 2) factor = 1728;
2420   else factor = 2304;
2421
2422
2423   // Fill the fCurrentCoefDetector
2424   for (Int_t k = 0; k < factor; k++) {
2425     fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
2426     // should be negative
2427     fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
2428   }
2429
2430
2431   //Put default opposite sign
2432   fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
2433   fCurrentCoefE    = 0.0;
2434   fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
2435   fCurrentCoefE2 = 0.0;
2436
2437   FillFillLinearFitter();
2438
2439   return kTRUE;
2440 }
2441
2442 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2443 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
2444 {
2445   //
2446   // Fill the coefficients found with the fits or other
2447   // methods from the Fit functions
2448   //
2449
2450   if (fDebugLevel != 1) {
2451
2452     Int_t factor = 0;
2453     if(GetChamber(fCountDet) == 2) factor = 12;
2454     else factor = 16;
2455
2456     for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2457       for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2458         fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2459       }
2460     }
2461
2462     FillFillCH(idect);
2463
2464   }
2465
2466   return kTRUE;
2467
2468 }
2469 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2470 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect)
2471 {
2472   //
2473   // Fill the coefficients found with the fits or other
2474   // methods from the Fit functions
2475   //
2476
2477   if (fDebugLevel != 1) {
2478
2479     Int_t factor = 0;
2480     if(GetChamber(fCountDet) == 2) factor = 12;
2481     else factor = 16;
2482
2483     for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2484       for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2485         fCurrentCoefDetector[(Int_t)(j*factor+k)]  = fCurrentCoef[0];
2486         fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
2487       }
2488     }
2489     FillFillPH(idect);
2490   }
2491   return kTRUE;
2492 }
2493 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2494 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
2495 {
2496   //
2497   // Fill the coefficients found with the fits or other
2498   // methods from the Fit functions
2499   //
2500
2501   if (fDebugLevel != 1) {
2502
2503     Int_t factor = 0;
2504     if(GetChamber(fCountDet) == 2) factor = 12;
2505     else factor = 16;
2506
2507     // Pointer to the branch
2508     for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2509       for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2510         fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2511       }
2512     }
2513     FillFillPRF(idect);
2514   }
2515
2516   return kTRUE;
2517
2518 }
2519 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2520 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
2521 {
2522   //
2523   // Fill the coefficients found with the fits or other
2524   // methods from the Fit functions
2525   //
2526
2527   Int_t factor = 0;
2528   if(GetChamber(fCountDet) == 2) factor = 1728;
2529   else factor = 2304;
2530
2531   // Pointer to the branch
2532   for (Int_t k = 0; k < factor; k++) {
2533     fCurrentCoefDetector[k]  = fCurrentCoef[0];
2534     fCurrentCoefDetector2[k] = fCurrentCoef2[0];
2535   }
2536
2537   FillFillLinearFitter();
2538
2539   return kTRUE;
2540
2541 }
2542 //________________________________________________________________________________
2543 void AliTRDCalibraFit::FillFillCH(Int_t idect)
2544 {
2545   //
2546   // DebugStream and fVectorFit
2547   //
2548
2549   // End of one detector
2550   if ((idect == (fCount-1))) {
2551     FillVectorFit();
2552     // Reset
2553     for (Int_t k = 0; k < 2304; k++) {
2554       fCurrentCoefDetector[k] = 0.0;
2555     }
2556   }
2557
2558   if(fDebugLevel > 1){
2559
2560     if ( !fDebugStreamer ) {
2561       //debug stream
2562       TDirectory *backup = gDirectory;
2563       fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2564       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2565     }
2566
2567     Int_t   detector   = fCountDet;
2568     Int_t   caligroup  = idect;
2569     Short_t rowmin     = fCalibraMode->GetRowMin(0);
2570     Short_t rowmax     = fCalibraMode->GetRowMax(0);
2571     Short_t colmin     = fCalibraMode->GetColMin(0);
2572     Short_t colmax     = fCalibraMode->GetColMax(0);
2573     Float_t gf         = fCurrentCoef[0];
2574     Float_t gfs        = fCurrentCoef[1];
2575     Float_t gfE        = fCurrentCoefE;
2576
2577     (*fDebugStreamer) << "FillFillCH" <<
2578       "detector=" << detector <<
2579       "caligroup=" << caligroup <<
2580       "rowmin=" << rowmin <<
2581       "rowmax=" << rowmax <<
2582       "colmin=" << colmin <<
2583       "colmax=" << colmax <<
2584       "gf=" << gf <<
2585       "gfs=" << gfs <<
2586       "gfE=" << gfE <<
2587       "\n";
2588
2589   }
2590 }
2591 //________________________________________________________________________________
2592 void AliTRDCalibraFit::FillFillPH(Int_t idect)
2593 {
2594   //
2595   // DebugStream and fVectorFit and fVectorFit2
2596   //
2597
2598   // End of one detector
2599     if ((idect == (fCount-1))) {
2600       FillVectorFit();
2601       FillVectorFit2();
2602       // Reset
2603       for (Int_t k = 0; k < 2304; k++) {
2604         fCurrentCoefDetector[k] = 0.0;
2605         fCurrentCoefDetector2[k] = 0.0;
2606       }
2607     }
2608
2609     if(fDebugLevel > 1){
2610
2611       if ( !fDebugStreamer ) {
2612         //debug stream
2613         TDirectory *backup = gDirectory;
2614         fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2615         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2616       }
2617
2618
2619       Int_t   detector     = fCountDet;
2620       Int_t   caligroup    = idect;
2621       Short_t rowmin       = fCalibraMode->GetRowMin(1);
2622       Short_t rowmax       = fCalibraMode->GetRowMax(1);
2623       Short_t colmin       = fCalibraMode->GetColMin(1);
2624       Short_t colmax       = fCalibraMode->GetColMax(1);
2625       Float_t vf           = fCurrentCoef[0];
2626       Float_t vs           = fCurrentCoef[1];
2627       Float_t vfE          = fCurrentCoefE;
2628       Float_t t0f          = fCurrentCoef2[0];
2629       Float_t t0s          = fCurrentCoef2[1];
2630       Float_t t0E          = fCurrentCoefE2;
2631
2632
2633
2634       (* fDebugStreamer) << "FillFillPH"<<
2635         "detector="<<detector<<
2636         "caligroup="<<caligroup<<
2637         "rowmin="<<rowmin<<
2638         "rowmax="<<rowmax<<
2639         "colmin="<<colmin<<
2640         "colmax="<<colmax<<
2641         "vf="<<vf<<
2642         "vs="<<vs<<
2643         "vfE="<<vfE<<
2644         "t0f="<<t0f<<
2645         "t0s="<<t0s<<
2646         "t0E="<<t0E<<
2647         "\n";
2648     }
2649
2650 }
2651 //________________________________________________________________________________
2652 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
2653 {
2654   //
2655   // DebugStream and fVectorFit
2656   //
2657
2658     // End of one detector
2659     if ((idect == (fCount-1))) {
2660       FillVectorFit();
2661       // Reset
2662       for (Int_t k = 0; k < 2304; k++) {
2663         fCurrentCoefDetector[k] = 0.0;
2664       }
2665     }
2666
2667
2668     if(fDebugLevel > 1){
2669
2670       if ( !fDebugStreamer ) {
2671         //debug stream
2672         TDirectory *backup = gDirectory;
2673         fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2674         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2675       }
2676
2677       Int_t   detector     = fCountDet;
2678       Int_t   plane        = GetPlane(fCountDet);
2679       Int_t   caligroup    = idect;
2680       Short_t rowmin       = fCalibraMode->GetRowMin(2);
2681       Short_t rowmax       = fCalibraMode->GetRowMax(2);
2682       Short_t colmin       = fCalibraMode->GetColMin(2);
2683       Short_t colmax       = fCalibraMode->GetColMax(2);
2684       Float_t widf         = fCurrentCoef[0];
2685       Float_t wids         = fCurrentCoef[1];
2686       Float_t widfE        = fCurrentCoefE;
2687
2688       (* fDebugStreamer) << "FillFillPRF"<<
2689         "detector="<<detector<<
2690         "plane="<<plane<<
2691         "caligroup="<<caligroup<<
2692         "rowmin="<<rowmin<<
2693         "rowmax="<<rowmax<<
2694         "colmin="<<colmin<<
2695         "colmax="<<colmax<<
2696         "widf="<<widf<<
2697         "wids="<<wids<<
2698         "widfE="<<widfE<<
2699         "\n";
2700     }
2701
2702 }
2703 //________________________________________________________________________________
2704 void AliTRDCalibraFit::FillFillLinearFitter()
2705 {
2706   //
2707   // DebugStream and fVectorFit
2708   //
2709
2710   // End of one detector
2711   FillVectorFit();
2712   FillVectorFit2();
2713
2714
2715   // Reset
2716   for (Int_t k = 0; k < 2304; k++) {
2717   fCurrentCoefDetector[k]  = 0.0;
2718   fCurrentCoefDetector2[k] = 0.0;
2719   }
2720
2721
2722   if(fDebugLevel > 1){
2723
2724     if ( !fDebugStreamer ) {
2725       //debug stream
2726       TDirectory *backup = gDirectory;
2727       fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2728       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2729     }
2730
2731     //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
2732     AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(fCountDet),GetChamber(fCountDet));
2733     Float_t rowmd            = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
2734     Float_t r                = AliTRDgeometry::GetTime0(GetPlane(fCountDet));
2735     Float_t tiltangle        = padplane->GetTiltingAngle();
2736     Int_t   detector         = fCountDet;
2737     Int_t   chamber          = GetChamber(fCountDet);
2738     Int_t   plane            = GetChamber(fCountDet);
2739     Float_t vf               = fCurrentCoef[0];
2740     Float_t vs               = fCurrentCoef[1];
2741     Float_t vfE              = fCurrentCoefE;
2742     Float_t lorentzangler    = fCurrentCoef2[0];
2743     Float_t elorentzangler   = fCurrentCoefE2;
2744     Float_t lorentzangles    = fCurrentCoef2[1];
2745
2746     (* fDebugStreamer) << "FillFillLinearFitter"<<
2747       "detector="<<detector<<
2748       "chamber="<<chamber<<
2749       "plane="<<plane<<
2750       "rowmd="<<rowmd<<
2751       "r="<<r<<
2752       "tiltangle="<<tiltangle<<
2753       "vf="<<vf<<
2754       "vs="<<vs<<
2755       "vfE="<<vfE<<
2756       "lorentzangler="<<lorentzangler<<
2757       "Elorentzangler="<<elorentzangler<<
2758       "lorentzangles="<<lorentzangles<<
2759       "\n";
2760   }
2761
2762 }
2763 //
2764 //____________Calcul Coef Mean_________________________________________________
2765 //
2766 //_____________________________________________________________________________
2767 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
2768 {
2769   //
2770   // For the detector Dect calcul the mean time 0
2771   // for the calibration group idect from the choosen database
2772   //
2773
2774   fCurrentCoef2[1] = 0.0;
2775   if(fDebugLevel != 1){
2776     if((fCalibraMode->GetNz(1) > 0) ||
2777        (fCalibraMode->GetNrphi(1) > 0)) {
2778       for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
2779         for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
2780           fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
2781         }
2782       }
2783       fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
2784     }
2785     else {
2786       if(!fAccCDB){
2787         fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
2788       }
2789       else{
2790         for(Int_t row = 0; row < fGeo->GetRowMax(GetPlane(fCountDet),GetChamber(fCountDet),GetSector(fCountDet)); row++){
2791           for(Int_t col = 0; col < fGeo->GetColMax(GetPlane(fCountDet)); col++){
2792             fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
2793           }
2794         }
2795         fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetPlane(fCountDet),GetChamber(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetPlane(fCountDet))));
2796       }
2797     }
2798   }
2799   return kTRUE;
2800 }
2801
2802 //_____________________________________________________________________________
2803 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
2804 {
2805   //
2806   // For the detector Dect calcul the mean gain factor
2807   // for the calibration group idect from the choosen database
2808   //
2809
2810   fCurrentCoef[1] = 0.0;
2811   if(fDebugLevel != 1){
2812     if ((fCalibraMode->GetNz(0)    > 0) ||
2813         (fCalibraMode->GetNrphi(0) > 0)) {
2814       for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
2815         for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
2816           fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2817           if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2818         }
2819       }
2820       fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
2821     }
2822     else {
2823       //Per detectors
2824       fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
2825       if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
2826     }
2827   }
2828   return kTRUE;
2829 }
2830 //_____________________________________________________________________________
2831 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
2832 {
2833   //
2834   // For the detector Dect calcul the mean sigma of pad response
2835   // function for the calibration group idect from the choosen database
2836   //
2837
2838   fCurrentCoef[1] = 0.0;
2839   if(fDebugLevel != 1){
2840     for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
2841       for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
2842         fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
2843       }
2844     }
2845     fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
2846   }
2847   return kTRUE;
2848 }
2849 //_____________________________________________________________________________
2850 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
2851 {
2852   //
2853   // For the detector dect calcul the mean drift velocity for the
2854   // calibration group idect from the choosen database
2855   //
2856
2857   fCurrentCoef[1] = 0.0;
2858   if(fDebugLevel != 1){
2859     if ((fCalibraMode->GetNz(1)    > 0) ||
2860         (fCalibraMode->GetNrphi(1) > 0)) {
2861       for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
2862         for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
2863           fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2864         }
2865       }
2866       fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
2867     }
2868     else {
2869       //per detectors
2870       fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
2871     }
2872   }
2873   return kTRUE;
2874 }
2875 //_____________________________________________________________________________
2876 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
2877 {
2878   //
2879   // For the detector fCountDet, mean drift velocity and tan lorentzangle
2880   //
2881
2882   fCurrentCoef[1]  = fCalDet->GetValue(fCountDet);
2883   fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
2884
2885   return kTRUE;
2886 }
2887 //_____________________________________________________________________________
2888 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t plane) const
2889 {
2890   //
2891   // Default width of the PRF if there is no database as reference
2892   //
2893   switch(plane)
2894     {
2895       // default database
2896       //case 0:  return 0.515;
2897       //case 1:  return 0.502;
2898       //case 2:  return 0.491;
2899       //case 3:  return 0.481;
2900       //case 4:  return 0.471;
2901       //case 5:  return 0.463;
2902       //default: return 0.0;
2903
2904       // fit database
2905     case 0:  return 0.538429;
2906     case 1:  return 0.524302;
2907     case 2:  return 0.511591;
2908     case 3:  return 0.500140;
2909     case 4:  return 0.489821;
2910     case 5:  return 0.480524;
2911     default: return 0.0;
2912   }
2913 }
2914 //________________________________________________________________________________
2915 void AliTRDCalibraFit::SetCalROC(Int_t i)
2916 {
2917   //
2918   // Set the calib object for fCountDet
2919   //
2920
2921   Float_t value = 0.0;
2922
2923   //Get the CalDet object
2924   if(fAccCDB){
2925     AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
2926     if (!cal) {
2927       AliInfo("Could not get calibDB");
2928       return;
2929     }
2930     switch (i)
2931       {
2932       case 0:
2933         if(fCalROC) delete fCalROC;
2934         fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
2935         break;
2936       case 1:
2937         if(fCalROC)  delete fCalROC;
2938         if(fCalROC2) delete fCalROC2;
2939         fCalROC  = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
2940         fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
2941         break;
2942       case 2:
2943         if(fCalROC) delete fCalROC;
2944         fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet))); break;
2945       default: return;
2946       }
2947   }
2948   else{
2949     switch (i)
2950       {
2951       case 0:
2952         if(fCalROC) delete fCalROC;
2953         fCalROC = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2954         for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2955           fCalROC->SetValue(k,1.0);
2956         }
2957         break;
2958       case 1:
2959         if(fCalROC)  delete fCalROC;
2960         if(fCalROC2) delete fCalROC2;
2961         fCalROC  = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2962         fCalROC2 = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2963         for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2964           fCalROC->SetValue(k,1.0);
2965           fCalROC2->SetValue(k,0.0);
2966         }
2967         break;
2968       case 2:
2969         if(fCalROC) delete fCalROC;
2970         value = GetPRFDefault(GetPlane(fCountDet));
2971         fCalROC = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2972         for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2973           fCalROC->SetValue(k,value);
2974         }
2975         break;
2976       default: return;
2977       }
2978   }
2979
2980 }
2981 //____________Fit Methods______________________________________________________
2982
2983 //_____________________________________________________________________________
2984 void AliTRDCalibraFit::FitPente(TH1* projPH)
2985 {
2986   //
2987   // Slope methode for the drift velocity
2988   //
2989
2990   // Constants
2991   const Float_t kDrWidth = AliTRDgeometry::DrThick();
2992   Int_t binmax           = 0;
2993   Int_t binmin           = 0;
2994   fPhd[0]                = 0.0;
2995   fPhd[1]                = 0.0;
2996   fPhd[2]                = 0.0;
2997   Int_t ju               = 0;
2998   fCurrentCoefE          = 0.0;
2999   fCurrentCoefE2         = 0.0;
3000   fCurrentCoef[0]        = 0.0;
3001   fCurrentCoef2[0]       = 0.0;
3002   TLine *line            = new TLine();
3003
3004   // Some variables
3005   TAxis   *xpph    = projPH->GetXaxis();
3006   Int_t    nbins   = xpph->GetNbins();
3007   Double_t lowedge = xpph->GetBinLowEdge(1);
3008   Double_t upedge  = xpph->GetBinUpEdge(xpph->GetNbins());
3009   Double_t widbins = (upedge - lowedge) / nbins;
3010   Double_t limit   = upedge + 0.5 * widbins;
3011   Bool_t put = kTRUE;
3012
3013   // Beginning of the signal
3014   TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
3015   for (Int_t k = 1; k <  projPH->GetNbinsX(); k++) {
3016     pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3017   }
3018   binmax = (Int_t) pentea->GetMaximumBin();
3019   if (binmax <= 1) {
3020     binmax = 2;
3021     AliInfo("Put the binmax from 1 to 2 to enable the fit");
3022   }
3023   if (binmax >= nbins) {
3024     binmax = nbins-1;
3025     put = kFALSE;
3026     AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3027   }
3028   pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
3029   Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
3030   Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
3031   Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
3032   Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
3033   if (l3P2am != 0) {
3034     fPhd[0] = -(l3P1am / (2 * l3P2am));
3035   }
3036   if(!fTakeTheMaxPH){
3037     if((l3P1am != 0.0) && (l3P2am != 0.0)){
3038       fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
3039     }
3040   }
3041   // Amplification region
3042   binmax = 0;
3043   ju     = 0;
3044   for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3045     if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3046       binmax = kbin;
3047       ju     = 1;
3048     }
3049   }
3050   if (binmax <= 1) {
3051     binmax = 2;
3052     AliInfo("Put the binmax from 1 to 2 to enable the fit");
3053   }
3054   if (binmax >= nbins) {
3055     binmax = nbins-1;
3056     put = kFALSE;
3057     AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3058   }
3059   projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
3060   Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
3061   Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
3062   Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
3063   Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
3064   if (l3P2amf != 0) {
3065     fPhd[1] = -(l3P1amf / (2 * l3P2amf));
3066   }
3067   if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
3068     fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
3069   }
3070   if(fTakeTheMaxPH){
3071     fCurrentCoefE2 = fCurrentCoefE;
3072   }
3073   // Drift region
3074   TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
3075   for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k <  projPH->GetNbinsX(); k++) {
3076     pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3077   }
3078   binmin = 0;
3079   if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
3080   if (binmin <= 1) {
3081     binmin = 2;
3082     AliInfo("Put the binmax from 1 to 2 to enable the fit");
3083   }
3084   if (binmin >= nbins) {
3085     binmin = nbins-1;
3086     put = kFALSE;
3087     AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3088   }
3089   pente->Fit("pol2"
3090             ,"0MR"
3091             ,""
3092             ,TMath::Max(pente->GetBinCenter(binmin-1),             0.0)
3093             ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
3094   Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
3095   Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
3096   Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
3097   Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
3098   if (l3P2dr != 0) {
3099     fPhd[2] = -(l3P1dr / (2 * l3P2dr));
3100   }
3101   if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
3102     fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
3103   }
3104   Float_t fPhdt0  = 0.0;
3105   Float_t t0Shift = 0.0;
3106   if(fTakeTheMaxPH) {
3107     fPhdt0 = fPhd[1];
3108     t0Shift = fT0Shift1;
3109   }
3110   else {
3111     fPhdt0 = fPhd[0];
3112     t0Shift = fT0Shift0;
3113   }
3114
3115   if ((fPhd[2] > fPhd[0]) &&
3116       (fPhd[2] > fPhd[1]) &&
3117       (fPhd[1] > fPhd[0]) &&
3118       (put)) {
3119     fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
3120     fNumberFitSuccess++;
3121
3122     if (fPhdt0 >= 0.0) {
3123       fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
3124       if (fCurrentCoef2[0] < -1.0) {
3125         fCurrentCoef2[0] = fCurrentCoef2[1];
3126       }
3127     }
3128     else {
3129       fCurrentCoef2[0] = fCurrentCoef2[1];
3130     }
3131
3132   }
3133   else {
3134     fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
3135     fCurrentCoef2[0] = fCurrentCoef2[1];
3136   }
3137
3138   if (fDebugLevel == 1) {
3139     TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
3140     cpentei->cd();
3141     projPH->Draw();
3142     line->SetLineColor(2);
3143     line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
3144     line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
3145     line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
3146     AliInfo(Form("fPhd[0] (beginning of the signal): %f"                  ,(Float_t) fPhd[0]));
3147     AliInfo(Form("fPhd[1] (end of the amplification region): %f"          ,(Float_t) fPhd[1]));
3148     AliInfo(Form("fPhd[2] (end of the drift region): %f"                  ,(Float_t) fPhd[2]));
3149     AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
3150     TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
3151     cpentei2->cd();
3152     pentea->Draw();
3153     TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
3154     cpentei3->cd();
3155     pente->Draw();
3156   }
3157   else {
3158     delete pentea;
3159     delete pente;
3160   }
3161 }
3162 //_____________________________________________________________________________
3163 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
3164 {
3165   //
3166   // Slope methode but with polynomes de Lagrange
3167   //
3168
3169   // Constants
3170   const Float_t kDrWidth = AliTRDgeometry::DrThick();
3171   Int_t binmax      = 0;
3172   Int_t binmin      = 0;
3173   Double_t    *x    = new Double_t[5];
3174   Double_t    *y    = new Double_t[5];
3175   x[0]              = 0.0;
3176   x[1]              = 0.0;
3177   x[2]              = 0.0;
3178   x[3]              = 0.0;
3179   x[4]              = 0.0;
3180   y[0]              = 0.0;
3181   y[1]              = 0.0;
3182   y[2]              = 0.0;
3183   y[3]              = 0.0;
3184   y[4]              = 0.0;
3185   fPhd[0]           = 0.0;
3186   fPhd[1]           = 0.0;
3187   fPhd[2]           = 0.0;
3188   Int_t ju          = 0;
3189   fCurrentCoefE     = 0.0;
3190   fCurrentCoefE2    = 1.0;
3191   fCurrentCoef[0]   = 0.0;
3192   fCurrentCoef2[0]  = 0.0;
3193   TLine *line = new TLine();
3194   TF1 * polynome = 0x0;
3195   TF1 * polynomea = 0x0;
3196   TF1 * polynomeb = 0x0;
3197   Double_t *c = 0x0;
3198
3199   // Some variables
3200   TAxis   *xpph    = projPH->GetXaxis();
3201   Int_t    nbins   = xpph->GetNbins();
3202   Double_t lowedge = xpph->GetBinLowEdge(1);
3203   Double_t upedge  = xpph->GetBinUpEdge(xpph->GetNbins());
3204   Double_t widbins = (upedge - lowedge) / nbins;
3205   Double_t limit   = upedge + 0.5 * widbins;
3206
3207
3208   Bool_t put = kTRUE;
3209
3210   // Beginning of the signal
3211   TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
3212   for (Int_t k = 1; k <  projPH->GetNbinsX(); k++) {
3213     pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3214   }
3215
3216   binmax = (Int_t) pentea->GetMaximumBin();
3217
3218   Double_t minnn = 0.0;
3219   Double_t maxxx = 0.0;
3220
3221   Int_t kase = nbins-binmax;
3222
3223   switch(kase)
3224     {
3225     case 0:
3226       put = kFALSE;
3227       break;
3228     case 1:
3229       minnn = pentea->GetBinCenter(binmax-2);
3230       maxxx = pentea->GetBinCenter(binmax);
3231       x[0] = pentea->GetBinCenter(binmax-2);
3232       x[1] = pentea->GetBinCenter(binmax-1);
3233       x[2] = pentea->GetBinCenter(binmax);
3234       y[0] = pentea->GetBinContent(binmax-2);
3235       y[1] = pentea->GetBinContent(binmax-1);
3236       y[2] = pentea->GetBinContent(binmax);
3237       c = CalculPolynomeLagrange2(x,y);
3238       AliInfo("At the limit for beginning!");
3239       break;
3240     case 2:
3241       minnn = pentea->GetBinCenter(binmax-2);
3242       maxxx = pentea->GetBinCenter(binmax+1);
3243       x[0] = pentea->GetBinCenter(binmax-2);
3244       x[1] = pentea->GetBinCenter(binmax-1);
3245       x[2] = pentea->GetBinCenter(binmax);
3246       x[3] = pentea->GetBinCenter(binmax+1);
3247       y[0] = pentea->GetBinContent(binmax-2);
3248       y[1] = pentea->GetBinContent(binmax-1);
3249       y[2] = pentea->GetBinContent(binmax);
3250       y[3] = pentea->GetBinContent(binmax+1);
3251       c = CalculPolynomeLagrange3(x,y);
3252       break;
3253     default:
3254       switch(binmax){
3255       case 0:
3256         put = kFALSE;
3257         break;
3258       case 1:
3259         minnn = pentea->GetBinCenter(binmax);
3260         maxxx = pentea->GetBinCenter(binmax+2);
3261         x[0] = pentea->GetBinCenter(binmax);
3262         x[1] = pentea->GetBinCenter(binmax+1);
3263         x[2] = pentea->GetBinCenter(binmax+2);
3264         y[0] = pentea->GetBinContent(binmax);
3265         y[1] = pentea->GetBinContent(binmax+1);
3266         y[2] = pentea->GetBinContent(binmax+2);
3267         c = CalculPolynomeLagrange2(x,y);
3268         break;
3269       case 2:
3270         minnn = pentea->GetBinCenter(binmax-1);
3271         maxxx = pentea->GetBinCenter(binmax+2);
3272         x[0] = pentea->GetBinCenter(binmax-1);
3273         x[1] = pentea->GetBinCenter(binmax);
3274         x[2] = pentea->GetBinCenter(binmax+1);
3275         x[3] = pentea->GetBinCenter(binmax+2);
3276         y[0] = pentea->GetBinContent(binmax-1);
3277         y[1] = pentea->GetBinContent(binmax);
3278         y[2] = pentea->GetBinContent(binmax+1);
3279         y[3] = pentea->GetBinContent(binmax+2);
3280         c = CalculPolynomeLagrange3(x,y);
3281         break;
3282       default:
3283         minnn = pentea->GetBinCenter(binmax-2);
3284         maxxx = pentea->GetBinCenter(binmax+2);
3285         x[0] = pentea->GetBinCenter(binmax-2);
3286         x[1] = pentea->GetBinCenter(binmax-1);
3287         x[2] = pentea->GetBinCenter(binmax);
3288         x[3] = pentea->GetBinCenter(binmax+1);
3289         x[4] = pentea->GetBinCenter(binmax+2);
3290         y[0] = pentea->GetBinContent(binmax-2);
3291         y[1] = pentea->GetBinContent(binmax-1);
3292         y[2] = pentea->GetBinContent(binmax);
3293         y[3] = pentea->GetBinContent(binmax+1);
3294         y[4] = pentea->GetBinContent(binmax+2);
3295         c = CalculPolynomeLagrange4(x,y);
3296         break;
3297       }
3298       break;
3299     }
3300
3301
3302   if(put) {
3303     polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
3304     polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3305
3306     Double_t step = (maxxx-minnn)/10000;
3307     Double_t l = minnn;
3308     Double_t maxvalue = 0.0;
3309     Double_t placemaximum = minnn;
3310     for(Int_t o = 0; o < 10000; o++){
3311       if(o == 0) maxvalue = polynomeb->Eval(l);
3312       if(maxvalue < (polynomeb->Eval(l))){
3313         maxvalue = polynomeb->Eval(l);
3314         placemaximum = l;
3315       }
3316       l += step;
3317     }
3318     fPhd[0] = placemaximum;
3319   }
3320
3321   // Amplification region
3322   binmax = 0;
3323   ju     = 0;
3324   for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3325     if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3326       binmax = kbin;
3327       ju     = 1;
3328     }
3329   }
3330
3331   Double_t minn = 0.0;
3332   Double_t maxx = 0.0;
3333   x[0] = 0.0;
3334   x[1] = 0.0;
3335   x[2] = 0.0;
3336   x[3] = 0.0;
3337   x[4] = 0.0;
3338   y[0] = 0.0;
3339   y[1] = 0.0;
3340   y[2] = 0.0;
3341   y[3] = 0.0;
3342   y[4] = 0.0;
3343
3344   Int_t    kase1 = nbins - binmax;
3345
3346   //Determination of minn and maxx
3347   //case binmax = nbins
3348   //pol2
3349   switch(kase1)
3350     {
3351     case 0:
3352       minn = projPH->GetBinCenter(binmax-2);
3353       maxx = projPH->GetBinCenter(binmax);
3354       x[0] = projPH->GetBinCenter(binmax-2);
3355       x[1] = projPH->GetBinCenter(binmax-1);
3356       x[2] = projPH->GetBinCenter(binmax);
3357       y[0] = projPH->GetBinContent(binmax-2);
3358       y[1] = projPH->GetBinContent(binmax-1);
3359       y[2] = projPH->GetBinContent(binmax);
3360       c = CalculPolynomeLagrange2(x,y);
3361       //AliInfo("At the limit for the drift!");
3362       break;
3363     case 1:
3364       minn = projPH->GetBinCenter(binmax-2);
3365       maxx = projPH->GetBinCenter(binmax+1);
3366       x[0] = projPH->GetBinCenter(binmax-2);
3367       x[1] = projPH->GetBinCenter(binmax-1);
3368       x[2] = projPH->GetBinCenter(binmax);
3369       x[3] = projPH->GetBinCenter(binmax+1);
3370       y[0] = projPH->GetBinContent(binmax-2);
3371       y[1] = projPH->GetBinContent(binmax-1);
3372       y[2] = projPH->GetBinContent(binmax);
3373       y[3] = projPH->GetBinContent(binmax+1);
3374       c = CalculPolynomeLagrange3(x,y);
3375       break;
3376     default:
3377       switch(binmax)
3378         {
3379         case 0:
3380           put = kFALSE;
3381           break;
3382         case 1:
3383           minn = projPH->GetBinCenter(binmax);
3384           maxx = projPH->GetBinCenter(binmax+2);
3385           x[0] = projPH->GetBinCenter(binmax);
3386           x[1] = projPH->GetBinCenter(binmax+1);
3387           x[2] = projPH->GetBinCenter(binmax+2);
3388           y[0] = projPH->GetBinContent(binmax);
3389           y[1] = projPH->GetBinContent(binmax+1);
3390           y[2] = projPH->GetBinContent(binmax+2);
3391           c = CalculPolynomeLagrange2(x,y);
3392           break;
3393         case 2:
3394           minn = projPH->GetBinCenter(binmax-1);
3395           maxx = projPH->GetBinCenter(binmax+2);
3396           x[0] = projPH->GetBinCenter(binmax-1);
3397           x[1] = projPH->GetBinCenter(binmax);
3398           x[2] = projPH->GetBinCenter(binmax+1);
3399           x[3] = projPH->GetBinCenter(binmax+2);
3400           y[0] = projPH->GetBinContent(binmax-1);
3401           y[1] = projPH->GetBinContent(binmax);
3402           y[2] = projPH->GetBinContent(binmax+1);
3403           y[3] = projPH->GetBinContent(binmax+2);
3404           c = CalculPolynomeLagrange3(x,y);
3405           break;
3406         default:
3407           minn = projPH->GetBinCenter(binmax-2);
3408           maxx = projPH->GetBinCenter(binmax+2);
3409           x[0] = projPH->GetBinCenter(binmax-2);
3410           x[1] = projPH->GetBinCenter(binmax-1);
3411           x[2] = projPH->GetBinCenter(binmax);
3412           x[3] = projPH->GetBinCenter(binmax+1);
3413           x[4] = projPH->GetBinCenter(binmax+2);
3414           y[0] = projPH->GetBinContent(binmax-2);
3415           y[1] = projPH->GetBinContent(binmax-1);
3416           y[2] = projPH->GetBinContent(binmax);
3417           y[3] = projPH->GetBinContent(binmax+1);
3418           y[4] = projPH->GetBinContent(binmax+2);
3419           c = CalculPolynomeLagrange4(x,y);
3420           break;
3421         }
3422       break;
3423     }
3424
3425   if(put) {
3426     polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
3427     polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3428
3429     Double_t step = (maxx-minn)/1000;
3430     Double_t l = minn;
3431     Double_t maxvalue = 0.0;
3432     Double_t placemaximum = minn;
3433     for(Int_t o = 0; o < 1000; o++){
3434       if(o == 0) maxvalue = polynomea->Eval(l);
3435       if(maxvalue < (polynomea->Eval(l))){
3436         maxvalue = polynomea->Eval(l);
3437         placemaximum = l;
3438       }
3439       l += step;
3440     }
3441     fPhd[1] = placemaximum;
3442   }
3443
3444   // Drift region
3445   TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
3446   for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k <  projPH->GetNbinsX(); k++) {
3447     pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3448   }
3449   binmin = 0;
3450   if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
3451
3452   //should not happen
3453   if (binmin <= 1) {
3454     binmin = 2;
3455     put = 1;
3456     AliInfo("Put the binmax from 1 to 2 to enable the fit");
3457   }
3458
3459   //check
3460   if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) put = kFALSE;
3461   if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) put = kFALSE;
3462
3463   x[0] = 0.0;
3464   x[1] = 0.0;
3465   x[2] = 0.0;
3466   x[3] = 0.0;
3467   x[4] = 0.0;
3468   y[0] = 0.0;
3469   y[1] = 0.0;
3470   y[2] = 0.0;
3471   y[3] = 0.0;
3472   y[4] = 0.0;
3473   Double_t min = 0.0;
3474   Double_t max = 0.0;
3475   Bool_t case1 = kFALSE;
3476   Bool_t case2 = kFALSE;
3477   Bool_t case4 = kFALSE;
3478
3479   //Determination of min and max
3480   //case binmin <= nbins-3
3481   //pol4 case 3
3482   if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
3483     min = pente->GetBinCenter(binmin-2);
3484     max = pente->GetBinCenter(binmin+2);
3485     x[0] = pente->GetBinCenter(binmin-2);
3486     x[1] = pente->GetBinCenter(binmin-1);
3487     x[2] = pente->GetBinCenter(binmin);
3488     x[3] = pente->GetBinCenter(binmin+1);
3489     x[4] = pente->GetBinCenter(binmin+2);
3490     y[0] = pente->GetBinContent(binmin-2);
3491     y[1] = pente->GetBinContent(binmin-1);
3492     y[2] = pente->GetBinContent(binmin);
3493     y[3] = pente->GetBinContent(binmin+1);
3494     y[4] = pente->GetBinContent(binmin+2);
3495     //Calcul the polynome de Lagrange
3496     c = CalculPolynomeLagrange4(x,y);
3497     //richtung +/-
3498     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
3499        (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
3500     if(((binmin+3) <= (nbins-1)) &&
3501        (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
3502        ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
3503        (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) put = kFALSE;
3504     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
3505        (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) case1 = kTRUE;
3506     if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
3507        (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case4 = kTRUE;
3508   }
3509   //case binmin = nbins-2
3510   //pol3 case 1
3511   if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3512      (case1)){
3513     min = pente->GetBinCenter(binmin-2);
3514     max = pente->GetBinCenter(binmin+1);
3515     x[0] = pente->GetBinCenter(binmin-2);
3516     x[1] = pente->GetBinCenter(binmin-1);
3517     x[2] = pente->GetBinCenter(binmin);
3518     x[3] = pente->GetBinCenter(binmin+1);
3519     y[0] = pente->GetBinContent(binmin-2);
3520     y[1] = pente->GetBinContent(binmin-1);
3521     y[2] = pente->GetBinContent(binmin);
3522     y[3] = pente->GetBinContent(binmin+1);
3523     //Calcul the polynome de Lagrange
3524     c = CalculPolynomeLagrange3(x,y);
3525     //richtung +: nothing
3526     //richtung -
3527     if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case2 = kTRUE;
3528   }
3529   //pol3 case 4
3530   if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3531      (case4)){
3532     min = pente->GetBinCenter(binmin-1);
3533     max = pente->GetBinCenter(binmin+2);
3534     x[0] = pente->GetBinCenter(binmin-1);
3535     x[1] = pente->GetBinCenter(binmin);
3536     x[2] = pente->GetBinCenter(binmin+1);
3537     x[3] = pente->GetBinCenter(binmin+2);
3538     y[0] = pente->GetBinContent(binmin-1);
3539     y[1] = pente->GetBinContent(binmin);
3540     y[2] = pente->GetBinContent(binmin+1);
3541     y[3] = pente->GetBinContent(binmin+2);
3542     //Calcul the polynome de Lagrange
3543     c = CalculPolynomeLagrange3(x,y);
3544     //richtung +
3545     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) case2 = kTRUE;
3546   }
3547   //pol2 case 5
3548   if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
3549     min = pente->GetBinCenter(binmin);
3550     max = pente->GetBinCenter(binmin+2);
3551     x[0] = pente->GetBinCenter(binmin);
3552     x[1] = pente->GetBinCenter(binmin+1);
3553     x[2] = pente->GetBinCenter(binmin+2);
3554     y[0] = pente->GetBinContent(binmin);
3555     y[1] = pente->GetBinContent(binmin+1);
3556     y[2] = pente->GetBinContent(binmin+2);
3557     //Calcul the polynome de Lagrange
3558     c = CalculPolynomeLagrange2(x,y);
3559     //richtung +
3560     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) put = kFALSE;
3561   }
3562   //pol2 case 2
3563   if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3564      (case2)){
3565     min = pente->GetBinCenter(binmin-1);
3566     max = pente->GetBinCenter(binmin+1);
3567     x[0] = pente->GetBinCenter(binmin-1);
3568     x[1] = pente->GetBinCenter(binmin);
3569     x[2] = pente->GetBinCenter(binmin+1);
3570     y[0] = pente->GetBinContent(binmin-1);
3571     y[1] = pente->GetBinContent(binmin);
3572     y[2] = pente->GetBinContent(binmin+1);
3573     //Calcul the polynome de Lagrange
3574     c = CalculPolynomeLagrange2(x,y);
3575     //richtung +: nothing
3576     //richtung -: nothing
3577   }
3578   //case binmin = nbins-1
3579   //pol2 case 0
3580   if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
3581     min = pente->GetBinCenter(binmin-2);
3582     max = pente->GetBinCenter(binmin);
3583     x[0] = pente->GetBinCenter(binmin-2);
3584     x[1] = pente->GetBinCenter(binmin-1);
3585     x[2] = pente->GetBinCenter(binmin);
3586     y[0] = pente->GetBinContent(binmin-2);
3587     y[1] = pente->GetBinContent(binmin-1);
3588     y[2] = pente->GetBinContent(binmin);
3589     //Calcul the polynome de Lagrange
3590     c = CalculPolynomeLagrange2(x,y);
3591     //AliInfo("At the limit for the drift!");
3592     //fluctuation too big!
3593     //richtung +: nothing
3594     //richtung -
3595     if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
3596   }
3597   if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
3598     put = kFALSE;
3599     AliInfo("At the limit for the drift and not usable!");
3600   }
3601
3602   //pass
3603   if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
3604     put = kFALSE;
3605     AliInfo("For the drift...problem!");
3606   }
3607   //pass but should not happen
3608   if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+4, projPH->GetNbinsX()))){
3609     put = kFALSE;
3610     AliInfo("For the drift...problem!");
3611   }
3612
3613   if(put) {
3614     polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
3615     polynome->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3616     //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
3617     Double_t step = (max-min)/1000;
3618     Double_t l = min;
3619     Double_t minvalue = 0.0;
3620     Double_t placeminimum = min;
3621     for(Int_t o = 0; o < 1000; o++){
3622       if(o == 0) minvalue = polynome->Eval(l);
3623       if(minvalue > (polynome->Eval(l))){
3624         minvalue = polynome->Eval(l);
3625         placeminimum = l;
3626       }
3627       l += step;
3628     }
3629     fPhd[2] = placeminimum;
3630   }
3631
3632   Float_t fPhdt0  = 0.0;
3633   Float_t t0Shift = 0.0;
3634   if(fTakeTheMaxPH) {
3635     fPhdt0 = fPhd[1];
3636     t0Shift = fT0Shift1;
3637   }
3638   else {
3639     fPhdt0 = fPhd[0];
3640     t0Shift = fT0Shift0;
3641   }
3642
3643   if ((fPhd[2] > fPhd[0]) &&
3644       (fPhd[2] > fPhd[1]) &&
3645       (fPhd[1] > fPhd[0]) &&
3646       (put)) {
3647     fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
3648     fNumberFitSuccess++;
3649     if (fPhdt0 >= 0.0) {
3650       fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
3651       if (fCurrentCoef2[0] < -1.0) {
3652         fCurrentCoef2[0] = fCurrentCoef2[1];
3653       }
3654     }
3655     else {
3656       fCurrentCoef2[0] = fCurrentCoef2[1];
3657     }
3658   }
3659   else {
3660     fCurrentCoef[0]      = -TMath::Abs(fCurrentCoef[1]);
3661     fCurrentCoef2[0]     = fCurrentCoef2[1];
3662     //printf("Fit failed!\n");
3663   }
3664
3665   if (fDebugLevel == 1) {
3666     TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
3667     cpentei->cd();
3668     projPH->Draw();
3669     line->SetLineColor(2);
3670     line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
3671     line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
3672     line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
3673     AliInfo(Form("fPhd[0] (beginning of the signal): %f"                  ,(Float_t) fPhd[0]));
3674     AliInfo(Form("fPhd[1] (end of the amplification region): %f"          ,(Float_t) fPhd[1]));
3675     AliInfo(Form("fPhd[2] (end of the drift region): %f"                  ,(Float_t) fPhd[2]));
3676     AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
3677     TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
3678     cpentei2->cd();
3679     pentea->Draw();
3680     TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
3681     cpentei3->cd();
3682     pente->Draw();
3683   }
3684   else {
3685     delete pentea;
3686     delete pente;
3687     delete polynome;
3688     delete polynomea;
3689     delete polynomeb;
3690   }
3691
3692   projPH->SetDirectory(0);
3693
3694 }
3695
3696 //_____________________________________________________________________________
3697 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
3698 {
3699   //
3700   // Fit methode for the drift velocity
3701   //
3702
3703   // Constants
3704   const Float_t kDrWidth = AliTRDgeometry::DrThick();
3705
3706   // Some variables
3707   TAxis   *xpph   = projPH->GetXaxis();
3708   Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
3709
3710   TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
3711   fPH->SetParameter(0,0.469);     // Scaling
3712   fPH->SetParameter(1,0.18);      // Start
3713   fPH->SetParameter(2,0.0857325); // AR
3714   fPH->SetParameter(3,1.89);      // DR
3715   fPH->SetParameter(4,0.08);      // QA/QD
3716   fPH->SetParameter(5,0.0);       // Baseline
3717
3718   TLine *line = new TLine();
3719
3720   fCurrentCoef[0]     = 0.0;
3721   fCurrentCoef2[0]    = 0.0;
3722   fCurrentCoefE       = 0.0;
3723   fCurrentCoefE2      = 0.0;
3724
3725   if (idect%fFitPHPeriode == 0) {
3726
3727     AliInfo(Form("The detector %d will be fitted",idect));
3728     fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
3729     fPH->SetParameter(1,fPhd[0] - 0.1);                                                                 // Start
3730     fPH->SetParameter(2,fPhd[1] - fPhd[0]);                                                             // AR
3731     fPH->SetParameter(3,fPhd[2] - fPhd[1]);                                                             // DR
3732     fPH->SetParameter(4,0.225);                                                                         // QA/QD
3733     fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
3734
3735     if (fDebugLevel != 1) {
3736       projPH->Fit(fPH,"0M","",0.0,upedge);
3737     }
3738     else {
3739       TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
3740       cpente->cd();
3741       projPH->Fit(fPH,"M+","",0.0,upedge);
3742       projPH->Draw("E0");
3743       line->SetLineColor(4);
3744       line->DrawLine(fPH->GetParameter(1)
3745                     ,0
3746                     ,fPH->GetParameter(1)
3747                     ,projPH->GetMaximum());
3748       line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
3749                     ,0
3750                     ,fPH->GetParameter(1)+fPH->GetParameter(2)
3751                     ,projPH->GetMaximum());
3752       line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
3753                     ,0
3754                     ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
3755                     ,projPH->GetMaximum());
3756     }
3757
3758     if (fPH->GetParameter(3) != 0) {
3759       fNumberFitSuccess++;
3760       fCurrentCoef[0]    = kDrWidth / (fPH->GetParameter(3));
3761       fCurrentCoefE      = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
3762       fCurrentCoef2[0]   = fPH->GetParameter(1);
3763       fCurrentCoefE2     = fPH->GetParError(1);
3764     }
3765     else {
3766       fCurrentCoef[0]     = -TMath::Abs(fCurrentCoef[1]);
3767       fCurrentCoef2[0]    = fCurrentCoef2[1];
3768     }
3769
3770   }
3771   else {
3772
3773     // Put the default value
3774     fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
3775     fCurrentCoef2[0] = fCurrentCoef2[1];
3776   }
3777
3778   if (fDebugLevel != 1) {
3779     delete fPH;
3780   }
3781
3782 }
3783 //_____________________________________________________________________________
3784 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
3785 {
3786   //
3787   // Fit methode for the sigma of the pad response function
3788   //
3789
3790   TVectorD param(3);
3791
3792   fCurrentCoef[0]  = 0.0;
3793   fCurrentCoefE = 0.0;
3794
3795   Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,&param);
3796
3797   if(ret == -4){
3798     fCurrentCoef[0] = -fCurrentCoef[1];
3799     return kFALSE;
3800   }
3801   else {
3802     fNumberFitSuccess++;
3803     fCurrentCoef[0] = param[2];
3804     fCurrentCoefE   = ret;
3805     return kTRUE;
3806   }
3807 }
3808 //_____________________________________________________________________________
3809 Double_t AliTRDCalibraFit::FitGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, Bool_t kError)
3810 {
3811   //
3812   // Fit methode for the sigma of the pad response function
3813   //
3814
3815   //We should have at least 3 points
3816   if(nBins <=3) return -4.0;
3817
3818   TLinearFitter fitter(3,"pol2");
3819   fitter.StoreData(kFALSE);
3820   fitter.ClearPoints();
3821   TVectorD  par(3);
3822   Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
3823   Float_t entries = 0;
3824   Int_t   nbbinwithentries = 0;
3825   for (Int_t i=0; i<nBins; i++){
3826     entries+=arraye[i];
3827     if(arraye[i] > 15) nbbinwithentries++;
3828     //printf("entries for i %d: %f\n",i,arraye[i]);
3829   }
3830   if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
3831   //printf("entries %f\n",entries);
3832   //printf("nbbinwithentries %d\n",nbbinwithentries);
3833
3834   Int_t npoints=0;
3835   Float_t errorm = 0.0;
3836   Float_t errorn = 0.0;
3837   Float_t error  = 0.0;
3838
3839   //
3840   for (Int_t ibin=0;ibin<nBins; ibin++){
3841       Float_t entriesI = arraye[ibin];
3842       Float_t valueI   = arraym[ibin];
3843       Double_t xcenter = 0.0;
3844       Float_t  val     = 0.0;
3845       if ((entriesI>15) && (valueI>0.0)){
3846         xcenter = xMin+(ibin+0.5)*binWidth;
3847         errorm   = 0.0;
3848         errorn   = 0.0;
3849         error    = 0.0;
3850         if(!kError){
3851           if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
3852           if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
3853           error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
3854         }
3855         else{
3856           if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
3857           if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
3858           error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
3859         }
3860         if(error == 0.0) continue;
3861         val      = TMath::Log(Float_t(valueI));
3862         fitter.AddPoint(&xcenter,val,error);
3863         npoints++;
3864       }
3865
3866       if(fDebugLevel > 1){
3867
3868       if ( !fDebugStreamer ) {
3869         //debug stream
3870         TDirectory *backup = gDirectory;
3871         fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
3872         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3873       }
3874
3875       Int_t    detector     = fCountDet;
3876       Int_t    plane        = GetPlane(fCountDet);
3877       Int_t    group        = ibin;
3878
3879       (* fDebugStreamer) << "FitGausMIFill"<<
3880         "detector="<<detector<<
3881         "plane="<<plane<<
3882         "nbins="<<nBins<<
3883         "group="<<group<<
3884         "entriesI="<<entriesI<<
3885         "valueI="<<valueI<<
3886         "val="<<val<<
3887         "xcenter="<<xcenter<<
3888         "errorm="<<errorm<<
3889         "errorn="<<errorn<<
3890         "error="<<error<<
3891         "kError="<<kError<<
3892         "\n";
3893     }
3894
3895   }
3896
3897   if(npoints <=3) return -4.0;
3898
3899   Double_t chi2 = 0;
3900   if (npoints>3){
3901     fitter.Eval();
3902     fitter.GetParameters(par);
3903     chi2 = fitter.GetChisquare()/Float_t(npoints);
3904
3905
3906     if (!param)  param  = new TVectorD(3);
3907     if(par[2] == 0.0) return -4.0;
3908     Double_t  x      = TMath::Sqrt(TMath::Abs(-2*par[2]));
3909     Double_t deltax = (fitter.GetParError(2))/x;
3910     Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
3911     chi2 = errorparam2;
3912
3913     (*param)[1] = par[1]/(-2.*par[2]);
3914     (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
3915     Double_t lnparam0 = par[0]+ par[1]* (*param)[1] +  par[2]*(*param)[1]*(*param)[1];
3916     if ( lnparam0>307 ) return -4;
3917     (*param)[0] = TMath::Exp(lnparam0);
3918
3919     if(fDebugLevel > 1){
3920
3921       if ( !fDebugStreamer ) {
3922         //debug stream
3923         TDirectory *backup = gDirectory;
3924         fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
3925         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3926       }
3927
3928       Int_t    detector     = fCountDet;
3929       Int_t    plane        = GetPlane(fCountDet);
3930
3931
3932       (* fDebugStreamer) << "FitGausMIFit"<<
3933         "detector="<<detector<<
3934         "plane="<<plane<<
3935         "nbins="<<nBins<<
3936         "errorsigma="<<chi2<<
3937         "mean="<<(*param)[1]<<
3938         "sigma="<<(*param)[2]<<
3939         "constant="<<(*param)[0]<<
3940         "\n";
3941     }
3942   }
3943
3944   if((chi2/(*param)[2]) > 0.1){
3945     if(kError){
3946       chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
3947     }
3948     else return -4.0;
3949   }
3950
3951   if(fDebugLevel == 1){
3952     TString name("PRF");
3953     name += (Int_t)xMin;
3954     name += (Int_t)xMax;
3955     TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
3956     c1->cd();
3957     name += "histo";
3958     TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
3959     for(Int_t k = 0; k < nBins; k++){
3960       histo->SetBinContent(k+1,arraym[k]);
3961       histo->SetBinError(k+1,arrayme[k]);
3962     }
3963     histo->Draw();
3964     name += "functionf";
3965     TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
3966     f1->SetParameter(0, (*param)[0]);
3967     f1->SetParameter(1, (*param)[1]);
3968     f1->SetParameter(2, (*param)[2]);
3969     f1->Draw("same");
3970   }
3971
3972
3973   return chi2;
3974
3975 }
3976 //_____________________________________________________________________________
3977 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
3978 {
3979   //
3980   // Fit methode for the sigma of the pad response function
3981   //
3982
3983   fCurrentCoef[0]  = 0.0;
3984   fCurrentCoefE = 0.0;
3985
3986   if (fDebugLevel != 1) {
3987     projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
3988   }
3989   else {
3990     TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
3991     cfit->cd();
3992     projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
3993     projPRF->Draw();
3994   }
3995   fCurrentCoef[0]  = projPRF->GetFunction("gaus")->GetParameter(2);
3996   fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
3997   if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
3998   else {
3999     fNumberFitSuccess++;
4000   }
4001 }
4002 //_____________________________________________________________________________
4003 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
4004 {
4005   //
4006   // Fit methode for the sigma of the pad response function
4007   //
4008   fCurrentCoef[0]   = 0.0;
4009   fCurrentCoefE  = 0.0;
4010   if (fDebugLevel == 1) {
4011     TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
4012     cfit->cd();
4013     projPRF->Draw();
4014   }
4015   fCurrentCoef[0] = projPRF->GetRMS();
4016   if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
4017   else {
4018     fNumberFitSuccess++;
4019   }
4020 }
4021 //_____________________________________________________________________________
4022 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
4023 {
4024   //
4025   // Fit methode for the sigma of the pad response function with 2*nbg tan bins
4026   //
4027
4028   TLinearFitter linearfitter = TLinearFitter(3,"pol2");
4029
4030
4031   Int_t   nbins    = (Int_t)(nybins/(2*nbg));
4032   Float_t lowedge  = -3.0*nbg;
4033   Float_t upedge   = lowedge + 3.0;
4034   Int_t   offset   = 0;
4035   Int_t   npoints  = 0;
4036   Double_t xvalues = -0.2*nbg+0.1;
4037   Double_t y       = 0.0;
4038   Int_t   total    = 2*nbg;
4039
4040
4041   for(Int_t k = 0; k < total; k++){
4042     if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
4043       npoints++;
4044       y = fCurrentCoef[0]*fCurrentCoef[0];
4045       linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
4046     }
4047
4048     if(fDebugLevel > 1){
4049
4050       if ( !fDebugStreamer ) {
4051         //debug stream
4052         TDirectory *backup = gDirectory;
4053         fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
4054         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
4055       }
4056
4057       Int_t    detector     = fCountDet;
4058       Int_t    plane        = GetPlane(fCountDet);
4059       Int_t    nbtotal      = total;
4060       Int_t    group        = k;
4061       Float_t  low          = lowedge;
4062       Float_t  up           = upedge;
4063       Float_t  tnp          = xvalues;
4064       Float_t  wid          = fCurrentCoef[0];
4065       Float_t  widfE        = fCurrentCoefE;
4066
4067       (* fDebugStreamer) << "FitTnpRange0"<<
4068         "detector="<<detector<<
4069         "plane="<<plane<<
4070         "nbtotal="<<nbtotal<<
4071         "group="<<group<<
4072         "low="<<low<<
4073         "up="<<up<<
4074         "offset="<<offset<<
4075         "tnp="<<tnp<<
4076         "wid="<<wid<<
4077         "widfE="<<widfE<<
4078         "\n";
4079     }
4080
4081     offset  += nbins;
4082     lowedge += 3.0;
4083     upedge  += 3.0;
4084     xvalues += 0.2;
4085
4086   }
4087
4088   fCurrentCoefE = 0.0;
4089   fCurrentCoef[0] = 0.0;
4090
4091   //printf("npoints\n",npoints);
4092
4093   if(npoints < 3){
4094     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4095   }
4096   else{
4097
4098     TVectorD pars0;
4099     linearfitter.Eval();
4100     linearfitter.GetParameters(pars0);
4101     Double_t pointError0  =  TMath::Sqrt(linearfitter.GetChisquare()/npoints);
4102     Double_t errorsx0     =  linearfitter.GetParError(2)*pointError0;
4103     Double_t min0         = 0.0;
4104     Double_t ermin0       = 0.0;
4105     //Double_t prfe0      = 0.0;
4106     Double_t prf0         = 0.0;
4107     if((pars0[2] > 0.0) && (pars0[1] != 0.0)) {
4108       min0 = -pars0[1]/(2*pars0[2]);
4109       ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
4110       prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
4111       if(prf0 > 0.0) {
4112         /*
4113           prfe0 = linearfitter->GetParError(0)*pointError0
4114           +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
4115           +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
4116           prfe0 = prfe0/(2*TMath::Sqrt(prf0));
4117           fCurrentCoefE   = (Float_t) prfe0;
4118         */
4119         fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
4120       }
4121       else{
4122         fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4123       }
4124     }
4125     else {
4126       fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4127     }
4128
4129     if(fDebugLevel > 1){
4130
4131       if ( !fDebugStreamer ) {
4132         //debug stream
4133         TDirectory *backup = gDirectory;
4134         fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
4135         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
4136       }
4137
4138       Int_t    detector     = fCountDet;
4139       Int_t    plane        = GetPlane(fCountDet);
4140       Int_t    nbtotal      = total;
4141       Double_t colsize[6]   = {0.635,0.665,0.695,0.725,0.755,0.785};
4142       Double_t sigmax       = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[plane];
4143
4144       (* fDebugStreamer) << "FitTnpRange1"<<
4145         "detector="<<detector<<
4146         "plane="<<plane<<
4147         "nbtotal="<<nbtotal<<
4148         "par0="<<pars0[0]<<
4149         "par1="<<pars0[1]<<
4150         "par2="<<pars0[2]<<
4151         "npoints="<<npoints<<
4152         "sigmax="<<sigmax<<
4153         "tan="<<min0<<
4154         "sigmaprf="<<fCurrentCoef[0]<<
4155         "sigprf="<<fCurrentCoef[1]<<
4156         "\n";
4157     }
4158
4159   }
4160
4161 }
4162 //_____________________________________________________________________________
4163 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
4164 {
4165   //
4166   // Only mean methode for the gain factor
4167   //
4168
4169   fCurrentCoef[0] = mean;
4170   fCurrentCoefE   = 0.0;
4171   if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
4172   if (fDebugLevel == 1) {
4173     TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
4174     cpmean->cd();
4175     projch->Draw();
4176   }
4177   CalculChargeCoefMean(kTRUE);
4178   fNumberFitSuccess++;
4179 }
4180 //_____________________________________________________________________________
4181 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
4182 {
4183   //
4184   // mean w methode for the gain factor
4185   //
4186
4187   //Number of bins
4188   Int_t nybins = projch->GetNbinsX();
4189
4190   //The weight function
4191   Double_t a = 0.00228515;
4192   Double_t b = -0.00231487;
4193   Double_t c = 0.00044298;
4194   Double_t d = -0.00379239;
4195   Double_t e = 0.00338349;
4196
4197 //   0 |0.00228515
4198 //    1 |-0.00231487
4199 //    2 |0.00044298
4200 //    3 |-0.00379239
4201 //    4 |0.00338349
4202
4203
4204
4205   //A arbitrary error for the moment
4206   fCurrentCoefE = 0.0;
4207   fCurrentCoef[0] = 0.0;
4208
4209   //Calcul
4210   Double_t sumw = 0.0;
4211   Double_t sum = 0.0;
4212   Float_t sumAll   = (Float_t) nentries;
4213   Int_t sumCurrent = 0;
4214   for(Int_t k = 0; k <nybins; k++){
4215     Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
4216     if (fraction>0.95) break;
4217     Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
4218       e*fraction*fraction*fraction*fraction;
4219     sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
4220     sum  += weight*projch->GetBinContent(k+1);
4221     sumCurrent += (Int_t) projch->GetBinContent(k+1);
4222     //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
4223   }
4224   if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
4225
4226   if (fDebugLevel == 1) {
4227     TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
4228     cpmeanw->cd();
4229     projch->Draw();
4230   }
4231   fNumberFitSuccess++;
4232   CalculChargeCoefMean(kTRUE);
4233 }
4234 //_____________________________________________________________________________
4235 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
4236 {
4237   //
4238   // mean w methode for the gain factor
4239   //
4240
4241   //Number of bins
4242   Int_t nybins = projch->GetNbinsX();
4243
4244   //The weight function
4245   Double_t a = 0.00228515;
4246   Double_t b = -0.00231487;
4247   Double_t c = 0.00044298;
4248   Double_t d = -0.00379239;
4249   Double_t e = 0.00338349;
4250
4251 //   0 |0.00228515
4252 //    1 |-0.00231487
4253 //    2 |0.00044298
4254 //    3 |-0.00379239
4255 //    4 |0.00338349
4256
4257
4258
4259   //A arbitrary error for the moment
4260   fCurrentCoefE = 0.0;
4261   fCurrentCoef[0] = 0.0;
4262
4263   //Calcul
4264   Double_t sumw = 0.0;
4265   Double_t sum = 0.0;
4266   Int_t sumCurrent = 0;
4267   for(Int_t k = 0; k <nybins; k++){
4268     Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
4269     if (fraction>0.95) break;
4270     Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
4271       e*fraction*fraction*fraction*fraction;
4272     sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
4273     sum  += weight*projch->GetBinContent(k+1);
4274     sumCurrent += (Int_t) projch->GetBinContent(k+1);
4275     //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
4276   }
4277   if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
4278
4279   if (fDebugLevel == 1) {
4280     TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
4281     cpmeanw->cd();
4282     projch->Draw();
4283   }
4284   fNumberFitSuccess++;
4285 }
4286 //_____________________________________________________________________________
4287 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
4288 {
4289   //
4290   // Fit methode for the gain factor
4291   //
4292
4293   fCurrentCoef[0]  = 0.0;
4294   fCurrentCoefE    = 0.0;
4295   Double_t chisqrl = 0.0;
4296   Double_t chisqrg = 0.0;
4297   Double_t chisqr  = 0.0;
4298   TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
4299
4300   projch->Fit("landau","0",""
4301              ,(Double_t) mean/fBeginFitCharge
4302              ,projch->GetBinCenter(projch->GetNbinsX()));
4303   Double_t l3P0         = projch->GetFunction("landau")->GetParameter(0);
4304   Double_t l3P1         = projch->GetFunction("landau")->GetParameter(1);
4305   Double_t l3P2         = projch->GetFunction("landau")->GetParameter(2);
4306   chisqrl = projch->GetFunction("landau")->GetChisquare();
4307
4308   projch->Fit("gaus","0",""
4309               ,(Double_t) mean/fBeginFitCharge
4310               ,projch->GetBinCenter(projch->GetNbinsX()));
4311   Double_t g3P0         = projch->GetFunction("gaus")->GetParameter(0);
4312   Double_t g3P2         = projch->GetFunction("gaus")->GetParameter(2);
4313   chisqrg = projch->GetFunction("gaus")->GetChisquare();
4314
4315   fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
4316   if (fDebugLevel != 1) {
4317     projch->Fit("fLandauGaus","0",""
4318                 ,(Double_t) mean/fBeginFitCharge
4319                 ,projch->GetBinCenter(projch->GetNbinsX()));
4320     chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
4321   }
4322   else  {
4323     TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
4324     cp->cd();
4325     projch->Fit("fLandauGaus","+",""
4326                 ,(Double_t) mean/fBeginFitCharge
4327                 ,projch->GetBinCenter(projch->GetNbinsX()));
4328     chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
4329     projch->Draw();
4330     fLandauGaus->Draw("same");
4331   }
4332
4333   if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
4334     //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
4335     fNumberFitSuccess++;
4336     CalculChargeCoefMean(kTRUE);
4337     fCurrentCoef[0]  = projch->GetFunction("fLandauGaus")->GetParameter(1);
4338     fCurrentCoefE    = projch->GetFunction("fLandauGaus")->GetParError(1);
4339   }
4340   else {
4341     CalculChargeCoefMean(kFALSE);
4342     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4343   }
4344
4345   if (fDebugLevel != 1) {
4346     delete fLandauGaus;
4347   }
4348
4349 }
4350 //_____________________________________________________________________________
4351 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
4352 {
4353   //
4354   // Fit methode for the gain factor more time consuming
4355   //
4356
4357
4358   //Some parameters to initialise
4359   Double_t widthLandau, widthGaus, mPV, integral;
4360   Double_t chisquarel = 0.0;
4361   Double_t chisquareg = 0.0;
4362   projch->Fit("landau","0M+",""
4363               ,(Double_t) mean/6
4364               ,projch->GetBinCenter(projch->GetNbinsX()));
4365   widthLandau  = projch->GetFunction("landau")->GetParameter(2);
4366   chisquarel = projch->GetFunction("landau")->GetChisquare();
4367   projch->Fit("gaus","0M+",""
4368               ,(Double_t) mean/6
4369               ,projch->GetBinCenter(projch->GetNbinsX()));
4370   widthGaus    = projch->GetFunction("gaus")->GetParameter(2);
4371   chisquareg = projch->GetFunction("gaus")->GetChisquare();
4372
4373   mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
4374   integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
4375
4376   // Setting fit range and start values
4377   Double_t fr[2];
4378   //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
4379   //Double_t sv[4]   = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
4380   Double_t sv[4]   = { widthLandau, mPV, integral, widthGaus};
4381   Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
4382   Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
4383   Double_t fp[4]   = { 1.0, 1.0, 1.0, 1.0 };
4384   Double_t fpe[4]  = { 1.0, 1.0, 1.0, 1.0 };
4385   fr[0]            = 0.3 * mean;
4386   fr[1]            = 3.0 * mean;
4387   fCurrentCoef[0]  = 0.0;
4388   fCurrentCoefE    = 0.0;
4389
4390   Double_t chisqr;
4391   Int_t    ndf;
4392   TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
4393                                 ,&pllo[0],&plhi[0]
4394                                 ,&fp[0],&fpe[0]
4395                                 ,&chisqr,&ndf);
4396
4397   Double_t projchPeak;
4398   Double_t projchFWHM;
4399   LanGauPro(fp,projchPeak,projchFWHM);
4400
4401   if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
4402     //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
4403     fNumberFitSuccess++;
4404     CalculChargeCoefMean(kTRUE);
4405     fCurrentCoef[0]  = fp[1];
4406     fCurrentCoefE = fpe[1];
4407     //chargeCoefE2 = chisqr;
4408   }
4409   else {
4410     CalculChargeCoefMean(kFALSE);
4411     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4412   }
4413   if (fDebugLevel == 1) {
4414     AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
4415     TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
4416     cpy->cd();
4417     projch->Draw();
4418     fitsnr->Draw("same");
4419   }
4420   else {
4421     delete fitsnr;
4422   }
4423 }
4424 //_____________________________________________________________________________
4425 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(Double_t *x, Double_t *y) const
4426 {
4427   //
4428   // Calcul the coefficients of the polynome passant par ces trois points de degre 2
4429   //
4430   Double_t *c = new Double_t[5];
4431   Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
4432   Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
4433   Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
4434
4435   c[4] = 0.0;
4436   c[3] = 0.0;
4437   c[2] = x0+x1+x2;
4438   c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
4439   c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
4440
4441   return c;
4442
4443
4444 }
4445
4446 //_____________________________________________________________________________
4447 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(Double_t *x, Double_t *y) const
4448 {
4449   //
4450   // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
4451   //
4452   Double_t *c = new Double_t[5];
4453   Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
4454   Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
4455   Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
4456   Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
4457
4458   c[4] = 0.0;
4459   c[3] = x0+x1+x2+x3;
4460   c[2] = -(x0*(x[1]+x[2]+x[3])
4461            +x1*(x[0]+x[2]+x[3])
4462            +x2*(x[0]+x[1]+x[3])
4463            +x3*(x[0]+x[1]+x[2]));
4464   c[1] = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
4465           +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
4466           +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
4467           +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
4468
4469   c[0] = -(x0*x[1]*x[2]*x[3]
4470           +x1*x[0]*x[2]*x[3]
4471           +x2*x[0]*x[1]*x[3]
4472           +x3*x[0]*x[1]*x[2]);
4473
4474
4475   return c;
4476
4477
4478 }
4479
4480 //_____________________________________________________________________________
4481 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(Double_t *x, Double_t *y) const
4482 {
4483   //
4484   // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
4485   //
4486   Double_t *c = new Double_t[5];
4487   Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
4488   Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
4489   Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
4490   Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
4491   Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
4492
4493
4494   c[4] = x0+x1+x2+x3+x4;
4495   c[3] = -(x0*(x[1]+x[2]+x[3]+x[4])
4496            +x1*(x[0]+x[2]+x[3]+x[4])
4497            +x2*(x[0]+x[1]+x[3]+x[4])
4498            +x3*(x[0]+x[1]+x[2]+x[4])
4499            +x4*(x[0]+x[1]+x[2]+x[3]));
4500   c[2] = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
4501           +x1*(x[0]*x[2]+x[0]*x[3]+x[0]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
4502           +x2*(x[0]*x[1]+x[0]*x[3]+x[0]*x[4]+x[1]*x[3]+x[1]*x[4]+x[3]*x[4])
4503           +x3*(x[0]*x[1]+x[0]*x[2]+x[0]*x[4]+x[1]*x[2]+x[1]*x[4]+x[2]*x[4])
4504           +x4*(x[0]*x[1]+x[0]*x[2]+x[0]*x[3]+x[1]*x[2]+x[1]*x[3]+x[2]*x[3]));
4505
4506   c[1] = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
4507           +x1*(x[0]*x[2]*x[3]+x[0]*x[2]*x[4]+x[0]*x[3]*x[4]+x[2]*x[3]*x[4])
4508           +x2*(x[0]*x[1]*x[3]+x[0]*x[1]*x[4]+x[0]*x[3]*x[4]+x[1]*x[3]*x[4])
4509           +x3*(x[0]*x[1]*x[2]+x[0]*x[1]*x[4]+x[0]*x[2]*x[4]+x[1]*x[2]*x[4])
4510           +x4*(x[0]*x[1]*x[2]+x[0]*x[1]*x[3]+x[0]*x[2]*x[3]+x[1]*x[2]*x[3]));
4511
4512   c[0] = (x0*x[1]*x[2]*x[3]*x[4]
4513           +x1*x[0]*x[2]*x[3]*x[4]
4514           +x2*x[0]*x[1]*x[3]*x[4]
4515           +x3*x[0]*x[1]*x[2]*x[4]
4516           +x4*x[0]*x[1]*x[2]*x[3]);
4517
4518   return c;
4519
4520
4521 }
4522 //_____________________________________________________________________________
4523 void AliTRDCalibraFit::NormierungCharge()
4524 {
4525   //
4526   // Normalisation of the gain factor resulting for the fits
4527   //
4528
4529   // Calcul of the mean of choosen method by fFitChargeNDB
4530   Double_t sum         = 0.0;
4531   //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
4532   for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
4533     Int_t    total    = 0;
4534     Int_t    detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
4535     Float_t *coef     = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
4536     //printf("detector %d coef[0] %f\n",detector,coef[0]);
4537     if (GetChamber(detector) == 2) {
4538       total = 1728;
4539     }
4540     if (GetChamber(detector) != 2) {
4541       total = 2304;
4542     }
4543     for (Int_t j = 0; j < total; j++) {
4544       if (coef[j] >= 0) {
4545         sum += coef[j];
4546       }
4547     }
4548   }
4549
4550   if (sum > 0) {
4551     fScaleFitFactor = fScaleFitFactor / sum;
4552   }
4553   else {
4554     fScaleFitFactor = 1.0;
4555   }
4556
4557   //methode de boeuf mais bon...
4558   Double_t scalefactor = fScaleFitFactor;
4559
4560   if(fDebugLevel > 1){
4561
4562     if ( !fDebugStreamer ) {
4563       //debug stream
4564       TDirectory *backup = gDirectory;
4565       fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
4566       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
4567     }
4568     (* fDebugStreamer) << "NormierungCharge"<<
4569       "scalefactor="<<scalefactor<<
4570       "\n";
4571     }
4572 }
4573 //_____________________________________________________________________________
4574 TH1I *AliTRDCalibraFit::ReBin(TH1I *hist) const
4575 {
4576   //
4577   // Rebin of the 1D histo for the gain calibration if needed.
4578   // you have to choose fRebin, divider of fNumberBinCharge
4579   //
4580
4581  TAxis *xhist  = hist->GetXaxis();
4582  TH1I  *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
4583                                         ,xhist->GetBinLowEdge(1)
4584                                         ,xhist->GetBinUpEdge(xhist->GetNbins()));
4585
4586  AliInfo(Form("fRebin: %d",fRebin));
4587  Int_t i = 1;
4588  for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
4589    Double_t sum = 0.0;
4590    for (Int_t ji = i; ji < i+fRebin; ji++) {
4591      sum += hist->GetBinContent(ji);
4592    }
4593    sum = sum / fRebin;
4594    rehist->SetBinContent(k,sum);
4595    i += fRebin;
4596  }
4597
4598  return rehist;
4599
4600 }
4601
4602 //_____________________________________________________________________________
4603 TH1F *AliTRDCalibraFit::ReBin(TH1F *hist) const
4604 {
4605   //
4606   // Rebin of the 1D histo for the gain calibration if needed
4607   // you have to choose fRebin divider of fNumberBinCharge
4608   //
4609
4610   TAxis *xhist  = hist->GetXaxis();
4611   TH1F  *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
4612                                          ,xhist->GetBinLowEdge(1)
4613                                          ,xhist->GetBinUpEdge(xhist->GetNbins()));
4614
4615   AliInfo(Form("fRebin: %d",fRebin));
4616   Int_t i = 1;
4617   for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
4618     Double_t sum = 0.0;
4619     for (Int_t ji = i; ji < i+fRebin; ji++) {
4620       sum += hist->GetBinContent(ji);
4621     }
4622     sum = sum/fRebin;
4623     rehist->SetBinContent(k,sum);
4624     i += fRebin;
4625   }
4626
4627   return rehist;
4628
4629 }
4630
4631 //_____________________________________________________________________________
4632 TH1F *AliTRDCalibraFit::CorrectTheError(TGraphErrors *hist)
4633 {
4634   //
4635   // In the case of the vectors method the trees contains TGraphErrors for PH and PRF
4636   // to be able to add them after
4637   // We convert it to a TH1F to be able to applied the same fit function method
4638   // After having called this function you can not add the statistics anymore
4639   //
4640
4641   TH1F *rehist = 0x0;
4642
4643   Int_t nbins       = hist->GetN();
4644   Double_t *x       = hist->GetX();
4645   Double_t *entries = hist->GetEX();
4646   Double_t *mean    = hist->GetY();
4647   Double_t *square  = hist->GetEY();
4648   fEntriesCurrent   = 0;
4649
4650   if (nbins < 2) {
4651     return rehist;
4652   }
4653
4654   Double_t step     = x[1] - x[0];
4655   Double_t minvalue = x[0] - step/2;
4656   Double_t maxvalue = x[(nbins-1)] + step/2;
4657
4658   rehist = new TH1F("projcorrecterror","",nbins,minvalue,maxvalue);
4659
4660   for (Int_t k = 0; k < nbins; k++) {
4661     rehist->SetBinContent(k+1,mean[k]);
4662     if (entries[k] > 0.0) {
4663       fEntriesCurrent += (Int_t) entries[k];
4664       Double_t d = TMath::Abs(square[k] - (mean[k]*mean[k]));
4665       rehist->SetBinError(k+1,TMath::Sqrt(d/entries[k]));
4666     }
4667     else {
4668       rehist->SetBinError(k+1,0.0);
4669     }
4670   }
4671
4672   if(fEntriesCurrent > 0) fNumberEnt++;
4673
4674   return rehist;
4675
4676 }
4677 //
4678 //____________Some basic geometry function_____________________________________
4679 //
4680
4681 //_____________________________________________________________________________
4682 Int_t AliTRDCalibraFit::GetPlane(Int_t d) const
4683 {
4684   //
4685   // Reconstruct the plane number from the detector number
4686   //
4687
4688   return ((Int_t) (d % 6));
4689
4690 }
4691
4692 //_____________________________________________________________________________
4693 Int_t AliTRDCalibraFit::GetChamber(Int_t d) const
4694 {
4695   //
4696   // Reconstruct the chamber number from the detector number
4697   //
4698   Int_t fgkNplan = 6;
4699
4700   return ((Int_t) (d % 30) / fgkNplan);
4701
4702 }
4703
4704 //_____________________________________________________________________________
4705 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
4706 {
4707   //
4708   // Reconstruct the sector number from the detector number
4709   //
4710   Int_t fg = 30;
4711
4712   return ((Int_t) (d / fg));
4713
4714 }
4715
4716 //
4717 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
4718 //
4719 //_______________________________________________________________________________
4720 void AliTRDCalibraFit::ResetVectorFit()
4721 {
4722   //
4723   // Reset the VectorFits
4724   //
4725
4726   fVectorFit.SetOwner();
4727   fVectorFit.Clear();
4728   fVectorFit2.SetOwner();
4729   fVectorFit2.Clear();
4730
4731 }
4732 //
4733 //____________Private Functions________________________________________________
4734 //
4735
4736 //_____________________________________________________________________________
4737 Double_t AliTRDCalibraFit::PH(Double_t *x, Double_t *par)
4738 {
4739   //
4740   // Function for the fit
4741   //
4742
4743   //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
4744
4745   //PARAMETERS FOR FIT PH
4746   // PASAv.4
4747   //fAsymmGauss->SetParameter(0,0.113755);
4748   //fAsymmGauss->SetParameter(1,0.350706);
4749   //fAsymmGauss->SetParameter(2,0.0604244);
4750   //fAsymmGauss->SetParameter(3,7.65596);
4751   //fAsymmGauss->SetParameter(4,1.00124);
4752   //fAsymmGauss->SetParameter(5,0.870597);  // No tail cancelation
4753
4754   Double_t xx = x[0];
4755
4756   if (xx < par[1]) {
4757     return par[5];
4758   }
4759
4760   Double_t dx       = 0.005;
4761   Double_t xs       = par[1];
4762   Double_t ss       = 0.0;
4763   Double_t paras[2] = { 0.0, 0.0 };
4764
4765   while (xs < xx) {
4766     if ((xs >= par[1]) &&
4767         (xs < (par[1]+par[2]))) {
4768       //fAsymmGauss->SetParameter(0,par[0]);
4769       //fAsymmGauss->SetParameter(1,xs);
4770       //ss += fAsymmGauss->Eval(xx);
4771       paras[0] = par[0];
4772       paras[1] = xs;
4773       ss += AsymmGauss(&xx,paras);
4774     }
4775     if ((xs >= (par[1]+par[2])) &&
4776         (xs <  (par[1]+par[2]+par[3]))) {
4777       //fAsymmGauss->SetParameter(0,par[0]*par[4]);
4778       //fAsymmGauss->SetParameter(1,xs);
4779       //ss += fAsymmGauss->Eval(xx);
4780       paras[0] = par[0]*par[4];
4781       paras[1] = xs;
4782       ss += AsymmGauss(&xx,paras);
4783     }
4784     xs += dx;
4785   }
4786
4787   return ss + par[5];
4788
4789 }
4790
4791 //_____________________________________________________________________________
4792 Double_t AliTRDCalibraFit::AsymmGauss(Double_t *x, Double_t *par)
4793 {
4794   //
4795   // Function for the fit
4796   //
4797
4798   //par[0] = normalization
4799   //par[1] = mean
4800   //par[2] = sigma
4801   //norm0  = 1
4802   //par[3] = lambda0
4803   //par[4] = norm1
4804   //par[5] = lambda1
4805
4806   Double_t par1save = par[1];
4807   //Double_t par2save = par[2];
4808   Double_t par2save = 0.0604244;
4809   //Double_t par3save = par[3];
4810   Double_t par3save = 7.65596;
4811   //Double_t par5save = par[5];
4812   Double_t par5save = 0.870597;
4813   Double_t dx       = x[0] - par1save;
4814
4815   Double_t  sigma2  = par2save*par2save;
4816   Double_t  sqrt2   = TMath::Sqrt(2.0);
4817   Double_t  exp1    = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
4818                                * (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
4819   Double_t  exp2    = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
4820                                * (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
4821
4822   //return par[0]*(exp1+par[4]*exp2);
4823   return par[0] * (exp1 + 1.00124 * exp2);
4824
4825 }
4826
4827 //_____________________________________________________________________________
4828 Double_t AliTRDCalibraFit::FuncLandauGaus(Double_t *x, Double_t *par)
4829 {
4830   //
4831   // Sum Landau + Gaus with identical mean
4832   //
4833
4834   Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
4835   //Double_t valGaus   = par[3] * TMath::Gaus(x[0],par[4],par[5]);
4836   Double_t valGaus   = par[3] * TMath::Gaus(x[0],par[1],par[4]);
4837   Double_t val       = valLandau + valGaus;
4838
4839   return val;
4840
4841 }
4842
4843 //_____________________________________________________________________________
4844 Double_t AliTRDCalibraFit::LanGauFun(Double_t *x, Double_t *par)
4845 {
4846   //
4847   // Function for the fit
4848   //
4849   // Fit parameters:
4850   // par[0]=Width (scale) parameter of Landau density
4851   // par[1]=Most Probable (MP, location) parameter of Landau density
4852   // par[2]=Total area (integral -inf to inf, normalization constant)
4853   // par[3]=Width (sigma) of convoluted Gaussian function
4854   //
4855   // In the Landau distribution (represented by the CERNLIB approximation),
4856   // the maximum is located at x=-0.22278298 with the location parameter=0.
4857   // This shift is corrected within this function, so that the actual
4858   // maximum is identical to the MP parameter.
4859   //
4860
4861   // Numeric constants
4862   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
4863   Double_t mpshift  = -0.22278298;       // Landau maximum location
4864
4865   // Control constants
4866   Double_t np       = 100.0;             // Number of convolution steps
4867   Double_t sc       =   5.0;             // Convolution extends to +-sc Gaussian sigmas
4868
4869   // Variables
4870   Double_t xx;
4871   Double_t mpc;
4872   Double_t fland;
4873   Double_t sum = 0.0;
4874   Double_t xlow;
4875   Double_t xupp;
4876   Double_t step;
4877   Double_t i;
4878
4879   // MP shift correction
4880   mpc = par[1] - mpshift * par[0];
4881
4882   // Range of convolution integral
4883   xlow = x[0] - sc * par[3];
4884   xupp = x[0] + sc * par[3];
4885
4886   step = (xupp - xlow) / np;
4887
4888   // Convolution integral of Landau and Gaussian by sum
4889   for (i = 1.0; i <= np/2; i++) {
4890
4891     xx    = xlow + (i-.5) * step;
4892     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
4893     sum  += fland * TMath::Gaus(x[0],xx,par[3]);
4894
4895     xx    = xupp - (i-.5) * step;
4896     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
4897     sum  += fland * TMath::Gaus(x[0],xx,par[3]);
4898
4899   }
4900
4901   return (par[2] * step * sum * invsq2pi / par[3]);
4902
4903 }
4904 //_____________________________________________________________________________
4905 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, Double_t *fitrange, Double_t *startvalues
4906                                       , Double_t *parlimitslo, Double_t *parlimitshi
4907                                       , Double_t *fitparams, Double_t *fiterrors
4908                                       , Double_t *chiSqr, Int_t *ndf) const
4909 {
4910   //
4911   // Function for the fit
4912   //
4913
4914   Int_t i;
4915   Char_t funname[100];
4916
4917   TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
4918   if (ffitold) {
4919     delete ffitold;
4920   }
4921
4922   TF1 *ffit    = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
4923   ffit->SetParameters(startvalues);
4924   ffit->SetParNames("Width","MP","Area","GSigma");
4925
4926   for (i = 0; i < 4; i++) {
4927     ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
4928   }
4929
4930   his->Fit(funname,"RB0");                   // Fit within specified range, use ParLimits, do not plot
4931
4932   ffit->GetParameters(fitparams);            // Obtain fit parameters
4933   for (i = 0; i < 4; i++) {
4934     fiterrors[i] = ffit->GetParError(i);     // Obtain fit parameter errors
4935   }
4936   chiSqr[0] = ffit->GetChisquare();          // Obtain chi^2
4937   ndf[0]    = ffit->GetNDF();                // Obtain ndf
4938
4939   return (ffit);                             // Return fit function
4940
4941 }
4942
4943 //_____________________________________________________________________________
4944 Int_t AliTRDCalibraFit::LanGauPro(Double_t *params, Double_t &maxx, Double_t &fwhm)
4945 {
4946   //
4947   // Function for the fit
4948   //
4949
4950   Double_t p;
4951   Double_t x;
4952   Double_t fy;
4953   Double_t fxr;
4954   Double_t fxl;
4955   Double_t step;
4956   Double_t l;
4957   Double_t lold;
4958
4959   Int_t    i        = 0;
4960   Int_t    maxcalls = 10000;
4961
4962   // Search for maximum
4963   p    = params[1] - 0.1 * params[0];
4964   step = 0.05 * params[0];
4965   lold = -2.0;
4966   l    = -1.0;
4967
4968   while ((l != lold) && (i < maxcalls)) {
4969     i++;
4970     lold = l;
4971     x    = p + step;
4972     l    = LanGauFun(&x,params);
4973     if (l < lold) {
4974       step = -step / 10.0;
4975     }
4976     p += step;
4977   }
4978
4979   if (i == maxcalls) {
4980     return (-1);
4981   }
4982   maxx = x;
4983   fy = l / 2.0;
4984
4985   // Search for right x location of fy
4986   p    = maxx + params[0];
4987   step = params[0];
4988   lold = -2.0;
4989   l    = -1e300;
4990   i    = 0;
4991
4992   while ( (l != lold) && (i < maxcalls) ) {
4993     i++;
4994
4995     lold = l;
4996     x = p + step;
4997     l = TMath::Abs(LanGauFun(&x,params) - fy);
4998
4999     if (l > lold)
5000       step = -step/10;
5001
5002     p += step;
5003   }
5004
5005   if (i == maxcalls)
5006     return (-2);
5007
5008   fxr = x;
5009
5010
5011   // Search for left x location of fy
5012
5013   p = maxx - 0.5 * params[0];
5014   step = -params[0];
5015   lold = -2.0;
5016   l    = -1.0e300;
5017   i    = 0;
5018
5019   while ((l != lold) && (i < maxcalls)) {
5020     i++;
5021     lold = l;
5022     x    = p + step;
5023     l    = TMath::Abs(LanGauFun(&x,params) - fy);
5024     if (l > lold) {
5025       step = -step / 10.0;
5026     }
5027     p += step;
5028   }
5029
5030   if (i == maxcalls) {
5031     return (-3);
5032   }
5033
5034   fxl  = x;
5035   fwhm = fxr - fxl;
5036
5037   return (0);
5038 }
5039 //_____________________________________________________________________________
5040 Double_t AliTRDCalibraFit::GausConstant(Double_t *x, Double_t *par)
5041 {
5042   //
5043   // Gaus with identical mean
5044   //
5045
5046   Double_t gauss   = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];
5047
5048   return gauss;
5049
5050 }