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