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