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