TRD/AliTRDCalibraFillHisto.cxx

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