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