TRD/AliTRDCalibraFillHisto.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 // AliTRDCalibraFillHisto
  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 fills histos or vectors.
  24 // It can be used for the calibration per chamber but also per group of pads and eventually per pad.
  25 // The user has to choose with the functions SetNz and SetNrphi the precision of the calibration (see AliTRDCalibraMode).
  26 // 2D Histograms (Histo2d) or vectors (Vector2d), then converted in Trees, will be filled
  27 // from RAW DATA in a run or from reconstructed TRD tracks during the offline tracking
  28 // in the function "FollowBackProlongation" (AliTRDtracker)
  29 // Per default the functions to fill are off.
  30 //
  31 // Authors:
  32 //   R. Bailhache (R.Bailhache@gsi.de, rbailhache@ikf.uni-frankfurt.de)
  33 //   J. Book (jbook@ikf.uni-frankfurt.de)
  34 //
  35 //////////////////////////////////////////////////////////////////////////////////////
  36
  37 #include <TProfile2D.h>
  38 #include <TProfile.h>
  39 #include <TFile.h>
  40 #include <TStyle.h>
  41 #include <TCanvas.h>
  42 #include <TObjArray.h>
  43 #include <TObject.h>
  44 #include <TH1F.h>
  45 #include <TH2I.h>
  46 #include <TH2.h>
  47 #include <TStopwatch.h>
  48 #include <TMath.h>
  49 #include <TDirectory.h>
  50 #include <TTreeStream.h>
  51 #include <TVectorD.h>
  52 #include <TLinearFitter.h>
  53
  54 #include "AliLog.h"
  55
  56 #include "AliESDtrack.h"
  57 #include "AliTRDCalibraFillHisto.h"
  58 #include "AliTRDCalibraMode.h"
  59 #include "AliTRDCalibraVector.h"
  60 #include "AliTRDCalibraVdriftLinearFit.h"
  61 #include "AliTRDCalibraExbAltFit.h"
  62 #include "AliTRDcalibDB.h"
  63 #include "AliTRDCommonParam.h"
  64 #include "AliTRDpadPlane.h"
  65 #include "AliTRDcluster.h"
  66 #include "AliTRDtrackV1.h"
  67 #include "AliRawReader.h"
  68 #include "AliRawReaderDate.h"
  69 #include "AliTRDgeometry.h"
  70 #include "AliTRDfeeParam.h"
  71 #include "./Cal/AliTRDCalROC.h"
  72 #include "./Cal/AliTRDCalPad.h"
  73 #include "./Cal/AliTRDCalDet.h"
  74
  75 #include "AliTRDdigitsManager.h"
  76 #include "AliTRDdigitsParam.h"
  77 #include "AliTRDSignalIndex.h"
  78 #include "AliTRDarrayADC.h"
  79
  80 #include "AliTRDrawStream.h"
  81
  82 #include "AliCDBEntry.h"
  83 #include "AliCDBManager.h"
  84
  85 #ifdef ALI_DATE
  86 #include "event.h"
  87 #endif
  88
  89
  90 ClassImp(AliTRDCalibraFillHisto)
  91
  92 AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
  93 Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
  94
  95 //_____________singleton implementation_________________________________________________
  96 AliTRDCalibraFillHisto *AliTRDCalibraFillHisto::Instance()
  97 {
  98   //
  99   // Singleton implementation
 100   //
 101
 102   if (fgTerminated != kFALSE) {
 103     return 0;
 104   }
 105
 106   if (fgInstance == 0) {
 107     fgInstance = new AliTRDCalibraFillHisto();
 108   }
 109
 110   return fgInstance;
 111
 112 }
 113
 114 //______________________________________________________________________________________
 115 void AliTRDCalibraFillHisto::Terminate()
 116 {
 117   //
 118   // Singleton implementation
 119   // Deletes the instance of this class
 120   //
 121
 122   fgTerminated = kTRUE;
 123
 124   if (fgInstance != 0) {
 125     delete fgInstance;
 126     fgInstance = 0;
 127   }
 128
 129 }
 130
 131 //______________________________________________________________________________________
 132 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
 133   :TObject()
 134   ,fGeo(0)
 135   ,fCalibDB(0)
 136   ,fIsHLT(kFALSE)
 137   ,fCH2dOn(kFALSE)
 138   ,fPH2dOn(kFALSE)
 139   ,fPRF2dOn(kFALSE)
 140   ,fHisto2d(kFALSE)
 141   ,fVector2d(kFALSE)
 142   ,fLinearFitterOn(kFALSE)
 143   ,fLinearFitterDebugOn(kFALSE)
 144   ,fExbAltFitOn(kFALSE)
 145   ,fScaleWithTPCSignal(kFALSE)
 146   ,fRelativeScale(0)
 147   ,fThresholdClusterPRF2(15.0)
 148   ,fLimitChargeIntegration(kFALSE)
 149   ,fFillWithZero(kFALSE)
 150   ,fNormalizeNbOfCluster(kFALSE)
 151   ,fMaxCluster(0)
 152   ,fNbMaxCluster(0)
 153   ,fCutWithVdriftCalib(kFALSE)
 154   ,fMinNbTRDtracklets(0)
 155   ,fMinTRDMomentum(0.0)
 156   ,fTakeSnapshot(kTRUE)
 157   ,fFirstRunGain(0)
 158   ,fVersionGainUsed(0)
 159   ,fSubVersionGainUsed(0)
 160   ,fFirstRunGainLocal(0)
 161   ,fVersionGainLocalUsed(0)
 162   ,fSubVersionGainLocalUsed(0)
 163   ,fFirstRunVdrift(0)
 164   ,fVersionVdriftUsed(0)
 165   ,fSubVersionVdriftUsed(0)
 166   ,fFirstRunExB(0)
 167   ,fVersionExBUsed(0)
 168   ,fSubVersionExBUsed(0)
 169   ,fCalibraMode(new AliTRDCalibraMode())
 170   ,fDebugStreamer(0)
 171   ,fDebugLevel(0)
 172   ,fDetectorPreviousTrack(-1)
 173   ,fMCMPrevious(-1)
 174   ,fROBPrevious(-1)
 175   ,fNumberClusters(1)
 176   ,fNumberClustersf(30)
 177   ,fNumberClustersProcent(0.5)
 178   ,fThresholdClustersDAQ(120.0)
 179   ,fNumberRowDAQ(2)
 180   ,fNumberColDAQ(4)
 181   ,fProcent(6.0)
 182   ,fDifference(17)
 183   ,fNumberTrack(0)
 184   ,fTimeMax(0)
 185   ,fSf(10.0)
 186   ,fRangeHistoCharge(150)
 187   ,fNumberBinCharge(50)
 188   ,fNumberBinPRF(10)
 189   ,fNgroupprf(3)
 190   ,fAmpTotal(0x0)
 191   ,fPHPlace(0x0)
 192   ,fPHValue(0x0)
 193   ,fGoodTracklet(kTRUE)
 194   ,fLinearFitterTracklet(0x0)
 195   ,fEntriesCH(0x0)
 196   ,fEntriesLinearFitter(0x0)
 197   ,fCalibraVector(0x0)
 198   ,fPH2d(0x0)
 199   ,fPRF2d(0x0)
 200   ,fCH2d(0x0)
 201   ,fLinearFitterArray(540)
 202   ,fLinearVdriftFit(0x0)
 203   ,fExbAltFit(0x0)
 204   ,fCalDetGain(0x0)
 205   ,fCalROCGain(0x0)
 206 {
 207   //
 208   // Default constructor
 209   //
 210
 211   //
 212   // Init some default values
 213   //
 214
 215   fNumberUsedCh[0]       = 0;
 216   fNumberUsedCh[1]       = 0;
 217   fNumberUsedPh[0]       = 0;
 218   fNumberUsedPh[1]       = 0;
 219
 220   fGeo = new AliTRDgeometry();
 221   fCalibDB = AliTRDcalibDB::Instance();
 222 }
 223
 224 //______________________________________________________________________________________
 225 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
 226   :TObject(c)
 227   ,fGeo(0)
 228   ,fCalibDB(0)
 229   ,fIsHLT(c.fIsHLT)
 230   ,fCH2dOn(c.fCH2dOn)
 231   ,fPH2dOn(c.fPH2dOn)
 232   ,fPRF2dOn(c.fPRF2dOn)
 233   ,fHisto2d(c.fHisto2d)
 234   ,fVector2d(c.fVector2d)
 235   ,fLinearFitterOn(c.fLinearFitterOn)
 236   ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
 237   ,fExbAltFitOn(c.fExbAltFitOn)
 238   ,fScaleWithTPCSignal(c.fScaleWithTPCSignal)
 239   ,fRelativeScale(c.fRelativeScale)
 240   ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
 241   ,fLimitChargeIntegration(c.fLimitChargeIntegration)
 242   ,fFillWithZero(c.fFillWithZero)
 243   ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster)
 244   ,fMaxCluster(c.fMaxCluster)
 245   ,fNbMaxCluster(c.fNbMaxCluster)
 246   ,fCutWithVdriftCalib(c.fCutWithVdriftCalib)
 247   ,fMinNbTRDtracklets(c.fMinNbTRDtracklets)
 248   ,fMinTRDMomentum(c.fMinTRDMomentum)
 249   ,fTakeSnapshot(c.fTakeSnapshot)
 250   ,fFirstRunGain(c.fFirstRunGain)
 251   ,fVersionGainUsed(c.fVersionGainUsed)
 252   ,fSubVersionGainUsed(c.fSubVersionGainUsed)
 253   ,fFirstRunGainLocal(c.fFirstRunGainLocal)
 254   ,fVersionGainLocalUsed(c.fVersionGainLocalUsed)
 255   ,fSubVersionGainLocalUsed(c.fSubVersionGainLocalUsed)
 256   ,fFirstRunVdrift(c.fFirstRunVdrift)
 257   ,fVersionVdriftUsed(c.fVersionVdriftUsed)
 258   ,fSubVersionVdriftUsed(c.fSubVersionVdriftUsed)
 259   ,fFirstRunExB(c.fFirstRunExB)
 260   ,fVersionExBUsed(c.fVersionExBUsed)
 261   ,fSubVersionExBUsed(c.fSubVersionExBUsed)
 262   ,fCalibraMode(0x0)
 263   ,fDebugStreamer(0)
 264   ,fDebugLevel(c.fDebugLevel)
 265   ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
 266   ,fMCMPrevious(c.fMCMPrevious)
 267   ,fROBPrevious(c.fROBPrevious)
 268   ,fNumberClusters(c.fNumberClusters)
 269   ,fNumberClustersf(c.fNumberClustersf)
 270   ,fNumberClustersProcent(c.fNumberClustersProcent)
 271   ,fThresholdClustersDAQ(c.fThresholdClustersDAQ)
 272   ,fNumberRowDAQ(c.fNumberRowDAQ)
 273   ,fNumberColDAQ(c.fNumberColDAQ)
 274   ,fProcent(c.fProcent)
 275   ,fDifference(c.fDifference)
 276   ,fNumberTrack(c.fNumberTrack)
 277   ,fTimeMax(c.fTimeMax)
 278   ,fSf(c.fSf)
 279   ,fRangeHistoCharge(c.fRangeHistoCharge)
 280   ,fNumberBinCharge(c.fNumberBinCharge)
 281   ,fNumberBinPRF(c.fNumberBinPRF)
 282   ,fNgroupprf(c.fNgroupprf)
 283   ,fAmpTotal(0x0)
 284   ,fPHPlace(0x0)
 285   ,fPHValue(0x0)
 286   ,fGoodTracklet(c.fGoodTracklet)
 287   ,fLinearFitterTracklet(0x0)
 288   ,fEntriesCH(0x0)
 289   ,fEntriesLinearFitter(0x0)
 290   ,fCalibraVector(0x0)
 291   ,fPH2d(0x0)
 292   ,fPRF2d(0x0)
 293   ,fCH2d(0x0)
 294   ,fLinearFitterArray(540)
 295   ,fLinearVdriftFit(0x0)
 296   ,fExbAltFit(0x0)
 297   ,fCalDetGain(0x0)
 298   ,fCalROCGain(0x0)
 299 {
 300   //
 301   // Copy constructor
 302   //
 303   if(c.fCalibraMode)   fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
 304   if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
 305   if(c.fPH2d) {
 306     fPH2d = new TProfile2D(*c.fPH2d);
 307     fPH2d->SetDirectory(0);
 308   }
 309   if(c.fPRF2d) {
 310     fPRF2d = new TProfile2D(*c.fPRF2d);
 311     fPRF2d->SetDirectory(0);
 312   }
 313   if(c.fCH2d) {
 314     fCH2d = new TH2I(*c.fCH2d);
 315     fCH2d->SetDirectory(0);
 316   }
 317   if(c.fLinearVdriftFit){
 318     fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
 319   }
 320   if(c.fExbAltFit){
 321     fExbAltFit = new AliTRDCalibraExbAltFit(*c.fExbAltFit);
 322   }
 323
 324   if(c.fCalDetGain)  fCalDetGain   = new AliTRDCalDet(*c.fCalDetGain);
 325   if(c.fCalROCGain)  fCalROCGain   = new AliTRDCalROC(*c.fCalROCGain);
 326
 327   if (fGeo) {
 328     delete fGeo;
 329   }
 330   fGeo = new AliTRDgeometry();
 331   fCalibDB = AliTRDcalibDB::Instance();
 332
 333   fNumberUsedCh[0]       = 0;
 334   fNumberUsedCh[1]       = 0;
 335   fNumberUsedPh[0]       = 0;
 336   fNumberUsedPh[1]       = 0;
 337
 338 }
 339
 340 //____________________________________________________________________________________
 341 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
 342 {
 343   //
 344   // AliTRDCalibraFillHisto destructor
 345   //
 346
 347   ClearHistos();
 348   if ( fDebugStreamer ) delete fDebugStreamer;
 349
 350   if ( fCalDetGain )  delete fCalDetGain;
 351   if ( fCalROCGain )  delete fCalROCGain;
 352
 353   if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
 354
 355   delete [] fPHPlace;
 356   delete [] fPHValue;
 357   delete [] fEntriesCH;
 358   delete [] fEntriesLinearFitter;
 359   delete [] fAmpTotal;
 360
 361   for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
 362     TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
 363     if(f) { delete f;}
 364   }
 365   if(fLinearVdriftFit) delete fLinearVdriftFit;
 366   if(fExbAltFit) delete fExbAltFit;
 367   if (fGeo) {
 368     delete fGeo;
 369   }
 370
 371 }
 372 //_____________________________________________________________________________
 373 void AliTRDCalibraFillHisto::Destroy()
 374 {
 375   //
 376   // Delete instance
 377   //
 378
 379   if (fgInstance) {
 380     delete fgInstance;
 381     fgInstance = 0x0;
 382   }
 383 }
 384 //_____________________________________________________________________________
 385 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
 386 {
 387   //
 388   // Delete DebugStreamer
 389   //
 390
 391   if ( fDebugStreamer ) delete fDebugStreamer;
 392   fDebugStreamer = 0x0;
 393
 394 }
 395 //_____________________________________________________________________________
 396 void AliTRDCalibraFillHisto::ClearHistos()
 397 {
 398   //
 399   // Delete the histos
 400   //
 401
 402   if (fPH2d) {
 403     delete fPH2d;
 404     fPH2d  = 0x0;
 405   }
 406   if (fCH2d) {
 407     delete fCH2d;
 408     fCH2d  = 0x0;
 409   }
 410   if (fPRF2d) {
 411     delete fPRF2d;
 412     fPRF2d = 0x0;
 413   }
 414
 415 }
 416 //////////////////////////////////////////////////////////////////////////////////
 417 // calibration with AliTRDtrackV1: Init, Update
 418 //////////////////////////////////////////////////////////////////////////////////
 419 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
 420 Bool_t AliTRDCalibraFillHisto::Init2Dhistos(Int_t nboftimebin)
 421 {
 422   //
 423   // Init the histograms and stuff to be filled
 424   //
 425
 426   // DB Setting
 427   // Get cal
 428   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 429   if (!cal) {
 430     AliInfo("Could not get calibDB");
 431     return kFALSE;
 432   }
 433   AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
 434   if (!parCom) {
 435     AliInfo("Could not get CommonParam");
 436     return kFALSE;
 437   }
 438
 439   // Some parameters
 440   if(nboftimebin > 0) fTimeMax = nboftimebin;
 441   else fTimeMax = cal->GetNumberOfTimeBinsDCS();
 442   if(fTimeMax <= 0) fTimeMax = 30;
 443   printf("////////////////////////////////////////////\n");
 444   printf("Number of time bins in calibration component %d\n",fTimeMax);
 445   printf("////////////////////////////////////////////\n");
 446   fSf                 = parCom->GetSamplingFrequency();
 447   if(!fNormalizeNbOfCluster) fRelativeScale = 20.0;
 448   else fRelativeScale = 1.18;
 449   fNumberClustersf    = fTimeMax;
 450   fNumberClusters     = (Int_t)(fNumberClustersProcent*fTimeMax);
 451
 452   // Init linear fitter
 453   if(!fLinearFitterTracklet) {
 454     fLinearFitterTracklet = new TLinearFitter(2,"pol1");
 455     fLinearFitterTracklet->StoreData(kTRUE);
 456   }
 457
 458   // Calcul Xbins Chambd0, Chamb2
 459   Int_t ntotal0 = CalculateTotalNumberOfBins(0);
 460   Int_t ntotal1 = CalculateTotalNumberOfBins(1);
 461   Int_t ntotal2 = CalculateTotalNumberOfBins(2);
 462
 463   // If vector method On initialised all the stuff
 464   if(fVector2d){
 465     fCalibraVector = new AliTRDCalibraVector();
 466     fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
 467     fCalibraVector->SetTimeMax(fTimeMax);
 468     if(fNgroupprf != 0) {
 469       fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
 470       fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
 471     }
 472     else {
 473       fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
 474       fCalibraVector->SetPRFRange(1.5);
 475     }
 476     for(Int_t k = 0; k < 3; k++){
 477       fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
 478       fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
 479     }
 480     fCalibraVector->SetNzNrphi(0,fCalibraMode->GetNz(0),fCalibraMode->GetNrphi(0));
 481     fCalibraVector->SetNzNrphi(1,fCalibraMode->GetNz(1),fCalibraMode->GetNrphi(1));
 482     fCalibraVector->SetNzNrphi(2,fCalibraMode->GetNz(2),fCalibraMode->GetNrphi(2));
 483     fCalibraVector->SetNbGroupPRF(fNgroupprf);
 484   }
 485
 486   // Create the 2D histos corresponding to the pad groupCalibration mode
 487   if (fCH2dOn) {
 488
 489     AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
 490                 ,fCalibraMode->GetNz(0)
 491                 ,fCalibraMode->GetNrphi(0)));
 492
 493     // Create the 2D histo
 494     if (fHisto2d) {
 495       CreateCH2d(ntotal0);
 496     }
 497     // Variable
 498     fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
 499     for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
 500       fAmpTotal[k] = 0.0;
 501     }
 502     //Statistics
 503     fEntriesCH = new Int_t[ntotal0];
 504     for(Int_t k = 0; k < ntotal0; k++){
 505       fEntriesCH[k] = 0;
 506     }
 507
 508   }
 509   if (fPH2dOn) {
 510
 511     AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
 512                 ,fCalibraMode->GetNz(1)
 513                 ,fCalibraMode->GetNrphi(1)));
 514
 515     // Create the 2D histo
 516     if (fHisto2d) {
 517       CreatePH2d(ntotal1);
 518     }
 519     // Variable
 520     fPHPlace = new Short_t[fTimeMax];
 521     for (Int_t k = 0; k < fTimeMax; k++) {
 522       fPHPlace[k] = -1;
 523     }
 524     fPHValue = new Float_t[fTimeMax];
 525     for (Int_t k = 0; k < fTimeMax; k++) {
 526       fPHValue[k] = 0.0;
 527     }
 528   }
 529   if (fLinearFitterOn) {
 530     if(fLinearFitterDebugOn) {
 531       fLinearFitterArray.SetName("ArrayLinearFitters");
 532       fEntriesLinearFitter = new Int_t[540];
 533       for(Int_t k = 0; k < 540; k++){
 534         fEntriesLinearFitter[k] = 0;
 535       }
 536     }
 537     fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
 538     TString nameee("Ver");
 539     nameee += fVersionExBUsed;
 540     nameee += "Subver";
 541     nameee += fSubVersionExBUsed;
 542     nameee += "FirstRun";
 543     nameee += fFirstRunExB;
 544     nameee += "Nz";
 545     fLinearVdriftFit->SetNameCalibUsed(nameee);
 546   }
 547   if(fExbAltFitOn){
 548     fExbAltFit = new AliTRDCalibraExbAltFit();
 549   }
 550
 551   if (fPRF2dOn) {
 552
 553     AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
 554                 ,fCalibraMode->GetNz(2)
 555                 ,fCalibraMode->GetNrphi(2)));
 556     // Create the 2D histo
 557     if (fHisto2d) {
 558       CreatePRF2d(ntotal2);
 559     }
 560   }
 561
 562   return kTRUE;
 563
 564 }
 565 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
 566 Bool_t AliTRDCalibraFillHisto::InitCalDet()
 567 {
 568   //
 569   // Init the Gain Cal Det
 570   //
 571
 572   // DB Setting
 573   // Get cal
 574   AliCDBEntry *entry = 0x0;
 575   if(fTakeSnapshot) {
 576     entry = AliCDBManager::Instance()->Get("TRD/Calib/ChamberGainFactor");
 577   }
 578   else {
 579     entry = AliCDBManager::Instance()->Get("TRD/Calib/ChamberGainFactor",fFirstRunGain,fVersionGainUsed,fSubVersionGainUsed);
 580   }
 581   if(!entry) {
 582     AliError("No gain det calibration entry found");
 583     return kFALSE;
 584   }
 585   AliTRDCalDet *calDet = (AliTRDCalDet *)entry->GetObject();
 586   if(!calDet) {
 587     AliError("No calDet gain found");
 588     return kFALSE;
 589   }
 590
 591
 592   if( fCalDetGain ){
 593     fCalDetGain->~AliTRDCalDet();
 594     new(fCalDetGain) AliTRDCalDet(*(calDet));
 595   }else fCalDetGain = new AliTRDCalDet(*(calDet));
 596
 597
 598   // title CH2d
 599   TString name("Ver");
 600   name += fVersionGainUsed;
 601   name += "Subver";
 602   name += fSubVersionGainUsed;
 603   name += "FirstRun";
 604   name += fFirstRunGain;
 605   name += "Nz";
 606   name += fCalibraMode->GetNz(0);
 607   name += "Nrphi";
 608   name += fCalibraMode->GetNrphi(0);
 609
 610   fCH2d->SetTitle(name);
 611
 612   // title PH2d
 613   TString namee("Ver");
 614   namee += fVersionVdriftUsed;
 615   namee += "Subver";
 616   namee += fSubVersionVdriftUsed;
 617   namee += "FirstRun";
 618   namee += fFirstRunVdrift;
 619   namee += "Nz";
 620   namee += fCalibraMode->GetNz(1);
 621   namee += "Nrphi";
 622   namee += fCalibraMode->GetNrphi(1);
 623
 624   fPH2d->SetTitle(namee);
 625
 626   // title AliTRDCalibraVdriftLinearFit
 627   TString nameee("Ver");
 628   nameee += fVersionExBUsed;
 629   nameee += "Subver";
 630   nameee += fSubVersionExBUsed;
 631   nameee += "FirstRun";
 632   nameee += fFirstRunExB;
 633   nameee += "Nz";
 634
 635
 636   fLinearVdriftFit->SetNameCalibUsed(nameee);
 637
 638
 639
 640   return kTRUE;
 641
 642 }
 643 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
 644 Bool_t AliTRDCalibraFillHisto::InitCalPad(Int_t detector)
 645 {
 646   //
 647   // Init the Gain Cal Pad
 648   //
 649
 650   // DB Setting
 651   // Get cal
 652   AliCDBEntry *entry = 0x0;
 653   if(fTakeSnapshot) {
 654     entry = AliCDBManager::Instance()->Get("TRD/Calib/LocalGainFactor");
 655   }
 656   else {
 657     entry = AliCDBManager::Instance()->Get("TRD/Calib/LocalGainFactor",fFirstRunGain,fVersionGainUsed,fSubVersionGainUsed);
 658   }
 659   if(!entry) {
 660     AliError("No gain pad calibration entry found");
 661     return kFALSE;
 662   }
 663   AliTRDCalPad *calPad = (AliTRDCalPad *)entry->GetObject();
 664   if(!calPad) {
 665     AliError("No calPad gain found");
 666     return kFALSE;
 667   }
 668   AliTRDCalROC *calRoc = (AliTRDCalROC *)calPad->GetCalROC(detector);
 669   if(!calRoc) {
 670     AliError("No calRoc gain found");
 671     return kFALSE;
 672   }
 673
 674   if( fCalROCGain ){
 675     fCalROCGain->~AliTRDCalROC();
 676     new(fCalROCGain) AliTRDCalROC(*(calRoc));
 677   }else fCalROCGain = new AliTRDCalROC(*(calRoc));
 678
 679
 680
 681
 682
 683   return kTRUE;
 684
 685 }
 686 //____________Offline tracking in the AliTRDtracker____________________________
 687 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(const AliTRDtrackV1 *t,const AliESDtrack *esdtrack)
 688 {
 689   //
 690   // Use AliTRDtrackV1 for the calibration
 691   //
 692
 693
 694   const AliTRDseedV1 *tracklet = 0x0;          // tracklet per plane
 695   AliTRDcluster *cl      = 0x0;                // cluster attached now to the tracklet
 696   AliTRDcluster *cls     = 0x0;                // shared cluster attached now to the tracklet
 697   Bool_t         newtr   = kTRUE;              // new track
 698
 699
 700   //
 701   // Cut on the number of TRD tracklets
 702   //
 703   Int_t numberoftrdtracklets = t->GetNumberOfTracklets();
 704   if(numberoftrdtracklets < fMinNbTRDtracklets) return kFALSE;
 705
 706   Double_t tpcsignal = 1.0;
 707   if(esdtrack) tpcsignal = esdtrack->GetTPCsignal()/50.0;
 708   if(fScaleWithTPCSignal && tpcsignal <0.00001) return kFALSE;
 709
 710   //
 711   if (!fCalibDB) {
 712     AliInfo("Could not get calibDB");
 713     return kFALSE;
 714   }
 715
 716
 717   ///////////////////////////
 718   // loop over the tracklet
 719   ///////////////////////////
 720   for(Int_t itr = 0; itr < 6; itr++){
 721
 722     if(!(tracklet = t->GetTracklet(itr))) continue;
 723     if(!tracklet->IsOK()) continue;
 724     fNumberTrack++;
 725     ResetfVariablestracklet();
 726     Float_t momentum = t->GetMomentum(itr);
 727     if(TMath::Abs(momentum) < fMinTRDMomentum)  continue;
 728
 729
 730     //////////////////////////////////////////
 731     // localisation of the tracklet and dqdl
 732     //////////////////////////////////////////
 733     Int_t layer    = tracklet->GetPlane();
 734     Int_t ic = 0;
 735     while(!(cl = tracklet->GetClusters(ic++))) continue;
 736     Int_t detector = cl->GetDetector();
 737     if (detector != fDetectorPreviousTrack) {
 738       // if not a new track
 739       if(!newtr){
 740         // don't use the rest of this track if in the same plane
 741         if (layer == GetLayer(fDetectorPreviousTrack)) {
 742           //printf("bad tracklet, same layer for detector %d\n",detector);
 743           break;
 744         }
 745       }
 746       //Localise the detector bin
 747       LocalisationDetectorXbins(detector);
 748       // Get calib objects
 749       if(!fIsHLT) InitCalPad(detector);
 750
 751       // reset
 752       fDetectorPreviousTrack = detector;
 753     }
 754     newtr = kFALSE;
 755
 756     ////////////////////////////
 757     // loop over the clusters
 758     ////////////////////////////
 759     Double_t chargeQ = 0.0;
 760     Int_t nbclusters = 0;
 761     for(int jc=0; jc<AliTRDseedV1::kNtb; jc++){
 762       if(!(cl = tracklet->GetClusters(jc))) continue;
 763       nbclusters++;
 764
 765       // Store the info bis of the tracklet
 766       Int_t row = cl->GetPadRow();
 767       Int_t col = cl->GetPadCol();
 768       CheckGoodTrackletV1(cl);
 769       Int_t     group[2] = {0,0};
 770       if(fCH2dOn)  group[0]  = CalculateCalibrationGroup(0,row,col);
 771       if(fPH2dOn)  group[1]  = CalculateCalibrationGroup(1,row,col);
 772       // Add the charge if shared cluster
 773       cls = tracklet->GetClusters(jc+AliTRDseedV1::kNtb);
 774       //
 775       //Scale with TPC signal or not
 776       if(!fScaleWithTPCSignal) {
 777         chargeQ += StoreInfoCHPHtrack(cl, tracklet->GetQperTB(jc),group,row,col,cls); //tracklet->GetdQdl(jc)
 778         //printf("Do not scale now\n");
 779       }
 780       else {
 781         chargeQ += StoreInfoCHPHtrack(cl, tracklet->GetQperTB(jc)/tpcsignal,group,row,col,cls); //tracklet->GetdQdl(jc)
 782       }
 783
 784     }
 785
 786     ////////////////////////////////////////
 787     // Fill the stuffs if a good tracklet
 788     ////////////////////////////////////////
 789     if (fGoodTracklet) {
 790
 791       // drift velocity unables to cut bad tracklets
 792       Bool_t  pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
 793
 794       //printf("pass %d and nbclusters %d\n",pass,nbclusters);
 795
 796       // Gain calibration
 797       if (fCH2dOn) {
 798         if(fCutWithVdriftCalib) {
 799           if(pass) FillTheInfoOfTheTrackCH(nbclusters);
 800         } else {
 801           FillTheInfoOfTheTrackCH(nbclusters);
 802         }
 803       }
 804
 805       // PH calibration
 806       if (fPH2dOn) {
 807         if(fCutWithVdriftCalib) {
 808           if(pass) FillTheInfoOfTheTrackPH();
 809         }
 810         else {
 811           FillTheInfoOfTheTrackPH();
 812         }
 813       }
 814
 815       if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
 816
 817
 818       /////////////////////////////////////////////////////////
 819       // Debug
 820       ////////////////////////////////////////////////////////
 821       if(fDebugLevel > 0){
 822         //printf("test\n");
 823         if ( !fDebugStreamer ) {
 824           //debug stream
 825           TDirectory *backup = gDirectory;
 826           fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
 827           if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 828         }
 829
 830         Int_t    stacke = AliTRDgeometry::GetStack(detector);
 831         Int_t    sme    = AliTRDgeometry::GetSector(detector);
 832         Int_t    layere  = AliTRDgeometry::GetLayer(detector);
 833         // Some variables
 834         Float_t b[2] = {0.0,0.0};
 835         Float_t bCov[3] = {0.0,0.0,0.0};
 836         if(esdtrack) esdtrack->GetImpactParameters(b,bCov);
 837         if (bCov[0]<=0 || bCov[2]<=0) {
 838           bCov[0]=0; bCov[2]=0;
 839         }
 840         Float_t dcaxy = b[0];
 841         Float_t dcaz = b[1];
 842         Int_t tpcnbclusters = 0;
 843         if(esdtrack) tpcnbclusters = esdtrack->GetTPCclusters(0);
 844         Double_t ttpcsignal = 0.0;
 845         if(esdtrack) ttpcsignal = esdtrack->GetTPCsignal();
 846         Int_t    cutvdriftlinear = 0;
 847         if(!pass) cutvdriftlinear = 1;
 848
 849         (* fDebugStreamer) << "FillCharge"<<
 850           "detector="<<detector<<
 851           "stack="<<stacke<<
 852           "sm="<<sme<<
 853           "layere="<<layere<<
 854           "dcaxy="<<dcaxy<<
 855           "dcaz="<<dcaz<<
 856           "nbtpccls="<<tpcnbclusters<<
 857           "tpcsignal="<<ttpcsignal<<
 858           "cutvdriftlinear="<<cutvdriftlinear<<
 859           "ptrd="<<momentum<<
 860           "nbtrdtracklet="<<numberoftrdtracklets<<
 861           "charge="<<chargeQ<<
 862           "\n";
 863       }
 864
 865
 866     } // if a good tracklet
 867   }
 868
 869   return kTRUE;
 870
 871 }
 872 ///////////////////////////////////////////////////////////////////////////////////
 873 // Routine inside the update with AliTRDtrack
 874 ///////////////////////////////////////////////////////////////////////////////////
 875 //____________Offine tracking in the AliTRDtracker_____________________________
 876 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
 877 {
 878   //
 879   // Drift velocity calibration:
 880   // Fit the clusters with a straight line
 881   // From the slope find the drift velocity
 882   //
 883
 884   ////////////////////////////////////////////////
 885   //Number of points: if less than 3 return kFALSE
 886   /////////////////////////////////////////////////
 887   if(nbclusters <= 2) return kFALSE;
 888
 889   ////////////
 890   //Variables
 891   ////////////
 892   // results of the linear fit
 893   Double_t dydt                       = 0.0;                                // dydt tracklet after straight line fit
 894   Double_t errorpar                   = 0.0;                                // error after straight line fit on dy/dt
 895   Double_t pointError                 = 0.0;                                // error after straight line fit
 896   // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
 897   Int_t    crossrow                   = 0;                                  // if it crosses a pad row
 898   Int_t    rowp                       = -1;                                 // if it crosses a pad row
 899   Float_t  tnt                        = tracklet->GetTilt();                // tan tiltingangle
 900   fLinearFitterTracklet->ClearPoints();
 901
 902
 903   ///////////////////////////////////////////
 904   // Take the parameters of the track
 905   //////////////////////////////////////////
 906   // take now the snp, tnp and tgl from the track
 907   Double_t snp = tracklet->GetSnp();             // sin dy/dx at the end of the chamber
 908   Double_t tnp = 0.0;                            // dy/dx at the end of the chamber
 909   if( TMath::Abs(snp) <  1.){
 910     tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
 911   }
 912   Double_t tgl  = tracklet->GetTgl();           // dz/dl
 913   Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp);   // dz/dx calculated from dz/dl
 914   // at the entrance
 915   //Double_t tnp = tracklet->GetYref(1);      // dy/dx at the entrance of the chamber
 916   //Double_t tgl = tracklet->GetZref(1);      // dz/dl at the entrance of the chamber
 917   //Double_t dzdx = tgl;                      //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
 918   // at the end with correction due to linear fit
 919   //Double_t tnp = tracklet->GetYfit(1);      // dy/dx at the end of the chamber after fit correction
 920   //Double_t tgl = tracklet->GetZfit(1);      // dz/dl at the end of the chamber after fit correction
 921
 922
 923   ////////////////////////////
 924   // loop over the clusters
 925   ////////////////////////////
 926   Int_t  nbli = 0;
 927   AliTRDcluster *cl                   = 0x0;
 928   //////////////////////////////
 929   // Check no shared clusters
 930   //////////////////////////////
 931   for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
 932     cl = tracklet->GetClusters(icc);
 933     if(cl)  crossrow = 1;
 934   }
 935   //////////////////////////////////
 936   // Loop clusters
 937   //////////////////////////////////
 938
 939   Float_t sigArr[AliTRDfeeParam::GetNcol()];
 940   memset(sigArr, 0, AliTRDfeeParam::GetNcol()*sizeof(sigArr[0]));
 941   Int_t ncl=0, tbf=0, tbl=0;
 942
 943   for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
 944     if(!(cl = tracklet->GetClusters(ic))) continue;
 945
 946     if(!tbf) tbf=ic;
 947     tbl=ic;
 948     ncl++;
 949     Int_t col = cl->GetPadCol();
 950     for(int ip=-1, jp=2; jp<5; ip++, jp++){
 951       Int_t idx=col+ip;
 952       if(idx<0 || idx>=AliTRDfeeParam::GetNcol()) continue;
 953       sigArr[idx]+=((Float_t)cl->GetSignals()[jp]);
 954     }
 955
 956     if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
 957
 958     Double_t ycluster                 = cl->GetY();
 959     Int_t time                        = cl->GetPadTime();
 960     Double_t timeis                   = time/fSf;
 961     //See if cross two pad rows
 962     Int_t    row                      = cl->GetPadRow();
 963     if(rowp==-1) rowp                 = row;
 964     if(row != rowp) crossrow          = 1;
 965
 966     fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
 967     nbli++;
 968
 969
 970   }
 971
 972   ////////////////////////////////////
 973   // Do the straight line fit now
 974   ///////////////////////////////////
 975   if(nbli <= 2){
 976     fLinearFitterTracklet->ClearPoints();
 977     return kFALSE;
 978   }
 979   TVectorD pars;
 980   fLinearFitterTracklet->Eval();
 981   fLinearFitterTracklet->GetParameters(pars);
 982   pointError  =  TMath::Sqrt(TMath::Abs(fLinearFitterTracklet->GetChisquare()/(nbli-2)));
 983   errorpar    =  fLinearFitterTracklet->GetParError(1)*pointError;
 984   dydt        = pars[1];
 985   //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
 986   fLinearFitterTracklet->ClearPoints();
 987
 988   ////////////////////////////////////
 989   // Calc the projection of the clusters on the y direction
 990   ///////////////////////////////////
 991
 992   Float_t signalSum(0.);
 993   Float_t mean = 0.0, rms = 0.0;
 994   Float_t dydx = tracklet->GetYref(1), tilt = tracklet->GetTilt(); // ,dzdx = tracklet->GetZref(1); (identical to the previous definition!)
 995   Float_t dz = dzdx*(tbl-tbf)/10;
 996   if(ncl>10){
 997     for(Int_t ip(0); ip<AliTRDfeeParam::GetNcol(); ip++){
 998       signalSum+=sigArr[ip];
 999       mean+=ip*sigArr[ip];
1000     }
1001     if(signalSum > 0.0) mean/=signalSum;
1002
1003     for(Int_t ip = 0; ip<AliTRDfeeParam::GetNcol(); ip++)
1004       rms+=sigArr[ip]*(ip-mean)*(ip-mean);
1005
1006     if(signalSum > 0.0) rms = TMath::Sqrt(TMath::Abs(rms/signalSum));
1007
1008     rms -= TMath::Abs(dz*tilt);
1009     dydx -= dzdx*tilt;
1010   }
1011
1012   ////////////////////////////////
1013   // Debug stuff
1014   ///////////////////////////////
1015
1016
1017   if(fDebugLevel > 1){
1018     if ( !fDebugStreamer ) {
1019       //debug stream
1020       TDirectory *backup = gDirectory;
1021       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1022       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1023     }
1024
1025     float xcoord = tnp-dzdx*tnt;
1026     float pt = tracklet->GetPt();
1027     Int_t layer = GetLayer(fDetectorPreviousTrack);
1028
1029     (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
1030       //"snpright="<<snpright<<
1031       "nbli="<<nbli<<
1032       "nbclusters="<<nbclusters<<
1033       "detector="<<fDetectorPreviousTrack<<
1034       "layer="<<layer<<
1035       "snp="<<snp<<
1036       "tnp="<<tnp<<
1037       "tgl="<<tgl<<
1038       "tnt="<<tnt<<
1039       "dydt="<<dydt<<
1040       "dzdx="<<dzdx<<
1041       "crossrow="<<crossrow<<
1042       "errorpar="<<errorpar<<
1043       "pointError="<<pointError<<
1044       "xcoord="<<xcoord<<
1045       "pt="<<pt<<
1046       "rms="<<rms<<
1047       "dydx="<<dydx<<
1048       "dz="<<dz<<
1049       "tilt="<<tilt<<
1050       "ncl="<<ncl<<
1051       "\n";
1052
1053   }
1054
1055   /////////////////////////
1056   // Cuts quality
1057   ////////////////////////
1058
1059   if(nbclusters < fNumberClusters) return kFALSE;
1060   if(nbclusters > fNumberClustersf) return kFALSE;
1061   if(pointError >= 0.3) return kFALSE;
1062   if(crossrow == 1) return kTRUE;
1063
1064   ///////////////////////
1065   // Fill
1066   //////////////////////
1067
1068   if(fLinearFitterOn){
1069     //Add to the linear fitter of the detector
1070     if( TMath::Abs(snp) <  1.){
1071       Double_t x = tnp-dzdx*tnt;
1072       if(fLinearFitterDebugOn) {
1073         (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1074         fEntriesLinearFitter[fDetectorPreviousTrack]++;
1075       }
1076       fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1077     }
1078   }
1079   if(fExbAltFitOn){
1080     fExbAltFit->Update(fDetectorPreviousTrack,dydx,rms);
1081   }
1082
1083   return kTRUE;
1084 }
1085 //____________Offine tracking in the AliTRDtracker_____________________________
1086 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1087 {
1088   //
1089   // PRF width calibration
1090   // Assume a Gaussian shape: determinate the position of the three pad clusters
1091   // Fit with a straight line
1092   // Take the fitted values for all the clusters (3 or 2 pad clusters)
1093   // Fill the PRF as function of angle of the track
1094   //
1095   //
1096
1097   //printf("begin\n");
1098   ///////////////////////////////////////////
1099   // Take the parameters of the track
1100   //////////////////////////////////////////
1101   // take now the snp, tnp and tgl from the track
1102   Double_t snp = tracklet->GetSnp();             // sin dy/dx at the end of the chamber
1103   Double_t tnp = 0.0;                            // dy/dx at the end of the chamber
1104   if( TMath::Abs(snp) <  1.){
1105     tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1106   }
1107   Double_t tgl  = tracklet->GetTgl();           // dz/dl
1108   Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp);   // dz/dx calculated from dz/dl
1109   // at the entrance
1110   //Double_t tnp = tracklet->GetYref(1);      // dy/dx at the entrance of the chamber
1111   //Double_t tgl = tracklet->GetZref(1);      // dz/dl at the entrance of the chamber
1112   //Double_t dzdx = tgl;                      //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1113   // at the end with correction due to linear fit
1114   //Double_t tnp = tracklet->GetYfit(1);      // dy/dx at the end of the chamber after fit correction
1115   //Double_t tgl = tracklet->GetZfit(1);      // dz/dl at the end of the chamber after fit correction
1116
1117   ///////////////////////////////
1118   // Calculate tnp group shift
1119   ///////////////////////////////
1120   Bool_t echec   = kFALSE;
1121   Double_t shift = 0.0;
1122   //Calculate the shift in x coresponding to this tnp
1123   if(fNgroupprf != 0.0){
1124     shift      = -3.0*(fNgroupprf-1)-1.5;
1125     Double_t limithigh  = -0.2*(fNgroupprf-1);
1126     if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1127     else{
1128       while(tnp > limithigh){
1129         limithigh += 0.2;
1130         shift += 3.0;
1131       }
1132     }
1133   }
1134   // do nothing if out of tnp range
1135   //printf("echec %d\n",(Int_t)echec);
1136   if(echec) return kFALSE;
1137
1138   ///////////////////////
1139   // Variables
1140   //////////////////////
1141
1142   Int_t nb3pc    = 0;              // number of three pads clusters used for fit
1143   // to see the difference between the fit and the 3 pad clusters position
1144   Double_t padPositions[AliTRDseedV1::kNtb];
1145   memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
1146   fLinearFitterTracklet->ClearPoints();
1147
1148   //printf("loop clusters \n");
1149   ////////////////////////////
1150   // loop over the clusters
1151   ////////////////////////////
1152   AliTRDcluster *cl                   = 0x0;
1153   for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1154     // reject shared clusters on pad row
1155     if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
1156       cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1157       if(cl) continue;
1158     }
1159     cl = tracklet->GetClusters(ic);
1160     if(!cl) continue;
1161     Double_t     time  = cl->GetPadTime();
1162     if((time<=7) || (time>=21)) continue;
1163     Short_t  *signals  = cl->GetSignals();
1164     Float_t xcenter    = 0.0;
1165     Bool_t  echec1      = kTRUE;
1166
1167     /////////////////////////////////////////////////////////////
1168     // Center 3 balanced: position with the center of the pad
1169     /////////////////////////////////////////////////////////////
1170     if ((((Float_t) signals[3]) > 0.0) &&
1171         (((Float_t) signals[2]) > 0.0) &&
1172         (((Float_t) signals[4]) > 0.0)) {
1173       echec1 = kFALSE;
1174       // Security if the denomiateur is 0
1175       if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1176            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1177         xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1178           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1179                         / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1180       }
1181       else {
1182         echec1 = kTRUE;
1183       }
1184     }
1185     if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1186     if(echec1) continue;
1187
1188     ////////////////////////////////////////////////////////
1189     //if no echec1: calculate with the position of the pad
1190     // Position of the cluster
1191     // fill the linear fitter
1192     ///////////////////////////////////////////////////////
1193     Double_t       padPosition = xcenter +  cl->GetPadCol();
1194     padPositions[ic]            = padPosition;
1195     nb3pc++;
1196     fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1197
1198
1199   }//clusters loop
1200
1201   //printf("Fin loop clusters \n");
1202   //////////////////////////////
1203   // fit with a straight line
1204   /////////////////////////////
1205   if(nb3pc < 3){
1206     fLinearFitterTracklet->ClearPoints();
1207     return kFALSE;
1208   }
1209   fLinearFitterTracklet->Eval();
1210   TVectorD line(2);
1211   fLinearFitterTracklet->GetParameters(line);
1212   Float_t  pointError  = -1.0;
1213   if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1214   pointError  =  TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1215   }
1216   fLinearFitterTracklet->ClearPoints();
1217
1218   //printf("PRF second loop \n");
1219   ////////////////////////////////////////////////
1220   // Fill the PRF: Second loop over clusters
1221   //////////////////////////////////////////////
1222   for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1223     // reject shared clusters on pad row
1224     cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1225     if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl)) continue;
1226     //
1227     cl = tracklet->GetClusters(ic);
1228     if(!cl) continue;
1229
1230     Short_t  *signals      = cl->GetSignals();              // signal
1231     Double_t     time      = cl->GetPadTime();         // time bin
1232     Float_t padPosTracklet = line[0]+line[1]*time;          // reconstruct position from fit
1233     Float_t padPos         = cl->GetPadCol();               // middle pad
1234     Double_t dpad          = padPosTracklet - padPos;       // reconstruct position relative to middle pad from fit
1235     Float_t ycenter        = 0.0;                           // relative center charge
1236     Float_t ymin           = 0.0;                           // relative left charge
1237     Float_t ymax           = 0.0;                           // relative right charge
1238
1239     ////////////////////////////////////////////////////////////////
1240     // Calculate ycenter, ymin and ymax even for two pad clusters
1241     ////////////////////////////////////////////////////////////////
1242     if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1243        ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1244       Float_t sum     = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1245       if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1246       if(sum > 0.0) ymin    = ((Float_t) signals[2])/ sum;
1247       if(sum > 0.0) ymax    = ((Float_t) signals[4])/ sum;
1248     }
1249
1250     /////////////////////////
1251     // Calibration group
1252     ////////////////////////
1253     Int_t     rowcl        = cl->GetPadRow();                           // row of cluster
1254     Int_t     colcl        = cl->GetPadCol();                           // col of cluster
1255     Int_t     grouplocal   = CalculateCalibrationGroup(2,rowcl,colcl);  // calcul the corresponding group
1256     Int_t     caligroup    = fCalibraMode->GetXbins(2)+ grouplocal;     // calcul the corresponding group
1257     Float_t   xcl          = cl->GetY();                                // y cluster
1258     Float_t   qcl          = cl->GetQ();                                // charge cluster
1259     Int_t     layer        = GetLayer(fDetectorPreviousTrack);          // layer
1260     Int_t     stack        = GetStack(fDetectorPreviousTrack);          // stack
1261     Double_t  xdiff        = dpad;                                      // reconstructed position constant
1262     Double_t  x            = dpad;                                      // reconstructed position moved
1263     Float_t   ep           = pointError;                                // error of fit
1264     Float_t   signal1      = (Float_t)signals[1];                       // signal at the border
1265     Float_t   signal3      = (Float_t)signals[3];                       // signal
1266     Float_t   signal2      = (Float_t)signals[2];                       // signal
1267     Float_t   signal4      = (Float_t)signals[4];                       // signal
1268     Float_t   signal5      = (Float_t)signals[5];                       // signal at the border
1269
1270
1271
1272     /////////////////////
1273     // Debug stuff
1274     ////////////////////
1275
1276     if(fDebugLevel > 1){
1277       if ( !fDebugStreamer ) {
1278         //debug stream
1279         TDirectory *backup = gDirectory;
1280         fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1281         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1282       }
1283
1284       x = xdiff;
1285       Int_t type=0;
1286       Float_t y = ycenter;
1287       (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1288         "caligroup="<<caligroup<<
1289         "detector="<<fDetectorPreviousTrack<<
1290         "layer="<<layer<<
1291         "stack="<<stack<<
1292         "npoints="<<nbclusters<<
1293         "Np="<<nb3pc<<
1294         "ep="<<ep<<
1295         "type="<<type<<
1296         "snp="<<snp<<
1297         "tnp="<<tnp<<
1298         "tgl="<<tgl<<
1299         "dzdx="<<dzdx<<
1300         "padPos="<<padPos<<
1301         "padPosition="<<padPositions[ic]<<
1302         "padPosTracklet="<<padPosTracklet<<
1303         "x="<<x<<
1304         "y="<<y<<
1305         "xcl="<<xcl<<
1306         "qcl="<<qcl<<
1307         "signal1="<<signal1<<
1308         "signal2="<<signal2<<
1309         "signal3="<<signal3<<
1310         "signal4="<<signal4<<
1311         "signal5="<<signal5<<
1312         "time="<<time<<
1313         "\n";
1314       x=-(xdiff+1);
1315       y = ymin;
1316       type=-1;
1317       (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1318         "caligroup="<<caligroup<<
1319         "detector="<<fDetectorPreviousTrack<<
1320         "layer="<<layer<<
1321         "stack="<<stack<<
1322         "npoints="<<nbclusters<<
1323         "Np="<<nb3pc<<
1324         "ep="<<ep<<
1325         "type="<<type<<
1326         "snp="<<snp<<
1327         "tnp="<<tnp<<
1328         "tgl="<<tgl<<
1329         "dzdx="<<dzdx<<
1330         "padPos="<<padPos<<
1331         "padPosition="<<padPositions[ic]<<
1332         "padPosTracklet="<<padPosTracklet<<
1333         "x="<<x<<
1334         "y="<<y<<
1335         "xcl="<<xcl<<
1336         "qcl="<<qcl<<
1337         "signal1="<<signal1<<
1338         "signal2="<<signal2<<
1339         "signal3="<<signal3<<
1340         "signal4="<<signal4<<
1341         "signal5="<<signal5<<
1342         "time="<<time<<
1343         "\n";
1344       x=1-xdiff;
1345       y = ymax;
1346       type=1;
1347       (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1348         "caligroup="<<caligroup<<
1349         "detector="<<fDetectorPreviousTrack<<
1350         "layer="<<layer<<
1351         "stack="<<stack<<
1352         "npoints="<<nbclusters<<
1353         "Np="<<nb3pc<<
1354         "ep="<<ep<<
1355         "type="<<type<<
1356         "snp="<<snp<<
1357         "tnp="<<tnp<<
1358         "tgl="<<tgl<<
1359         "dzdx="<<dzdx<<
1360         "padPos="<<padPos<<
1361         "padPosition="<<padPositions[ic]<<
1362         "padPosTracklet="<<padPosTracklet<<
1363         "x="<<x<<
1364         "y="<<y<<
1365         "xcl="<<xcl<<
1366         "qcl="<<qcl<<
1367         "signal1="<<signal1<<
1368         "signal2="<<signal2<<
1369         "signal3="<<signal3<<
1370         "signal4="<<signal4<<
1371         "signal5="<<signal5<<
1372         "time="<<time<<
1373         "\n";
1374
1375     }
1376
1377     /////////////////////
1378     // Cuts quality
1379     /////////////////////
1380     if(nbclusters < fNumberClusters) continue;
1381     if(nbclusters > fNumberClustersf) continue;
1382     if(nb3pc <= 5) continue;
1383     if((time >= 21) || (time < 7)) continue;
1384     if(TMath::Abs(qcl) < 80) continue;
1385     if( TMath::Abs(snp) >  1.) continue;
1386
1387
1388     ////////////////////////
1389     // Fill the histos
1390     ///////////////////////
1391     if (fHisto2d) {
1392       if(TMath::Abs(dpad) < 1.5) {
1393         fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1394         fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1395         //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1396       }
1397       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1398         fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1399         fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1400       }
1401       if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1402         fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1403         fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1404       }
1405     }
1406     // vector method
1407     if (fVector2d) {
1408       if(TMath::Abs(dpad) < 1.5) {
1409         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1410         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1411       }
1412       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1413         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1414         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1415       }
1416       if((ymax > 0.0)  && (TMath::Abs(dpad-1.0) < 1.5)) {
1417         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1418         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1419       }
1420     }
1421   } // second loop over clusters
1422
1423
1424   return kTRUE;
1425 }
1426 ///////////////////////////////////////////////////////////////////////////////////////
1427 // Pad row col stuff: see if masked or not
1428 ///////////////////////////////////////////////////////////////////////////////////////
1429 //_____________________________________________________________________________
1430 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(const AliTRDcluster *cl)
1431 {
1432   //
1433   // See if we are not near a masked pad
1434   //
1435
1436   if(cl->IsMasked()) fGoodTracklet = kFALSE;
1437
1438
1439 }
1440 //_____________________________________________________________________________
1441 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(const Int_t detector,const Int_t row,const Int_t col)
1442 {
1443   //
1444   // See if we are not near a masked pad
1445   //
1446
1447   if (!IsPadOn(detector, col, row)) {
1448     fGoodTracklet = kFALSE;
1449   }
1450
1451   if (col > 0) {
1452     if (!IsPadOn(detector, col-1, row)) {
1453       fGoodTracklet = kFALSE;
1454     }
1455   }
1456
1457   if (col < 143) {
1458     if (!IsPadOn(detector, col+1, row)) {
1459       fGoodTracklet = kFALSE;
1460     }
1461   }
1462
1463 }
1464 //_____________________________________________________________________________
1465 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1466 {
1467   //
1468   // Look in the choosen database if the pad is On.
1469   // If no the track will be "not good"
1470   //
1471
1472   // Get the parameter object
1473   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1474   if (!cal) {
1475     AliInfo("Could not get calibDB");
1476     return kFALSE;
1477   }
1478
1479   if (!cal->IsChamberGood(detector)     ||
1480        cal->IsChamberNoData(detector)        ||
1481        cal->IsPadMasked(detector,col,row)) {
1482     return kFALSE;
1483   }
1484   else {
1485     return kTRUE;
1486   }
1487
1488 }
1489 ///////////////////////////////////////////////////////////////////////////////////////
1490 // Calibration groups: calculate the number of groups, localise...
1491 ////////////////////////////////////////////////////////////////////////////////////////
1492 //_____________________________________________________________________________
1493 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1494 {
1495   //
1496   // Calculate the calibration group number for i
1497   //
1498
1499   // Row of the cluster and position in the pad groups
1500   Int_t posr = 0;
1501   if (fCalibraMode->GetNnZ(i) != 0) {
1502     posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1503   }
1504
1505
1506   // Col of the cluster and position in the pad groups
1507   Int_t posc = 0;
1508   if (fCalibraMode->GetNnRphi(i) != 0) {
1509     posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1510   }
1511
1512   return posc*fCalibraMode->GetNfragZ(i)+posr;
1513
1514 }
1515 //____________________________________________________________________________________
1516 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1517 {
1518   //
1519   // Calculate the total number of calibration groups
1520   //
1521
1522   Int_t ntotal = 0;
1523
1524   // All together
1525   if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1526     ntotal = 1;
1527     AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1528     return ntotal;
1529   }
1530
1531   // Per Supermodule
1532   if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1533     ntotal = 18;
1534     AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1535     return ntotal;
1536   }
1537
1538   // More
1539   fCalibraMode->ModePadCalibration(2,i);
1540   fCalibraMode->ModePadFragmentation(0,2,0,i);
1541   fCalibraMode->SetDetChamb2(i);
1542   ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1543   fCalibraMode->ModePadCalibration(0,i);
1544   fCalibraMode->ModePadFragmentation(0,0,0,i);
1545   fCalibraMode->SetDetChamb0(i);
1546   ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1547   AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1548   return ntotal;
1549
1550 }
1551 //_____________________________________________________________________________
1552 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1553 {
1554   //
1555   // Set the mode of calibration group in the z direction for the parameter i
1556   //
1557
1558   if ((Nz >= 0) &&
1559       (Nz <  5)) {
1560     fCalibraMode->SetNz(i, Nz);
1561   }
1562   else {
1563     AliInfo("You have to choose between 0 and 4");
1564   }
1565
1566 }
1567 //_____________________________________________________________________________
1568 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1569 {
1570   //
1571   // Set the mode of calibration group in the rphi direction for the parameter i
1572   //
1573
1574   if ((Nrphi >= 0) &&
1575       (Nrphi <  7)) {
1576     fCalibraMode->SetNrphi(i ,Nrphi);
1577   }
1578   else {
1579     AliInfo("You have to choose between 0 and 6");
1580   }
1581
1582 }
1583
1584 //_____________________________________________________________________________
1585 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1586 {
1587   //
1588   // Set the mode of calibration group all together
1589   //
1590   if(fVector2d == kTRUE) {
1591     AliInfo("Can not work with the vector method");
1592     return;
1593   }
1594   fCalibraMode->SetAllTogether(i);
1595
1596 }
1597
1598 //_____________________________________________________________________________
1599 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1600 {
1601   //
1602   // Set the mode of calibration group per supermodule
1603   //
1604   if(fVector2d == kTRUE) {
1605     AliInfo("Can not work with the vector method");
1606     return;
1607   }
1608   fCalibraMode->SetPerSuperModule(i);
1609
1610 }
1611
1612 //____________Set the pad calibration variables for the detector_______________
1613 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1614 {
1615   //
1616   // For the detector calcul the first Xbins and set the number of row
1617   // and col pads per calibration groups, the number of calibration
1618   // groups in the detector.
1619   //
1620
1621   // first Xbins of the detector
1622   if (fCH2dOn) {
1623     fCalibraMode->CalculXBins(detector,0);
1624   }
1625   if (fPH2dOn) {
1626     fCalibraMode->CalculXBins(detector,1);
1627   }
1628   if (fPRF2dOn) {
1629     fCalibraMode->CalculXBins(detector,2);
1630   }
1631
1632   // fragmentation of idect
1633   for (Int_t i = 0; i < 3; i++) {
1634     fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1635     fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1636                        , (Int_t) GetStack(detector)
1637                        , (Int_t) GetSector(detector),i);
1638   }
1639
1640   return kTRUE;
1641
1642 }
1643 //_____________________________________________________________________________
1644 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1645 {
1646   //
1647   // Should be between 0 and 6
1648   //
1649
1650   if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1651     AliInfo("The number of groups must be between 0 and 6!");
1652   }
1653   else {
1654     fNgroupprf = numberGroupsPRF;
1655   }
1656
1657 }
1658 ///////////////////////////////////////////////////////////////////////////////////////////
1659 // Per tracklet: store or reset the info, fill the histos with the info
1660 //////////////////////////////////////////////////////////////////////////////////////////
1661 //_____________________________________________________________________________
1662 Float_t AliTRDCalibraFillHisto::StoreInfoCHPHtrack(const AliTRDcluster *cl,const Double_t dqdl,const Int_t *group,const Int_t row,const Int_t col,const AliTRDcluster *cls)
1663 {
1664   //
1665   // Store the infos in fAmpTotal, fPHPlace and fPHValue
1666   // Correct from the gain correction before
1667   // cls is shared cluster if any
1668   // Return the charge
1669   //
1670   //
1671
1672   //printf("StoreInfoCHPHtrack\n");
1673
1674   // time bin of the cluster not corrected
1675   Int_t    time     = cl->GetPadTime();
1676   Float_t  charge   = TMath::Abs(cl->GetQ());
1677   if(cls) {
1678     charge += TMath::Abs(cls->GetQ());
1679     //printf("AliTRDCalibraFillHisto::Add the cluster charge");
1680   }
1681
1682   //printf("Store::time %d, amplitude %f\n",time,dqdl);
1683
1684   //Correct for the gain coefficient used in the database for reconstruction
1685   Float_t correctthegain = 1.0;
1686   if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1687   else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1688   Float_t correction    = 1.0;
1689   Float_t normalisation = 1.13; //org: 6.67; 1st: 1.056; 2nd: 1.13;
1690   // we divide with gain in AliTRDclusterizer::Transform...
1691   if( correctthegain > 0 ) normalisation /= correctthegain;
1692
1693
1694   // take dd/dl corrected from the angle of the track
1695   correction = dqdl / (normalisation);
1696
1697
1698   // Fill the fAmpTotal with the charge
1699   if (fCH2dOn) {
1700     if((!fLimitChargeIntegration) || (cl->IsInChamber())) {
1701       //printf("Store::group %d, amplitude %f\n",group[0],correction);
1702       fAmpTotal[(Int_t) group[0]] += correction;
1703     }
1704   }
1705
1706   // Fill the fPHPlace and value
1707   if (fPH2dOn) {
1708     if (time>=0 && time<fTimeMax) {
1709       fPHPlace[time] = group[1];
1710       fPHValue[time] = charge;
1711     }
1712   }
1713
1714   return correction;
1715
1716 }
1717 //____________Offine tracking in the AliTRDtracker_____________________________
1718 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1719 {
1720   //
1721   // Reset values per tracklet
1722   //
1723
1724   //Reset good tracklet
1725   fGoodTracklet = kTRUE;
1726
1727   // Reset the fPHValue
1728   if (fPH2dOn) {
1729     //Reset the fPHValue and fPHPlace
1730     for (Int_t k = 0; k < fTimeMax; k++) {
1731       fPHValue[k] = 0.0;
1732       fPHPlace[k] = -1;
1733     }
1734   }
1735
1736   // Reset the fAmpTotal where we put value
1737   if (fCH2dOn) {
1738     for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1739       fAmpTotal[k] = 0.0;
1740     }
1741   }
1742 }
1743 //____________Offine tracking in the AliTRDtracker_____________________________
1744 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1745 {
1746   //
1747   // For the offline tracking or mcm tracklets
1748   // This function will be called in the functions UpdateHistogram...
1749   // to fill the info of a track for the relativ gain calibration
1750   //
1751
1752   Int_t nb            =  0;   // Nombre de zones traversees
1753   Int_t fd            = -1;   // Premiere zone non nulle
1754   Float_t totalcharge = 0.0;  // Total charge for the supermodule histo
1755
1756   //printf("CH2d nbclusters %d, fNumberClusters %d, fNumberClustersf %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1757
1758   if(nbclusters < fNumberClusters) return;
1759   if(nbclusters > fNumberClustersf) return;
1760
1761
1762   // Normalize with the number of clusters
1763   Double_t normalizeCst = fRelativeScale;
1764   if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1765
1766   //printf("Number of groups in one detector %d\n",fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0));
1767
1768   // See if the track goes through different zones
1769   for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1770     //printf("fAmpTotal %f for %d\n",fAmpTotal[k],k);
1771     if (fAmpTotal[k] > 0.0) {
1772       totalcharge += fAmpTotal[k];
1773       nb++;
1774       if (nb == 1) {
1775         fd = k;
1776       }
1777     }
1778   }
1779
1780   //printf("CH2d: nb %d, fd %d, calibration group %d, amplitude %f, detector %d\n",nb,fd,fCalibraMode->GetXbins(0),fAmpTotal[fd]/normalizeCst,fDetectorPreviousTrack);
1781
1782   switch (nb)
1783     {
1784     case 1:
1785       fNumberUsedCh[0]++;
1786       fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1787       if (fHisto2d) {
1788         FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1789         //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1790       }
1791       if (fVector2d) {
1792         fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1793       }
1794       break;
1795     case 2:
1796       if ((fAmpTotal[fd]   > 0.0) &&
1797           (fAmpTotal[fd+1] > 0.0)) {
1798         // One of the two very big
1799         if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1800           if (fHisto2d) {
1801             FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1802             //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1803           }
1804           if (fVector2d) {
1805             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1806           }
1807           fNumberUsedCh[1]++;
1808           fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1809         }
1810         if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1811           if (fHisto2d) {
1812             FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
1813             //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
1814           }
1815           if (fVector2d) {
1816             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
1817           }
1818           fNumberUsedCh[1]++;
1819           fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1820         }
1821       }
1822       if (fCalibraMode->GetNfragZ(0) > 1) {
1823         if (fAmpTotal[fd] > 0.0) {
1824           if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1825             if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1826               // One of the two very big
1827               if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1828                 if (fHisto2d) {
1829                   FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1830                   //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1831                 }
1832                 if (fVector2d) {
1833                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1834                 }
1835                 fNumberUsedCh[1]++;
1836                 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1837               }
1838               if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1839                 if (fHisto2d) {
1840                   FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1841                   //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1842                 }
1843                 fNumberUsedCh[1]++;
1844                 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1845                 if (fVector2d) {
1846                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1847                 }
1848               }
1849             }
1850           }
1851         }
1852       }
1853       break;
1854     default: break;
1855   }
1856 }
1857 //____________Offine tracking in the AliTRDtracker_____________________________
1858 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1859 {
1860   //
1861   // For the offline tracking or mcm tracklets
1862   // This function will be called in the functions UpdateHistogram...
1863   // to fill the info of a track for the drift velocity  calibration
1864   //
1865
1866   Int_t nb  =  1; // Nombre de zones traversees 1, 2 ou plus de 3
1867   Int_t fd1 = -1; // Premiere zone non nulle
1868   Int_t fd2 = -1; // Deuxieme zone non nulle
1869   Int_t k1  = -1; // Debut de la premiere zone
1870   Int_t k2  = -1; // Debut de la seconde zone
1871   Int_t nbclusters = 0; // number of clusters
1872
1873
1874
1875   // See if the track goes through different zones
1876   for (Int_t k = 0; k < fTimeMax; k++) {
1877     if (fPHValue[k] > 0.0) {
1878       nbclusters++;
1879       if (fd1 == -1) {
1880         fd1 = fPHPlace[k];
1881         k1  = k;
1882       }
1883       if (fPHPlace[k] != fd1) {
1884         if (fd2 == -1) {
1885           k2  = k;
1886           fd2 = fPHPlace[k];
1887           nb  = 2;
1888         }
1889         if (fPHPlace[k] != fd2) {
1890           nb = 3;
1891         }
1892       }
1893     }
1894   }
1895
1896   // See if noise before and after
1897   if(fMaxCluster > 0) {
1898     if(fPHValue[0] > fMaxCluster) return;
1899     if(fTimeMax > fNbMaxCluster) {
1900       for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
1901         if(fPHValue[k] > fMaxCluster) return;
1902       }
1903     }
1904   }
1905
1906   //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1907
1908   if(nbclusters < fNumberClusters) return;
1909   if(nbclusters > fNumberClustersf) return;
1910
1911   switch(nb)
1912     {
1913     case 1:
1914       fNumberUsedPh[0]++;
1915       for (Int_t i = 0; i < fTimeMax; i++) {
1916         if (fHisto2d) {
1917           if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1918           else {
1919             if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1920               }
1921           //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
1922         }
1923         if (fVector2d) {
1924           if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1925           else {
1926             if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1927           }
1928         }
1929       }
1930       break;
1931     case 2:
1932       if ((fd1 == fd2+1) ||
1933           (fd2 == fd1+1)) {
1934         // One of the two fast all the think
1935         if (k2 > (k1+fDifference)) {
1936           //we choose to fill the fd1 with all the values
1937           fNumberUsedPh[1]++;
1938           for (Int_t i = 0; i < fTimeMax; i++) {
1939             if (fHisto2d) {
1940               if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1941               else {
1942                 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1943                   }
1944             }
1945             if (fVector2d) {
1946               if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1947               else {
1948                 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1949                   }
1950             }
1951           }
1952         }
1953         if ((k2+fDifference) < fTimeMax) {
1954           //we choose to fill the fd2 with all the values
1955           fNumberUsedPh[1]++;
1956           for (Int_t i = 0; i < fTimeMax; i++) {
1957             if (fHisto2d) {
1958               if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1959               else {
1960                 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1961               }
1962             }
1963           if (fVector2d) {
1964             if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1965             else {
1966               if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1967             }
1968           }
1969           }
1970         }
1971       }
1972       // Two zones voisines sinon rien!
1973       if (fCalibraMode->GetNfragZ(1) > 1) {
1974         // Case 2
1975         if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1976           if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1977             // One of the two fast all the think
1978             if (k2 > (k1+fDifference)) {
1979               //we choose to fill the fd1 with all the values
1980               fNumberUsedPh[1]++;
1981               for (Int_t i = 0; i < fTimeMax; i++) {
1982                 if (fHisto2d) {
1983                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1984                   else {
1985                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1986                   }
1987                 }
1988                 if (fVector2d) {
1989                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1990                   else {
1991                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1992                   }
1993                 }
1994               }
1995             }
1996             if ((k2+fDifference) < fTimeMax) {
1997               //we choose to fill the fd2 with all the values
1998               fNumberUsedPh[1]++;
1999               for (Int_t i = 0; i < fTimeMax; i++) {
2000                 if (fHisto2d) {
2001                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2002                   else {
2003                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2004                   }
2005                 }
2006                 if (fVector2d) {
2007                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2008                   else {
2009                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2010                   }
2011                 }
2012               }
2013             }
2014           }
2015         }
2016         // Two zones voisines sinon rien!
2017         // Case 3
2018         if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2019           if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2020             // One of the two fast all the think
2021             if (k2 > (k1 + fDifference)) {
2022               //we choose to fill the fd1 with all the values
2023               fNumberUsedPh[1]++;
2024               for (Int_t i = 0; i < fTimeMax; i++) {
2025                 if (fHisto2d) {
2026                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2027                   else {
2028                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2029                   }
2030                 }
2031                 if (fVector2d) {
2032                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2033                   else {
2034                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2035                   }
2036                 }
2037               }
2038             }
2039             if ((k2+fDifference) < fTimeMax) {
2040               //we choose to fill the fd2 with all the values
2041               fNumberUsedPh[1]++;
2042               for (Int_t i = 0; i < fTimeMax; i++) {
2043                 if (fHisto2d) {
2044                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2045                   else {
2046                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2047                   }
2048                 }
2049                 if (fVector2d) {
2050                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2051                   else {
2052                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2053                   }
2054                 }
2055               }
2056             }
2057           }
2058         }
2059       }
2060       break;
2061     default: break;
2062     }
2063 }
2064 //////////////////////////////////////////////////////////////////////////////////////////
2065 // DAQ process functions
2066 /////////////////////////////////////////////////////////////////////////////////////////
2067 //_____________________________________________________________________
2068 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader)
2069  { //main
2070   //
2071   // Event Processing loop - AliTRDrawStream
2072   //
2073   // 0 timebin problem
2074   // 1 no input
2075   // 2 input
2076   // Same algorithm as TestBeam but different reader
2077   //
2078
2079   AliTRDrawStream *rawStream = new AliTRDrawStream(rawReader);
2080
2081   AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
2082   digitsManager->CreateArrays();
2083
2084   Int_t withInput = 1;
2085
2086   Double_t phvalue[16][144][36];
2087   for(Int_t k = 0; k < 36; k++){
2088     for(Int_t j = 0; j < 16; j++){
2089       for(Int_t c = 0; c < 144; c++){
2090         phvalue[j][c][k] = 0.0;
2091       }
2092     }
2093   }
2094
2095   fDetectorPreviousTrack = -1;
2096   fMCMPrevious           = -1;
2097   fROBPrevious           = -1;
2098
2099   Int_t nbtimebin = 0;
2100   Int_t baseline  = 10;
2101
2102
2103     fTimeMax = 0;
2104
2105     Int_t det    = 0;
2106     while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
2107
2108       if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
2109         //      printf("there is ADC data on this chamber!\n");
2110
2111         AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
2112         if (digits->HasData()) { //array
2113
2114           AliTRDSignalIndex   *indexes = digitsManager->GetIndexes(det);
2115           if (indexes->IsAllocated() == kFALSE) {
2116             AliError("Indexes do not exist!");
2117             break;
2118           }
2119           Int_t iRow  = 0;
2120           Int_t iCol  = 0;
2121           indexes->ResetCounters();
2122
2123           while (indexes->NextRCIndex(iRow, iCol)) { //column,row
2124             //printf(" det %d \t row %d \t col %d \t digit\n",det,iRow,iCol);
2125             //while (rawStream->Next()) {
2126
2127             Int_t idetector = det;                                             //  current detector
2128             //Int_t imcm      = rawStream->GetMCM();                            //  current MCM
2129             //Int_t irob      = rawStream->GetROB();                            //  current ROB
2130
2131
2132             if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)) {
2133               // Fill
2134               withInput = TMath::Max(FillDAQ(phvalue),withInput);
2135
2136               // reset
2137               for(Int_t k = 0; k < 36; k++){
2138                 for(Int_t j = 0; j < 16; j++){
2139                   for(Int_t c = 0; c < 144; c++){
2140                     phvalue[j][c][k] = 0.0;
2141                   }
2142                 }
2143               }
2144             }
2145
2146             fDetectorPreviousTrack = idetector;
2147             //fMCMPrevious           = imcm;
2148             //fROBPrevious           = irob;
2149
2150             //    nbtimebin              = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
2151             AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
2152             nbtimebin              = digitParam->GetNTimeBins(det);              //  number of time bins read from data
2153             baseline               = digitParam->GetADCbaseline(det);            //  baseline
2154
2155             if(nbtimebin == 0) return 0;
2156             if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2157             fTimeMax          = nbtimebin;
2158
2159             fNumberClustersf    = fTimeMax;
2160             fNumberClusters     = (Int_t)(fNumberClustersProcent*fTimeMax);
2161
2162
2163             for(Int_t itime = 0; itime < nbtimebin; itime++) {
2164               //            phvalue[row][col][itime] = signal[itime]-baseline;
2165               phvalue[iRow][iCol][itime] = (Short_t)(digits->GetData(iRow,iCol,itime) - baseline);
2166               /*if(phvalue[iRow][iCol][itime] >= 20) {
2167                  printf("----------> phvalue[%d][%d][%d] %d  baseline %d \n",
2168                        iRow,
2169                        iCol,
2170                        itime,
2171                        (Short_t)(digits->GetData(iRow,iCol,itime)),
2172                        baseline);
2173                        }*/
2174             }
2175
2176           }//column,row
2177
2178           // fill the last one
2179           if(fDetectorPreviousTrack != -1){
2180
2181             // Fill
2182             withInput = TMath::Max(FillDAQ(phvalue),withInput);
2183             //      printf("\n ---> withinput %d\n\n",withInput);
2184             // reset
2185             for(Int_t k = 0; k < 36; k++){
2186               for(Int_t j = 0; j < 16; j++){
2187                 for(Int_t c = 0; c < 144; c++){
2188                   phvalue[j][c][k] = 0.0;
2189                 }
2190               }
2191             }
2192           }
2193
2194         }//array
2195       }//QA
2196       digitsManager->ClearArrays(det);
2197     }//idetector
2198     delete digitsManager;
2199
2200     delete rawStream;
2201     return withInput;
2202  }//main
2203 //_____________________________________________________________________
2204 //////////////////////////////////////////////////////////////////////////////
2205 // Routine inside the DAQ process
2206 /////////////////////////////////////////////////////////////////////////////
2207 //_______________________________________________________________________
2208 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2209
2210   //
2211   // Look for the maximum by collapsing over the time
2212   // Sum over four pad col and two pad row
2213   //
2214
2215   Int_t used = 0;
2216
2217
2218   Int_t idect = fDetectorPreviousTrack;
2219   //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2220   Double_t sum[36];
2221   for(Int_t tb = 0; tb < 36; tb++){
2222     sum[tb] = 0.0;
2223   }
2224
2225   //fGoodTracklet = kTRUE;
2226   //fDetectorPreviousTrack = 0;
2227
2228
2229   ///////////////////////////
2230   // look for maximum
2231   /////////////////////////
2232
2233   Int_t imaxRow = 0;
2234   Int_t imaxCol = 0;
2235   Double_t integralMax = -1;
2236
2237   for (Int_t ir = 1; ir <= 15; ir++)
2238     {
2239       for (Int_t ic = 2; ic <= 142; ic++)
2240         {
2241           Double_t integral = 0;
2242           for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2243             {
2244               for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2245                 {
2246                   if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2247                       ic + ishiftC >= 1 && ic + ishiftC <= 144)
2248                     {
2249
2250                       for(Int_t tb = 0; tb< fTimeMax; tb++){
2251                         integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2252                       }// addtb
2253                     } //addsignal
2254                 } //shiftC
2255             } // shiftR
2256           if (integralMax < integral)
2257             {
2258               imaxRow = ir;
2259               imaxCol = ic;
2260               integralMax = integral;
2261
2262             } // check max integral
2263         } //ic
2264     } // ir
2265
2266   //  printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2267   //if((imaxRow == 0) || (imaxRow >= 15) || (imaxCol <= 3) || (imaxCol >= 140)) {
2268   //  used=1;
2269   //  return used;
2270   // }
2271
2272   if(((imaxRow + fNumberRowDAQ) > 16) || (imaxRow == 0) || ((imaxCol - fNumberColDAQ) <= 1) || ((imaxCol + fNumberColDAQ) >= 144)) {
2273     used=1;
2274     return used;
2275   }
2276   //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2277   //if(!fGoodTracklet) used = 1;;
2278
2279   //  /////////////////////////////////////////////////////
2280   // sum ober 2 row and 4 pad cols for each time bins
2281   //  ////////////////////////////////////////////////////
2282
2283
2284
2285   for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2286     {
2287       for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2288         {
2289           if (imaxRow + ishiftR >= 1 && imaxRow + ishiftR <= 16 &&
2290               imaxCol + ishiftC >= 1 && imaxCol + ishiftC <= 144)
2291             {
2292               for(Int_t it = 0; it < fTimeMax; it++){
2293                 sum[it] += phvalue[imaxRow + ishiftR-1][imaxCol + ishiftC-1][it];
2294               }
2295             }
2296         } // col shift
2297     }// row shift
2298
2299   Int_t nbcl = 0;
2300   Double_t sumcharge = 0.0;
2301   for(Int_t it = 0; it < fTimeMax; it++){
2302     sumcharge += sum[it];
2303     if(sum[it] > fThresholdClustersDAQ)  nbcl++;
2304   }
2305
2306
2307   /////////////////////////////////////////////////////////
2308   // Debug
2309   ////////////////////////////////////////////////////////
2310   if(fDebugLevel > 1){
2311     if ( !fDebugStreamer ) {
2312       //debug stream
2313       TDirectory *backup = gDirectory;
2314       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2315       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2316     }
2317
2318     Double_t amph0 = sum[0];
2319     Double_t amphlast = sum[fTimeMax-1];
2320     Double_t rms      = TMath::RMS(fTimeMax,sum);
2321     Int_t    goodtracklet = (Int_t) fGoodTracklet;
2322     for(Int_t it = 0; it < fTimeMax; it++){
2323       Double_t clustera = sum[it];
2324
2325     (* fDebugStreamer) << "FillDAQa"<<
2326       "ampTotal="<<sumcharge<<
2327       "row="<<imaxRow<<
2328       "col="<<imaxCol<<
2329       "detector="<<idect<<
2330       "amph0="<<amph0<<
2331       "amphlast="<<amphlast<<
2332       "goodtracklet="<<goodtracklet<<
2333       "clustera="<<clustera<<
2334       "it="<<it<<
2335       "rms="<<rms<<
2336       "nbcl="<<nbcl<<
2337       "\n";
2338     }
2339   }
2340
2341   ////////////////////////////////////////////////////////
2342   // fill
2343   ///////////////////////////////////////////////////////
2344   //printf("fNumberClusters %d, fNumberClustersf %d\n",fNumberClusters,fNumberClustersf);
2345   if(sum[0] > 100.0) used = 1;
2346   if(nbcl < fNumberClusters) used = 1;
2347   if(nbcl > fNumberClustersf) used = 1;
2348
2349   //if(fDetectorPreviousTrack == 15){
2350   //  printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2351   //}
2352   //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2353   if(used == 0){
2354     for(Int_t it = 0; it < fTimeMax; it++){
2355       if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2356       else{
2357         if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2358       }
2359       //if(fFillWithZero) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2360       //else{
2361       // if(sum[it] > 0.0) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2362       //}
2363     }
2364
2365
2366     //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2367     used = 2;
2368     //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2369
2370   }
2371
2372   return used;
2373
2374 }
2375 //____________Online trackling in AliTRDtrigger________________________________
2376 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2377 {
2378   //
2379   // For the DAQ
2380   // Fill a simple average pulse height
2381   //
2382
2383
2384   ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2385
2386
2387   return kTRUE;
2388
2389 }
2390 //____________Write_____________________________________________________
2391 //_____________________________________________________________________
2392 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2393 {
2394   //
2395   //  Write infos to file
2396   //
2397
2398   //For debugging
2399   if ( fDebugStreamer ) {
2400     delete fDebugStreamer;
2401     fDebugStreamer = 0x0;
2402   }
2403
2404   AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2405                ,fNumberTrack
2406                ,fNumberUsedCh[0]
2407                ,fNumberUsedCh[1]
2408                ,fNumberUsedPh[0]
2409                ,fNumberUsedPh[1]));
2410
2411   TDirectory *backup = gDirectory;
2412   TString option;
2413
2414   if ( append )
2415     option = "update";
2416   else
2417     option = "recreate";
2418
2419   TFile f(filename,option.Data());
2420
2421   TStopwatch stopwatch;
2422   stopwatch.Start();
2423   if(fVector2d) {
2424     f.WriteTObject(fCalibraVector);
2425   }
2426
2427   if (fCH2dOn ) {
2428     if (fHisto2d) {
2429       f.WriteTObject(fCH2d);
2430     }
2431   }
2432   if (fPH2dOn ) {
2433     if (fHisto2d) {
2434       f.WriteTObject(fPH2d);
2435     }
2436   }
2437   if (fPRF2dOn) {
2438     if (fHisto2d) {
2439         f.WriteTObject(fPRF2d);
2440     }
2441   }
2442   if(fLinearFitterOn){
2443     if(fLinearFitterDebugOn) AnalyseLinearFitter();
2444     f.WriteTObject(fLinearVdriftFit);
2445   }
2446
2447   f.Close();
2448
2449   if ( backup ) backup->cd();
2450
2451   AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2452                ,stopwatch.RealTime(),stopwatch.CpuTime()));
2453 }
2454 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2455 // Stats stuff
2456 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2457 //___________________________________________probe the histos__________________________________________________
2458 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2459 {
2460   //
2461   // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2462   // debug mode with 2 for TH2I and 3 for TProfile2D
2463   // It gives a pointer to a Double_t[7] with the info following...
2464   // [0] : number of calibration groups with entries
2465   // [1] : minimal number of entries found
2466   // [2] : calibration group number of the min
2467   // [3] : maximal number of entries found
2468   // [4] : calibration group number of the max
2469   // [5] : mean number of entries found
2470   // [6] : mean relative error
2471   //
2472
2473   Double_t *info = new Double_t[7];
2474
2475   // Number of Xbins (detectors or groups of pads)
2476   Int_t    nbins   = h->GetNbinsY(); //number of calibration groups
2477   Int_t    nxbins  = h->GetNbinsX(); //number of bins per histo
2478
2479   // Initialise
2480   Double_t nbwe = 0; //number of calibration groups with entries
2481   Double_t minentries = 0; //minimal number of entries found
2482   Double_t maxentries = 0; //maximal number of entries found
2483   Double_t placemin = 0; //calibration group number of the min
2484   Double_t placemax = -1; //calibration group number of the max
2485   Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2486   Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2487
2488   Double_t counter = 0;
2489
2490   //Debug
2491   TH1F *nbEntries = 0x0;//distribution of the number of entries
2492   TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2493   TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2494
2495   // Beginning of the loop over the calibration groups
2496   for (Int_t idect = 0; idect < nbins; idect++) {
2497
2498     TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2499     projch->SetDirectory(0);
2500
2501     // Number of entries for this calibration group
2502     Double_t nentries = 0.0;
2503     if((i%2) == 0){
2504       for (Int_t k = 0; k < nxbins; k++) {
2505         nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2506       }
2507     }
2508     else{
2509       for (Int_t k = 0; k < nxbins; k++) {
2510         nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2511         if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2512           meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
2513           counter++;
2514         }
2515       }
2516     }
2517
2518     //Debug
2519     if(i > 1){
2520       if(nentries > 0){
2521         if(!((Bool_t)nbEntries)) nbEntries = new TH1F("Number of entries","Number of entries",100,(Int_t)nentries/2,nentries*2);
2522         nbEntries->SetDirectory(0);
2523         nbEntries->Fill(nentries);
2524         if(!((Bool_t)nbEntriesPerGroup)) nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group",nbins,0,nbins);
2525         nbEntriesPerGroup->SetDirectory(0);
2526         nbEntriesPerGroup->Fill(idect+0.5,nentries);
2527         if(!((Bool_t)nbEntriesPerSp)) nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule",(Int_t)(nbins/18),0,(Int_t)(nbins/18));
2528         nbEntriesPerSp->SetDirectory(0);
2529         nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
2530       }
2531     }
2532
2533     //min amd max
2534     if(nentries > maxentries){
2535       maxentries = nentries;
2536       placemax = idect;
2537     }
2538     if(idect == 0) {
2539       minentries = nentries;
2540     }
2541     if(nentries < minentries){
2542       minentries = nentries;
2543       placemin = idect;
2544     }
2545     //nbwe
2546     if(nentries > 0) {
2547       nbwe++;
2548       meanstats += nentries;
2549     }
2550   }//calibration groups loop
2551
2552   if(nbwe > 0) meanstats /= nbwe;
2553   if(counter > 0) meanrelativerror /= counter;
2554
2555   AliInfo(Form("There are %f calibration groups with entries",nbwe));
2556   AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
2557   AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
2558   AliInfo(Form("The mean number of entries is %f",meanstats));
2559   if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
2560
2561   info[0] = nbwe;
2562   info[1] = minentries;
2563   info[2] = placemin;
2564   info[3] = maxentries;
2565   info[4] = placemax;
2566   info[5] = meanstats;
2567   info[6] = meanrelativerror;
2568
2569   if(nbEntries && nbEntriesPerSp && nbEntriesPerGroup){
2570     gStyle->SetPalette(1);
2571     gStyle->SetOptStat(1111);
2572     gStyle->SetPadBorderMode(0);
2573     gStyle->SetCanvasColor(10);
2574     gStyle->SetPadLeftMargin(0.13);
2575     gStyle->SetPadRightMargin(0.01);
2576     TCanvas *stat = new TCanvas("stat","",50,50,600,800);
2577     stat->Divide(2,1);
2578     stat->cd(1);
2579     nbEntries->Draw("");
2580     stat->cd(2);
2581     nbEntriesPerSp->SetStats(0);
2582     nbEntriesPerSp->Draw("");
2583     TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
2584     stat1->cd();
2585     nbEntriesPerGroup->SetStats(0);
2586     nbEntriesPerGroup->Draw("");
2587   }
2588
2589   return info;
2590
2591 }
2592 //____________________________________________________________________________
2593 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
2594 {
2595   //
2596   // Return a Int_t[4] with:
2597   // 0 Mean number of entries
2598   // 1 median of number of entries
2599   // 2 rms of number of entries
2600   // 3 number of group with entries
2601   //
2602
2603   Double_t *stat = new Double_t[4];
2604   stat[3]             = 0.0;
2605
2606   Int_t    nbofgroups = CalculateTotalNumberOfBins(0);
2607
2608   Double_t *weight = new Double_t[nbofgroups];
2609   Double_t *nonul = new Double_t[nbofgroups];
2610
2611   for(Int_t k = 0; k < nbofgroups; k++){
2612     if(fEntriesCH[k] > 0) {
2613       weight[k] = 1.0;
2614       nonul[(Int_t)stat[3]] = fEntriesCH[k];
2615       stat[3]++;
2616     }
2617     else weight[k] = 0.0;
2618   }
2619   stat[0]          = TMath::Mean(nbofgroups,fEntriesCH,weight);
2620   stat[1]          = TMath::Median(nbofgroups,fEntriesCH,weight);
2621   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
2622
2623   delete [] weight;
2624   delete [] nonul;
2625
2626   return stat;
2627
2628 }
2629 //____________________________________________________________________________
2630 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
2631 {
2632   //
2633   // Return a Int_t[4] with:
2634   // 0 Mean number of entries
2635   // 1 median of number of entries
2636   // 2 rms of number of entries
2637   // 3 number of group with entries
2638   //
2639
2640   Double_t *stat      = new Double_t[4];
2641   stat[3]             = 0.0;
2642
2643   Int_t    nbofgroups = 540;
2644   Double_t *weight    = new Double_t[nbofgroups];
2645   Int_t    *nonul     = new Int_t[nbofgroups];
2646
2647   for(Int_t k = 0; k < nbofgroups; k++){
2648     if(fEntriesLinearFitter[k] > 0) {
2649       weight[k] = 1.0;
2650       nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
2651       stat[3]++;
2652     }
2653     else weight[k] = 0.0;
2654   }
2655   stat[0]          = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
2656   stat[1]          = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
2657   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
2658
2659   delete [] weight;
2660   delete [] nonul;
2661
2662   return stat;
2663
2664 }
2665 //////////////////////////////////////////////////////////////////////////////////////
2666 // Create Histos
2667 //////////////////////////////////////////////////////////////////////////////////////
2668 //_____________________________________________________________________________
2669 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
2670 {
2671   //
2672   // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
2673   // If fNgroupprf is zero then no binning in tnp
2674   //
2675
2676   TString name("Nz");
2677   name += fCalibraMode->GetNz(2);
2678   name += "Nrphi";
2679   name += fCalibraMode->GetNrphi(2);
2680   name += "Ngp";
2681   name += fNgroupprf;
2682
2683   if(fNgroupprf != 0){
2684
2685     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2686                             ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
2687     fPRF2d->SetYTitle("Det/pad groups");
2688     fPRF2d->SetXTitle("Position x/W [pad width units]");
2689     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2690     fPRF2d->SetStats(0);
2691   }
2692   else{
2693     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2694                             ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
2695     fPRF2d->SetYTitle("Det/pad groups");
2696     fPRF2d->SetXTitle("Position x/W [pad width units]");
2697     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2698     fPRF2d->SetStats(0);
2699   }
2700
2701 }
2702
2703 //_____________________________________________________________________________
2704 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
2705 {
2706   //
2707   // Create the 2D histos
2708   //
2709
2710   TString name("Ver");
2711   name += fVersionVdriftUsed;
2712   name += "Subver";
2713   name += fSubVersionVdriftUsed;
2714   name += "FirstRun";
2715   name += fFirstRunVdrift;
2716   name += "Nz";
2717   name += fCalibraMode->GetNz(1);
2718   name += "Nrphi";
2719   name += fCalibraMode->GetNrphi(1);
2720
2721   fPH2d = new TProfile2D("PH2d",(const Char_t *) name
2722                          ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
2723                          ,nn,0,nn);
2724   fPH2d->SetYTitle("Det/pad groups");
2725   fPH2d->SetXTitle("time [#mus]");
2726   fPH2d->SetZTitle("<PH> [a.u.]");
2727   fPH2d->SetStats(0);
2728
2729 }
2730 //_____________________________________________________________________________
2731 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
2732 {
2733   //
2734   // Create the 2D histos
2735   //
2736
2737   TString name("Ver");
2738   name += fVersionGainUsed;
2739   name += "Subver";
2740   name += fSubVersionGainUsed;
2741   name += "FirstRun";
2742   name += fFirstRunGain;
2743   name += "Nz";
2744   name += fCalibraMode->GetNz(0);
2745   name += "Nrphi";
2746   name += fCalibraMode->GetNrphi(0);
2747
2748   fCH2d = new TH2I("CH2d",(const Char_t *) name
2749                    ,(Int_t)fNumberBinCharge,0,fRangeHistoCharge,nn,0,nn);
2750   fCH2d->SetYTitle("Det/pad groups");
2751   fCH2d->SetXTitle("charge deposit [a.u]");
2752   fCH2d->SetZTitle("counts");
2753   fCH2d->SetStats(0);
2754   fCH2d->Sumw2();
2755
2756 }
2757 //////////////////////////////////////////////////////////////////////////////////
2758 // Set relative scale
2759 /////////////////////////////////////////////////////////////////////////////////
2760 //_____________________________________________________________________________
2761 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
2762 {
2763   //
2764   // Set the factor that will divide the deposited charge
2765   // to fit in the histo range [0,fRangeHistoCharge]
2766   //
2767
2768   if (RelativeScale > 0.0) {
2769     fRelativeScale = RelativeScale;
2770   }
2771   else {
2772     AliInfo("RelativeScale must be strict positif!");
2773   }
2774
2775 }
2776 //////////////////////////////////////////////////////////////////////////////////
2777 // Quick way to fill a histo
2778 //////////////////////////////////////////////////////////////////////////////////
2779 //_____________________________________________________________________
2780 void  AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
2781 {
2782   //
2783   // FillCH2d: Marian style
2784   //
2785
2786   //skip simply the value out of range
2787   if((y>=fRangeHistoCharge) || (y<0.0)) return;
2788   if(fRangeHistoCharge < 0.0) return;
2789
2790   //Calcul the y place
2791   Int_t yplace = (Int_t) (fNumberBinCharge*y/fRangeHistoCharge)+1;
2792   Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
2793
2794   //Fill
2795   fCH2d->GetArray()[place]++;
2796
2797 }
2798
2799 //////////////////////////////////////////////////////////////////////////////////
2800 // Geometrical functions
2801 ///////////////////////////////////////////////////////////////////////////////////
2802 //_____________________________________________________________________________
2803 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
2804 {
2805   //
2806   // Reconstruct the layer number from the detector number
2807   //
2808
2809   return ((Int_t) (d % 6));
2810
2811 }
2812
2813 //_____________________________________________________________________________
2814 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
2815 {
2816   //
2817   // Reconstruct the stack number from the detector number
2818   //
2819   const Int_t kNlayer = 6;
2820
2821   return ((Int_t) (d % 30) / kNlayer);
2822
2823 }
2824
2825 //_____________________________________________________________________________
2826 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2827 {
2828   //
2829   // Reconstruct the sector number from the detector number
2830   //
2831   Int_t fg = 30;
2832
2833   return ((Int_t) (d / fg));
2834
2835 }
2836 ///////////////////////////////////////////////////////////////////////////////////
2837 // Getter functions for DAQ of the CH2d and the PH2d
2838 //////////////////////////////////////////////////////////////////////////////////
2839 //_____________________________________________________________________
2840 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
2841 {
2842     //
2843     // return pointer to fPH2d TProfile2D
2844     // create a new TProfile2D if it doesn't exist allready
2845     //
2846     if ( fPH2d )
2847         return fPH2d;
2848
2849     // Some parameters
2850     fTimeMax = nbtimebin;
2851     fSf      = samplefrequency;
2852
2853     CreatePH2d(540);
2854
2855     return fPH2d;
2856 }
2857 //_____________________________________________________________________
2858 TH2I* AliTRDCalibraFillHisto::GetCH2d()
2859 {
2860     //
2861     // return pointer to fCH2d TH2I
2862     // create a new TH2I if it doesn't exist allready
2863     //
2864     if ( fCH2d )
2865         return fCH2d;
2866
2867     CreateCH2d(540);
2868
2869     return fCH2d;
2870 }
2871 ////////////////////////////////////////////////////////////////////////////////////////////
2872 // Drift velocity calibration
2873 ///////////////////////////////////////////////////////////////////////////////////////////
2874 //_____________________________________________________________________
2875 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
2876 {
2877     //
2878     // return pointer to TLinearFitter Calibration
2879     // if force is true create a new TLinearFitter if it doesn't exist allready
2880     //
2881
2882   if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
2883     return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
2884   }
2885
2886   // if we are forced and TLinearFitter doesn't yet exist create it
2887
2888   // new TLinearFitter
2889   TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
2890   fLinearFitterArray.AddAt(linearfitter,detector);
2891   return linearfitter;
2892 }
2893
2894 //____________________________________________________________________________
2895 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
2896 {
2897   //
2898   // Analyse array of linear fitter because can not be written
2899   // Store two arrays: one with the param the other one with the error param + number of entries
2900   //
2901
2902   for(Int_t k = 0; k < 540; k++){
2903     TLinearFitter *linearfitter = GetLinearFitter(k);
2904     if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
2905       TVectorD  *par  = new TVectorD(2);
2906       TVectorD   pare = TVectorD(2);
2907       TVectorD  *parE = new TVectorD(3);
2908       if((linearfitter->EvalRobust(0.8)==0)) {
2909         //linearfitter->Eval();
2910         linearfitter->GetParameters(*par);
2911         //linearfitter->GetErrors(pare);
2912         //Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
2913         //(*parE)[0] = pare[0]*ppointError;
2914         //(*parE)[1] = pare[1]*ppointError;
2915
2916         (*parE)[0] = 0.0;
2917         (*parE)[1] = 0.0;
2918         (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
2919         ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
2920         ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);
2921       }
2922     }
2923   }
2924 }
2925
2926
2927