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     if (time>=0 && time<fTimeMax) {
1695       fPHPlace[time] = group[1];
1696       fPHValue[time] = charge;
1697     }
1698   }
1699
1700   return correction;
1701
1702 }
1703 //____________Offine tracking in the AliTRDtracker_____________________________
1704 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1705 {
1706   //
1707   // Reset values per tracklet
1708   //
1709
1710   //Reset good tracklet
1711   fGoodTracklet = kTRUE;
1712
1713   // Reset the fPHValue
1714   if (fPH2dOn) {
1715     //Reset the fPHValue and fPHPlace
1716     for (Int_t k = 0; k < fTimeMax; k++) {
1717       fPHValue[k] = 0.0;
1718       fPHPlace[k] = -1;
1719     }
1720   }
1721
1722   // Reset the fAmpTotal where we put value
1723   if (fCH2dOn) {
1724     for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1725       fAmpTotal[k] = 0.0;
1726     }
1727   }
1728 }
1729 //____________Offine tracking in the AliTRDtracker_____________________________
1730 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1731 {
1732   //
1733   // For the offline tracking or mcm tracklets
1734   // This function will be called in the functions UpdateHistogram...
1735   // to fill the info of a track for the relativ gain calibration
1736   //
1737
1738   Int_t nb            =  0;   // Nombre de zones traversees
1739   Int_t fd            = -1;   // Premiere zone non nulle
1740   Float_t totalcharge = 0.0;  // Total charge for the supermodule histo
1741
1742   //printf("CH2d nbclusters %d, fNumberClusters %d, fNumberClustersf %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1743
1744   if(nbclusters < fNumberClusters) return;
1745   if(nbclusters > fNumberClustersf) return;
1746
1747
1748   // Normalize with the number of clusters
1749   Double_t normalizeCst = fRelativeScale;
1750   if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1751
1752   //printf("Number of groups in one detector %d\n",fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0));
1753
1754   // See if the track goes through different zones
1755   for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1756     //printf("fAmpTotal %f for %d\n",fAmpTotal[k],k);
1757     if (fAmpTotal[k] > 0.0) {
1758       totalcharge += fAmpTotal[k];
1759       nb++;
1760       if (nb == 1) {
1761         fd = k;
1762       }
1763     }
1764   }
1765
1766   //printf("CH2d: nb %d, fd %d, calibration group %d, amplitude %f, detector %d\n",nb,fd,fCalibraMode->GetXbins(0),fAmpTotal[fd]/normalizeCst,fDetectorPreviousTrack);
1767
1768   switch (nb)
1769     {
1770     case 1:
1771       fNumberUsedCh[0]++;
1772       fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1773       if (fHisto2d) {
1774         FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1775         //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1776       }
1777       if (fVector2d) {
1778         fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1779       }
1780       break;
1781     case 2:
1782       if ((fAmpTotal[fd]   > 0.0) &&
1783           (fAmpTotal[fd+1] > 0.0)) {
1784         // One of the two very big
1785         if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1786           if (fHisto2d) {
1787             FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1788             //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1789           }
1790           if (fVector2d) {
1791             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1792           }
1793           fNumberUsedCh[1]++;
1794           fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1795         }
1796         if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1797           if (fHisto2d) {
1798             FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
1799             //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
1800           }
1801           if (fVector2d) {
1802             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
1803           }
1804           fNumberUsedCh[1]++;
1805           fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1806         }
1807       }
1808       if (fCalibraMode->GetNfragZ(0) > 1) {
1809         if (fAmpTotal[fd] > 0.0) {
1810           if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1811             if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1812               // One of the two very big
1813               if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1814                 if (fHisto2d) {
1815                   FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1816                   //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1817                 }
1818                 if (fVector2d) {
1819                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1820                 }
1821                 fNumberUsedCh[1]++;
1822                 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1823               }
1824               if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1825                 if (fHisto2d) {
1826                   FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1827                   //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1828                 }
1829                 fNumberUsedCh[1]++;
1830                 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1831                 if (fVector2d) {
1832                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1833                 }
1834               }
1835             }
1836           }
1837         }
1838       }
1839       break;
1840     default: break;
1841   }
1842 }
1843 //____________Offine tracking in the AliTRDtracker_____________________________
1844 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1845 {
1846   //
1847   // For the offline tracking or mcm tracklets
1848   // This function will be called in the functions UpdateHistogram...
1849   // to fill the info of a track for the drift velocity  calibration
1850   //
1851
1852   Int_t nb  =  1; // Nombre de zones traversees 1, 2 ou plus de 3
1853   Int_t fd1 = -1; // Premiere zone non nulle
1854   Int_t fd2 = -1; // Deuxieme zone non nulle
1855   Int_t k1  = -1; // Debut de la premiere zone
1856   Int_t k2  = -1; // Debut de la seconde zone
1857   Int_t nbclusters = 0; // number of clusters
1858
1859
1860
1861   // See if the track goes through different zones
1862   for (Int_t k = 0; k < fTimeMax; k++) {
1863     if (fPHValue[k] > 0.0) {
1864       nbclusters++;
1865       if (fd1 == -1) {
1866         fd1 = fPHPlace[k];
1867         k1  = k;
1868       }
1869       if (fPHPlace[k] != fd1) {
1870         if (fd2 == -1) {
1871           k2  = k;
1872           fd2 = fPHPlace[k];
1873           nb  = 2;
1874         }
1875         if (fPHPlace[k] != fd2) {
1876           nb = 3;
1877         }
1878       }
1879     }
1880   }
1881
1882   // See if noise before and after
1883   if(fMaxCluster > 0) {
1884     if(fPHValue[0] > fMaxCluster) return;
1885     if(fTimeMax > fNbMaxCluster) {
1886       for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
1887         if(fPHValue[k] > fMaxCluster) return;
1888       }
1889     }
1890   }
1891
1892   //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1893
1894   if(nbclusters < fNumberClusters) return;
1895   if(nbclusters > fNumberClustersf) return;
1896
1897   switch(nb)
1898     {
1899     case 1:
1900       fNumberUsedPh[0]++;
1901       for (Int_t i = 0; i < fTimeMax; i++) {
1902         if (fHisto2d) {
1903           if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1904           else {
1905             if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1906               }
1907           //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
1908         }
1909         if (fVector2d) {
1910           if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1911           else {
1912             if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1913           }
1914         }
1915       }
1916       break;
1917     case 2:
1918       if ((fd1 == fd2+1) ||
1919           (fd2 == fd1+1)) {
1920         // One of the two fast all the think
1921         if (k2 > (k1+fDifference)) {
1922           //we choose to fill the fd1 with all the values
1923           fNumberUsedPh[1]++;
1924           for (Int_t i = 0; i < fTimeMax; i++) {
1925             if (fHisto2d) {
1926               if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1927               else {
1928                 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1929                   }
1930             }
1931             if (fVector2d) {
1932               if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1933               else {
1934                 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1935                   }
1936             }
1937           }
1938         }
1939         if ((k2+fDifference) < fTimeMax) {
1940           //we choose to fill the fd2 with all the values
1941           fNumberUsedPh[1]++;
1942           for (Int_t i = 0; i < fTimeMax; i++) {
1943             if (fHisto2d) {
1944               if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1945               else {
1946                 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1947               }
1948             }
1949           if (fVector2d) {
1950             if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1951             else {
1952               if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1953             }
1954           }
1955           }
1956         }
1957       }
1958       // Two zones voisines sinon rien!
1959       if (fCalibraMode->GetNfragZ(1) > 1) {
1960         // Case 2
1961         if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1962           if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1963             // One of the two fast all the think
1964             if (k2 > (k1+fDifference)) {
1965               //we choose to fill the fd1 with all the values
1966               fNumberUsedPh[1]++;
1967               for (Int_t i = 0; i < fTimeMax; i++) {
1968                 if (fHisto2d) {
1969                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1970                   else {
1971                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1972                   }
1973                 }
1974                 if (fVector2d) {
1975                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1976                   else {
1977                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1978                   }
1979                 }
1980               }
1981             }
1982             if ((k2+fDifference) < fTimeMax) {
1983               //we choose to fill the fd2 with all the values
1984               fNumberUsedPh[1]++;
1985               for (Int_t i = 0; i < fTimeMax; i++) {
1986                 if (fHisto2d) {
1987                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1988                   else {
1989                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1990                   }
1991                 }
1992                 if (fVector2d) {
1993                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1994                   else {
1995                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1996                   }
1997                 }
1998               }
1999             }
2000           }
2001         }
2002         // Two zones voisines sinon rien!
2003         // Case 3
2004         if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2005           if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2006             // One of the two fast all the think
2007             if (k2 > (k1 + fDifference)) {
2008               //we choose to fill the fd1 with all the values
2009               fNumberUsedPh[1]++;
2010               for (Int_t i = 0; i < fTimeMax; i++) {
2011                 if (fHisto2d) {
2012                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2013                   else {
2014                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2015                   }
2016                 }
2017                 if (fVector2d) {
2018                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2019                   else {
2020                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2021                   }
2022                 }
2023               }
2024             }
2025             if ((k2+fDifference) < fTimeMax) {
2026               //we choose to fill the fd2 with all the values
2027               fNumberUsedPh[1]++;
2028               for (Int_t i = 0; i < fTimeMax; i++) {
2029                 if (fHisto2d) {
2030                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2031                   else {
2032                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2033                   }
2034                 }
2035                 if (fVector2d) {
2036                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2037                   else {
2038                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2039                   }
2040                 }
2041               }
2042             }
2043           }
2044         }
2045       }
2046       break;
2047     default: break;
2048     }
2049 }
2050 //////////////////////////////////////////////////////////////////////////////////////////
2051 // DAQ process functions
2052 /////////////////////////////////////////////////////////////////////////////////////////
2053 //_____________________________________________________________________
2054 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader)
2055  { //main
2056   //
2057   // Event Processing loop - AliTRDrawStream
2058   //
2059   // 0 timebin problem
2060   // 1 no input
2061   // 2 input
2062   // Same algorithm as TestBeam but different reader
2063   //
2064
2065   AliTRDrawStream *rawStream = new AliTRDrawStream(rawReader);
2066
2067   AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
2068   digitsManager->CreateArrays();
2069
2070   Int_t withInput = 1;
2071
2072   Double_t phvalue[16][144][36];
2073   for(Int_t k = 0; k < 36; k++){
2074     for(Int_t j = 0; j < 16; j++){
2075       for(Int_t c = 0; c < 144; c++){
2076         phvalue[j][c][k] = 0.0;
2077       }
2078     }
2079   }
2080
2081   fDetectorPreviousTrack = -1;
2082   fMCMPrevious           = -1;
2083   fROBPrevious           = -1;
2084
2085   Int_t nbtimebin = 0;
2086   Int_t baseline  = 10;
2087
2088
2089     fTimeMax = 0;
2090
2091     Int_t det    = 0;
2092     while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
2093
2094       if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
2095         //      printf("there is ADC data on this chamber!\n");
2096
2097         AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
2098         if (digits->HasData()) { //array
2099
2100           AliTRDSignalIndex   *indexes = digitsManager->GetIndexes(det);
2101           if (indexes->IsAllocated() == kFALSE) {
2102             AliError("Indexes do not exist!");
2103             break;
2104           }
2105           Int_t iRow  = 0;
2106           Int_t iCol  = 0;
2107           indexes->ResetCounters();
2108
2109           while (indexes->NextRCIndex(iRow, iCol)) { //column,row
2110             //printf(" det %d \t row %d \t col %d \t digit\n",det,iRow,iCol);
2111             //while (rawStream->Next()) {
2112
2113             Int_t idetector = det;                                             //  current detector
2114             //Int_t imcm      = rawStream->GetMCM();                            //  current MCM
2115             //Int_t irob      = rawStream->GetROB();                            //  current ROB
2116
2117
2118             if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)) {
2119               // Fill
2120               withInput = TMath::Max(FillDAQ(phvalue),withInput);
2121
2122               // reset
2123               for(Int_t k = 0; k < 36; k++){
2124                 for(Int_t j = 0; j < 16; j++){
2125                   for(Int_t c = 0; c < 144; c++){
2126                     phvalue[j][c][k] = 0.0;
2127                   }
2128                 }
2129               }
2130             }
2131
2132             fDetectorPreviousTrack = idetector;
2133             //fMCMPrevious           = imcm;
2134             //fROBPrevious           = irob;
2135
2136             //    nbtimebin              = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
2137             AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
2138             nbtimebin              = digitParam->GetNTimeBins(det);              //  number of time bins read from data
2139             baseline               = digitParam->GetADCbaseline(det);            //  baseline
2140
2141             if(nbtimebin == 0) return 0;
2142             if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2143             fTimeMax          = nbtimebin;
2144
2145             fNumberClustersf    = fTimeMax;
2146             fNumberClusters     = (Int_t)(fNumberClustersProcent*fTimeMax);
2147
2148
2149             for(Int_t itime = 0; itime < nbtimebin; itime++) {
2150               //            phvalue[row][col][itime] = signal[itime]-baseline;
2151               phvalue[iRow][iCol][itime] = (Short_t)(digits->GetData(iRow,iCol,itime) - baseline);
2152               /*if(phvalue[iRow][iCol][itime] >= 20) {
2153                  printf("----------> phvalue[%d][%d][%d] %d  baseline %d \n",
2154                        iRow,
2155                        iCol,
2156                        itime,
2157                        (Short_t)(digits->GetData(iRow,iCol,itime)),
2158                        baseline);
2159                        }*/
2160             }
2161
2162           }//column,row
2163
2164           // fill the last one
2165           if(fDetectorPreviousTrack != -1){
2166
2167             // Fill
2168             withInput = TMath::Max(FillDAQ(phvalue),withInput);
2169             //      printf("\n ---> withinput %d\n\n",withInput);
2170             // reset
2171             for(Int_t k = 0; k < 36; k++){
2172               for(Int_t j = 0; j < 16; j++){
2173                 for(Int_t c = 0; c < 144; c++){
2174                   phvalue[j][c][k] = 0.0;
2175                 }
2176               }
2177             }
2178           }
2179
2180         }//array
2181       }//QA
2182       digitsManager->ClearArrays(det);
2183     }//idetector
2184     delete digitsManager;
2185
2186     delete rawStream;
2187     return withInput;
2188  }//main
2189 //_____________________________________________________________________
2190 //////////////////////////////////////////////////////////////////////////////
2191 // Routine inside the DAQ process
2192 /////////////////////////////////////////////////////////////////////////////
2193 //_______________________________________________________________________
2194 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2195
2196   //
2197   // Look for the maximum by collapsing over the time
2198   // Sum over four pad col and two pad row
2199   //
2200
2201   Int_t used = 0;
2202
2203
2204   Int_t idect = fDetectorPreviousTrack;
2205   //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2206   Double_t sum[36];
2207   for(Int_t tb = 0; tb < 36; tb++){
2208     sum[tb] = 0.0;
2209   }
2210
2211   //fGoodTracklet = kTRUE;
2212   //fDetectorPreviousTrack = 0;
2213
2214
2215   ///////////////////////////
2216   // look for maximum
2217   /////////////////////////
2218
2219   Int_t imaxRow = 0;
2220   Int_t imaxCol = 0;
2221   Double_t integralMax = -1;
2222
2223   for (Int_t ir = 1; ir <= 15; ir++)
2224     {
2225       for (Int_t ic = 2; ic <= 142; ic++)
2226         {
2227           Double_t integral = 0;
2228           for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2229             {
2230               for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2231                 {
2232                   if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2233                       ic + ishiftC >= 1 && ic + ishiftC <= 144)
2234                     {
2235
2236                       for(Int_t tb = 0; tb< fTimeMax; tb++){
2237                         integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2238                       }// addtb
2239                     } //addsignal
2240                 } //shiftC
2241             } // shiftR
2242           if (integralMax < integral)
2243             {
2244               imaxRow = ir;
2245               imaxCol = ic;
2246               integralMax = integral;
2247
2248             } // check max integral
2249         } //ic
2250     } // ir
2251
2252   //  printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2253   //if((imaxRow == 0) || (imaxRow >= 15) || (imaxCol <= 3) || (imaxCol >= 140)) {
2254   //  used=1;
2255   //  return used;
2256   // }
2257
2258   if(((imaxRow + fNumberRowDAQ) > 16) || (imaxRow == 0) || ((imaxCol - fNumberColDAQ) <= 1) || ((imaxCol + fNumberColDAQ) >= 144)) {
2259     used=1;
2260     return used;
2261   }
2262   //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2263   //if(!fGoodTracklet) used = 1;;
2264
2265   //  /////////////////////////////////////////////////////
2266   // sum ober 2 row and 4 pad cols for each time bins
2267   //  ////////////////////////////////////////////////////
2268
2269
2270
2271   for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2272     {
2273       for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2274         {
2275           if (imaxRow + ishiftR >= 1 && imaxRow + ishiftR <= 16 &&
2276               imaxCol + ishiftC >= 1 && imaxCol + ishiftC <= 144)
2277             {
2278               for(Int_t it = 0; it < fTimeMax; it++){
2279                 sum[it] += phvalue[imaxRow + ishiftR-1][imaxCol + ishiftC-1][it];
2280               }
2281             }
2282         } // col shift
2283     }// row shift
2284
2285   Int_t nbcl = 0;
2286   Double_t sumcharge = 0.0;
2287   for(Int_t it = 0; it < fTimeMax; it++){
2288     sumcharge += sum[it];
2289     if(sum[it] > fThresholdClustersDAQ)  nbcl++;
2290   }
2291
2292
2293   /////////////////////////////////////////////////////////
2294   // Debug
2295   ////////////////////////////////////////////////////////
2296   if(fDebugLevel > 1){
2297     if ( !fDebugStreamer ) {
2298       //debug stream
2299       TDirectory *backup = gDirectory;
2300       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2301       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2302     }
2303
2304     Double_t amph0 = sum[0];
2305     Double_t amphlast = sum[fTimeMax-1];
2306     Double_t rms      = TMath::RMS(fTimeMax,sum);
2307     Int_t    goodtracklet = (Int_t) fGoodTracklet;
2308     for(Int_t it = 0; it < fTimeMax; it++){
2309       Double_t clustera = sum[it];
2310
2311     (* fDebugStreamer) << "FillDAQa"<<
2312       "ampTotal="<<sumcharge<<
2313       "row="<<imaxRow<<
2314       "col="<<imaxCol<<
2315       "detector="<<idect<<
2316       "amph0="<<amph0<<
2317       "amphlast="<<amphlast<<
2318       "goodtracklet="<<goodtracklet<<
2319       "clustera="<<clustera<<
2320       "it="<<it<<
2321       "rms="<<rms<<
2322       "nbcl="<<nbcl<<
2323       "\n";
2324     }
2325   }
2326
2327   ////////////////////////////////////////////////////////
2328   // fill
2329   ///////////////////////////////////////////////////////
2330   //printf("fNumberClusters %d, fNumberClustersf %d\n",fNumberClusters,fNumberClustersf);
2331   if(sum[0] > 100.0) used = 1;
2332   if(nbcl < fNumberClusters) used = 1;
2333   if(nbcl > fNumberClustersf) used = 1;
2334
2335   //if(fDetectorPreviousTrack == 15){
2336   //  printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2337   //}
2338   //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2339   if(used == 0){
2340     for(Int_t it = 0; it < fTimeMax; it++){
2341       if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2342       else{
2343         if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2344       }
2345       //if(fFillWithZero) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2346       //else{
2347       // if(sum[it] > 0.0) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2348       //}
2349     }
2350
2351
2352     //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2353     used = 2;
2354     //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2355
2356   }
2357
2358   return used;
2359
2360 }
2361 //____________Online trackling in AliTRDtrigger________________________________
2362 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2363 {
2364   //
2365   // For the DAQ
2366   // Fill a simple average pulse height
2367   //
2368
2369
2370   ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2371
2372
2373   return kTRUE;
2374
2375 }
2376 //____________Write_____________________________________________________
2377 //_____________________________________________________________________
2378 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2379 {
2380   //
2381   //  Write infos to file
2382   //
2383
2384   //For debugging
2385   if ( fDebugStreamer ) {
2386     delete fDebugStreamer;
2387     fDebugStreamer = 0x0;
2388   }
2389
2390   AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2391                ,fNumberTrack
2392                ,fNumberUsedCh[0]
2393                ,fNumberUsedCh[1]
2394                ,fNumberUsedPh[0]
2395                ,fNumberUsedPh[1]));
2396
2397   TDirectory *backup = gDirectory;
2398   TString option;
2399
2400   if ( append )
2401     option = "update";
2402   else
2403     option = "recreate";
2404
2405   TFile f(filename,option.Data());
2406
2407   TStopwatch stopwatch;
2408   stopwatch.Start();
2409   if(fVector2d) {
2410     f.WriteTObject(fCalibraVector);
2411   }
2412
2413   if (fCH2dOn ) {
2414     if (fHisto2d) {
2415       f.WriteTObject(fCH2d);
2416     }
2417   }
2418   if (fPH2dOn ) {
2419     if (fHisto2d) {
2420       f.WriteTObject(fPH2d);
2421     }
2422   }
2423   if (fPRF2dOn) {
2424     if (fHisto2d) {
2425         f.WriteTObject(fPRF2d);
2426     }
2427   }
2428   if(fLinearFitterOn){
2429     if(fLinearFitterDebugOn) AnalyseLinearFitter();
2430     f.WriteTObject(fLinearVdriftFit);
2431   }
2432
2433   f.Close();
2434
2435   if ( backup ) backup->cd();
2436
2437   AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2438                ,stopwatch.RealTime(),stopwatch.CpuTime()));
2439 }
2440 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2441 // Stats stuff
2442 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2443 //___________________________________________probe the histos__________________________________________________
2444 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2445 {
2446   //
2447   // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2448   // debug mode with 2 for TH2I and 3 for TProfile2D
2449   // It gives a pointer to a Double_t[7] with the info following...
2450   // [0] : number of calibration groups with entries
2451   // [1] : minimal number of entries found
2452   // [2] : calibration group number of the min
2453   // [3] : maximal number of entries found
2454   // [4] : calibration group number of the max
2455   // [5] : mean number of entries found
2456   // [6] : mean relative error
2457   //
2458
2459   Double_t *info = new Double_t[7];
2460
2461   // Number of Xbins (detectors or groups of pads)
2462   Int_t    nbins   = h->GetNbinsY(); //number of calibration groups
2463   Int_t    nxbins  = h->GetNbinsX(); //number of bins per histo
2464
2465   // Initialise
2466   Double_t nbwe = 0; //number of calibration groups with entries
2467   Double_t minentries = 0; //minimal number of entries found
2468   Double_t maxentries = 0; //maximal number of entries found
2469   Double_t placemin = 0; //calibration group number of the min
2470   Double_t placemax = -1; //calibration group number of the max
2471   Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2472   Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2473
2474   Double_t counter = 0;
2475
2476   //Debug
2477   TH1F *nbEntries = 0x0;//distribution of the number of entries
2478   TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2479   TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2480
2481   // Beginning of the loop over the calibration groups
2482   for (Int_t idect = 0; idect < nbins; idect++) {
2483
2484     TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2485     projch->SetDirectory(0);
2486
2487     // Number of entries for this calibration group
2488     Double_t nentries = 0.0;
2489     if((i%2) == 0){
2490       for (Int_t k = 0; k < nxbins; k++) {
2491         nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2492       }
2493     }
2494     else{
2495       for (Int_t k = 0; k < nxbins; k++) {
2496         nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2497         if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2498           meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
2499           counter++;
2500         }
2501       }
2502     }
2503
2504     //Debug
2505     if(i > 1){
2506       if(nentries > 0){
2507         if(!((Bool_t)nbEntries)) nbEntries = new TH1F("Number of entries","Number of entries",100,(Int_t)nentries/2,nentries*2);
2508         nbEntries->SetDirectory(0);
2509         nbEntries->Fill(nentries);
2510         if(!((Bool_t)nbEntriesPerGroup)) nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group",nbins,0,nbins);
2511         nbEntriesPerGroup->SetDirectory(0);
2512         nbEntriesPerGroup->Fill(idect+0.5,nentries);
2513         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));
2514         nbEntriesPerSp->SetDirectory(0);
2515         nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
2516       }
2517     }
2518
2519     //min amd max
2520     if(nentries > maxentries){
2521       maxentries = nentries;
2522       placemax = idect;
2523     }
2524     if(idect == 0) {
2525       minentries = nentries;
2526     }
2527     if(nentries < minentries){
2528       minentries = nentries;
2529       placemin = idect;
2530     }
2531     //nbwe
2532     if(nentries > 0) {
2533       nbwe++;
2534       meanstats += nentries;
2535     }
2536   }//calibration groups loop
2537
2538   if(nbwe > 0) meanstats /= nbwe;
2539   if(counter > 0) meanrelativerror /= counter;
2540
2541   AliInfo(Form("There are %f calibration groups with entries",nbwe));
2542   AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
2543   AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
2544   AliInfo(Form("The mean number of entries is %f",meanstats));
2545   if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
2546
2547   info[0] = nbwe;
2548   info[1] = minentries;
2549   info[2] = placemin;
2550   info[3] = maxentries;
2551   info[4] = placemax;
2552   info[5] = meanstats;
2553   info[6] = meanrelativerror;
2554
2555   if(nbEntries && nbEntriesPerSp && nbEntriesPerGroup){
2556     gStyle->SetPalette(1);
2557     gStyle->SetOptStat(1111);
2558     gStyle->SetPadBorderMode(0);
2559     gStyle->SetCanvasColor(10);
2560     gStyle->SetPadLeftMargin(0.13);
2561     gStyle->SetPadRightMargin(0.01);
2562     TCanvas *stat = new TCanvas("stat","",50,50,600,800);
2563     stat->Divide(2,1);
2564     stat->cd(1);
2565     nbEntries->Draw("");
2566     stat->cd(2);
2567     nbEntriesPerSp->SetStats(0);
2568     nbEntriesPerSp->Draw("");
2569     TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
2570     stat1->cd();
2571     nbEntriesPerGroup->SetStats(0);
2572     nbEntriesPerGroup->Draw("");
2573   }
2574
2575   return info;
2576
2577 }
2578 //____________________________________________________________________________
2579 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
2580 {
2581   //
2582   // Return a Int_t[4] with:
2583   // 0 Mean number of entries
2584   // 1 median of number of entries
2585   // 2 rms of number of entries
2586   // 3 number of group with entries
2587   //
2588
2589   Double_t *stat = new Double_t[4];
2590   stat[3]             = 0.0;
2591
2592   Int_t    nbofgroups = CalculateTotalNumberOfBins(0);
2593
2594   Double_t *weight = new Double_t[nbofgroups];
2595   Double_t *nonul = new Double_t[nbofgroups];
2596
2597   for(Int_t k = 0; k < nbofgroups; k++){
2598     if(fEntriesCH[k] > 0) {
2599       weight[k] = 1.0;
2600       nonul[(Int_t)stat[3]] = fEntriesCH[k];
2601       stat[3]++;
2602     }
2603     else weight[k] = 0.0;
2604   }
2605   stat[0]          = TMath::Mean(nbofgroups,fEntriesCH,weight);
2606   stat[1]          = TMath::Median(nbofgroups,fEntriesCH,weight);
2607   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
2608
2609   delete [] weight;
2610   delete [] nonul;
2611
2612   return stat;
2613
2614 }
2615 //____________________________________________________________________________
2616 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
2617 {
2618   //
2619   // Return a Int_t[4] with:
2620   // 0 Mean number of entries
2621   // 1 median of number of entries
2622   // 2 rms of number of entries
2623   // 3 number of group with entries
2624   //
2625
2626   Double_t *stat      = new Double_t[4];
2627   stat[3]             = 0.0;
2628
2629   Int_t    nbofgroups = 540;
2630   Double_t *weight    = new Double_t[nbofgroups];
2631   Int_t    *nonul     = new Int_t[nbofgroups];
2632
2633   for(Int_t k = 0; k < nbofgroups; k++){
2634     if(fEntriesLinearFitter[k] > 0) {
2635       weight[k] = 1.0;
2636       nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
2637       stat[3]++;
2638     }
2639     else weight[k] = 0.0;
2640   }
2641   stat[0]          = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
2642   stat[1]          = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
2643   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
2644
2645   delete [] weight;
2646   delete [] nonul;
2647
2648   return stat;
2649
2650 }
2651 //////////////////////////////////////////////////////////////////////////////////////
2652 // Create Histos
2653 //////////////////////////////////////////////////////////////////////////////////////
2654 //_____________________________________________________________________________
2655 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
2656 {
2657   //
2658   // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
2659   // If fNgroupprf is zero then no binning in tnp
2660   //
2661
2662   TString name("Nz");
2663   name += fCalibraMode->GetNz(2);
2664   name += "Nrphi";
2665   name += fCalibraMode->GetNrphi(2);
2666   name += "Ngp";
2667   name += fNgroupprf;
2668
2669   if(fNgroupprf != 0){
2670
2671     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2672                             ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
2673     fPRF2d->SetYTitle("Det/pad groups");
2674     fPRF2d->SetXTitle("Position x/W [pad width units]");
2675     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2676     fPRF2d->SetStats(0);
2677   }
2678   else{
2679     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2680                             ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
2681     fPRF2d->SetYTitle("Det/pad groups");
2682     fPRF2d->SetXTitle("Position x/W [pad width units]");
2683     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2684     fPRF2d->SetStats(0);
2685   }
2686
2687 }
2688
2689 //_____________________________________________________________________________
2690 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
2691 {
2692   //
2693   // Create the 2D histos
2694   //
2695
2696   TString name("Ver");
2697   name += fVersionVdriftUsed;
2698   name += "Subver";
2699   name += fSubVersionVdriftUsed;
2700   name += "FirstRun";
2701   name += fFirstRunVdrift;
2702   name += "Nz";
2703   name += fCalibraMode->GetNz(1);
2704   name += "Nrphi";
2705   name += fCalibraMode->GetNrphi(1);
2706
2707   fPH2d = new TProfile2D("PH2d",(const Char_t *) name
2708                          ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
2709                          ,nn,0,nn);
2710   fPH2d->SetYTitle("Det/pad groups");
2711   fPH2d->SetXTitle("time [#mus]");
2712   fPH2d->SetZTitle("<PH> [a.u.]");
2713   fPH2d->SetStats(0);
2714
2715 }
2716 //_____________________________________________________________________________
2717 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
2718 {
2719   //
2720   // Create the 2D histos
2721   //
2722
2723   TString name("Ver");
2724   name += fVersionGainUsed;
2725   name += "Subver";
2726   name += fSubVersionGainUsed;
2727   name += "FirstRun";
2728   name += fFirstRunGain;
2729   name += "Nz";
2730   name += fCalibraMode->GetNz(0);
2731   name += "Nrphi";
2732   name += fCalibraMode->GetNrphi(0);
2733
2734   fCH2d = new TH2I("CH2d",(const Char_t *) name
2735                    ,(Int_t)fNumberBinCharge,0,fRangeHistoCharge,nn,0,nn);
2736   fCH2d->SetYTitle("Det/pad groups");
2737   fCH2d->SetXTitle("charge deposit [a.u]");
2738   fCH2d->SetZTitle("counts");
2739   fCH2d->SetStats(0);
2740   fCH2d->Sumw2();
2741
2742 }
2743 //////////////////////////////////////////////////////////////////////////////////
2744 // Set relative scale
2745 /////////////////////////////////////////////////////////////////////////////////
2746 //_____________________________________________________________________________
2747 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
2748 {
2749   //
2750   // Set the factor that will divide the deposited charge
2751   // to fit in the histo range [0,fRangeHistoCharge]
2752   //
2753
2754   if (RelativeScale > 0.0) {
2755     fRelativeScale = RelativeScale;
2756   }
2757   else {
2758     AliInfo("RelativeScale must be strict positif!");
2759   }
2760
2761 }
2762 //////////////////////////////////////////////////////////////////////////////////
2763 // Quick way to fill a histo
2764 //////////////////////////////////////////////////////////////////////////////////
2765 //_____________________________________________________________________
2766 void  AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
2767 {
2768   //
2769   // FillCH2d: Marian style
2770   //
2771
2772   //skip simply the value out of range
2773   if((y>=fRangeHistoCharge) || (y<0.0)) return;
2774   if(fRangeHistoCharge < 0.0) return;
2775
2776   //Calcul the y place
2777   Int_t yplace = (Int_t) (fNumberBinCharge*y/fRangeHistoCharge)+1;
2778   Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
2779
2780   //Fill
2781   fCH2d->GetArray()[place]++;
2782
2783 }
2784
2785 //////////////////////////////////////////////////////////////////////////////////
2786 // Geometrical functions
2787 ///////////////////////////////////////////////////////////////////////////////////
2788 //_____________________________________________________________________________
2789 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
2790 {
2791   //
2792   // Reconstruct the layer number from the detector number
2793   //
2794
2795   return ((Int_t) (d % 6));
2796
2797 }
2798
2799 //_____________________________________________________________________________
2800 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
2801 {
2802   //
2803   // Reconstruct the stack number from the detector number
2804   //
2805   const Int_t kNlayer = 6;
2806
2807   return ((Int_t) (d % 30) / kNlayer);
2808
2809 }
2810
2811 //_____________________________________________________________________________
2812 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2813 {
2814   //
2815   // Reconstruct the sector number from the detector number
2816   //
2817   Int_t fg = 30;
2818
2819   return ((Int_t) (d / fg));
2820
2821 }
2822 ///////////////////////////////////////////////////////////////////////////////////
2823 // Getter functions for DAQ of the CH2d and the PH2d
2824 //////////////////////////////////////////////////////////////////////////////////
2825 //_____________________________________________________________________
2826 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
2827 {
2828     //
2829     // return pointer to fPH2d TProfile2D
2830     // create a new TProfile2D if it doesn't exist allready
2831     //
2832     if ( fPH2d )
2833         return fPH2d;
2834
2835     // Some parameters
2836     fTimeMax = nbtimebin;
2837     fSf      = samplefrequency;
2838
2839     CreatePH2d(540);
2840
2841     return fPH2d;
2842 }
2843 //_____________________________________________________________________
2844 TH2I* AliTRDCalibraFillHisto::GetCH2d()
2845 {
2846     //
2847     // return pointer to fCH2d TH2I
2848     // create a new TH2I if it doesn't exist allready
2849     //
2850     if ( fCH2d )
2851         return fCH2d;
2852
2853     CreateCH2d(540);
2854
2855     return fCH2d;
2856 }
2857 ////////////////////////////////////////////////////////////////////////////////////////////
2858 // Drift velocity calibration
2859 ///////////////////////////////////////////////////////////////////////////////////////////
2860 //_____________________________________________________________________
2861 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
2862 {
2863     //
2864     // return pointer to TLinearFitter Calibration
2865     // if force is true create a new TLinearFitter if it doesn't exist allready
2866     //
2867
2868   if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
2869     return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
2870   }
2871
2872   // if we are forced and TLinearFitter doesn't yet exist create it
2873
2874   // new TLinearFitter
2875   TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
2876   fLinearFitterArray.AddAt(linearfitter,detector);
2877   return linearfitter;
2878 }
2879
2880 //____________________________________________________________________________
2881 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
2882 {
2883   //
2884   // Analyse array of linear fitter because can not be written
2885   // Store two arrays: one with the param the other one with the error param + number of entries
2886   //
2887
2888   for(Int_t k = 0; k < 540; k++){
2889     TLinearFitter *linearfitter = GetLinearFitter(k);
2890     if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
2891       TVectorD  *par  = new TVectorD(2);
2892       TVectorD   pare = TVectorD(2);
2893       TVectorD  *parE = new TVectorD(3);
2894       if((linearfitter->EvalRobust(0.8)==0)) {
2895         //linearfitter->Eval();
2896         linearfitter->GetParameters(*par);
2897         //linearfitter->GetErrors(pare);
2898         //Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
2899         //(*parE)[0] = pare[0]*ppointError;
2900         //(*parE)[1] = pare[1]*ppointError;
2901
2902         (*parE)[0] = 0.0;
2903         (*parE)[1] = 0.0;
2904         (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
2905         ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
2906         ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);
2907       }
2908     }
2909   }
2910 }
2911
2912
2913