TRD/AliTRDCalibraFillHisto.cxx

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