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   AliTRDrawStream *rawStream = new AliTRDrawStream(rawReader);
2677   rawStream->SetNoErrorWarning();
2678   rawStream->SetSharedPadReadout(kFALSE);
2679
2680   AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
2681   digitsManager->CreateArrays();
2682
2683   Int_t withInput = 1;
2684
2685   Double_t phvalue[16][144][36];
2686   for(Int_t k = 0; k < 36; k++){
2687     for(Int_t j = 0; j < 16; j++){
2688       for(Int_t c = 0; c < 144; c++){
2689         phvalue[j][c][k] = 0.0;
2690       }
2691     }
2692   }
2693
2694   fDetectorPreviousTrack = -1;
2695   fMCMPrevious           = -1;
2696   fROBPrevious           = -1;
2697
2698   Int_t nbtimebin = 0;
2699   Int_t baseline  = 10;
2700
2701
2702     fTimeMax = 0;
2703
2704     Int_t det    = 0;
2705     while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
2706
2707       if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
2708         //      printf("there is ADC data on this chamber!\n");
2709
2710         AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
2711         if (digits->HasData()) { //array
2712
2713           AliTRDSignalIndex   *indexes = digitsManager->GetIndexes(det);
2714           if (indexes->IsAllocated() == kFALSE) {
2715             AliError("Indexes do not exist!");
2716             break;
2717           }
2718           Int_t iRow  = 0;
2719           Int_t iCol  = 0;
2720           indexes->ResetCounters();
2721
2722           while (indexes->NextRCIndex(iRow, iCol)) { //column,row
2723             //printf(" det %d \t row %d \t col %d \t digit\n",det,iRow,iCol);
2724             //while (rawStream->Next()) {
2725
2726             Int_t idetector = det;                                             //  current detector
2727             //Int_t imcm      = rawStream->GetMCM();                            //  current MCM
2728             //Int_t irob      = rawStream->GetROB();                            //  current ROB
2729
2730
2731             if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)) {
2732               // Fill
2733               withInput = TMath::Max(FillDAQ(phvalue),withInput);
2734
2735               // reset
2736               for(Int_t k = 0; k < 36; k++){
2737                 for(Int_t j = 0; j < 16; j++){
2738                   for(Int_t c = 0; c < 144; c++){
2739                     phvalue[j][c][k] = 0.0;
2740                   }
2741                 }
2742               }
2743             }
2744
2745             fDetectorPreviousTrack = idetector;
2746             //fMCMPrevious           = imcm;
2747             //fROBPrevious           = irob;
2748
2749             //    nbtimebin              = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
2750             AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
2751             nbtimebin              = digitParam->GetNTimeBins(det);              //  number of time bins read from data
2752             baseline               = digitParam->GetADCbaseline(det);            //  baseline
2753
2754             if(nbtimebin == 0) return 0;
2755             if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2756             fTimeMax          = nbtimebin;
2757
2758             fNumberClustersf    = fTimeMax;
2759             fNumberClusters     = (Int_t)(fNumberClustersProcent*fTimeMax);
2760
2761
2762             for(Int_t itime = 0; itime < nbtimebin; itime++) {
2763               //            phvalue[row][col][itime] = signal[itime]-baseline;
2764               phvalue[iRow][iCol][itime] = (Short_t)(digits->GetData(iRow,iCol,itime) - baseline);
2765               /*if(phvalue[iRow][iCol][itime] >= 20) {
2766                  printf("----------> phvalue[%d][%d][%d] %d  baseline %d \n",
2767                        iRow,
2768                        iCol,
2769                        itime,
2770                        (Short_t)(digits->GetData(iRow,iCol,itime)),
2771                        baseline);
2772                        }*/
2773             }
2774
2775           }//column,row
2776
2777           // fill the last one
2778           if(fDetectorPreviousTrack != -1){
2779
2780             // Fill
2781             withInput = TMath::Max(FillDAQ(phvalue),withInput);
2782             //      printf("\n ---> withinput %d\n\n",withInput);
2783             // reset
2784             for(Int_t k = 0; k < 36; k++){
2785               for(Int_t j = 0; j < 16; j++){
2786                 for(Int_t c = 0; c < 144; c++){
2787                   phvalue[j][c][k] = 0.0;
2788                 }
2789               }
2790             }
2791           }
2792
2793         }//array
2794       }//QA
2795       digitsManager->ClearArrays(det);
2796     }//idetector
2797     delete digitsManager;
2798
2799     delete rawStream;
2800     return withInput;
2801  }//main
2802 //_____________________________________________________________________
2803 //////////////////////////////////////////////////////////////////////////////
2804 // Routine inside the DAQ process
2805 /////////////////////////////////////////////////////////////////////////////
2806 //_______________________________________________________________________
2807 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2808
2809   //
2810   // Look for the maximum by collapsing over the time
2811   // Sum over four pad col and two pad row
2812   //
2813
2814   Int_t used = 0;
2815
2816
2817   Int_t idect = fDetectorPreviousTrack;
2818   //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2819   Double_t sum[36];
2820   for(Int_t tb = 0; tb < 36; tb++){
2821     sum[tb] = 0.0;
2822   }
2823
2824   //fGoodTracklet = kTRUE;
2825   //fDetectorPreviousTrack = 0;
2826
2827
2828   ///////////////////////////
2829   // look for maximum
2830   /////////////////////////
2831
2832   Int_t imaxRow = 0;
2833   Int_t imaxCol = 0;
2834   Double_t integralMax = -1;
2835
2836   for (Int_t ir = 1; ir <= 15; ir++)
2837     {
2838       for (Int_t ic = 2; ic <= 142; ic++)
2839         {
2840           Double_t integral = 0;
2841           for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2842             {
2843               for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2844                 {
2845                   if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2846                       ic + ishiftC >= 1 && ic + ishiftC <= 144)
2847                     {
2848
2849                       for(Int_t tb = 0; tb< fTimeMax; tb++){
2850                         integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2851                       }// addtb
2852                     } //addsignal
2853                 } //shiftC
2854             } // shiftR
2855           if (integralMax < integral)
2856             {
2857               imaxRow = ir;
2858               imaxCol = ic;
2859               integralMax = integral;
2860
2861             } // check max integral
2862         } //ic
2863     } // ir
2864
2865   //  printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2866   //if((imaxRow == 0) || (imaxRow >= 15) || (imaxCol <= 3) || (imaxCol >= 140)) {
2867   //  used=1;
2868   //  return used;
2869   // }
2870
2871   if(((imaxRow + fNumberRowDAQ) > 16) || (imaxRow == 0) || ((imaxCol - fNumberColDAQ) <= 1) || ((imaxCol + fNumberColDAQ) >= 144)) {
2872     used=1;
2873     return used;
2874   }
2875   //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2876   //if(!fGoodTracklet) used = 1;;
2877
2878   //  /////////////////////////////////////////////////////
2879   // sum ober 2 row and 4 pad cols for each time bins
2880   //  ////////////////////////////////////////////////////
2881
2882
2883
2884   for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2885     {
2886       for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2887         {
2888           if (imaxRow + ishiftR >= 1 && imaxRow + ishiftR <= 16 &&
2889               imaxCol + ishiftC >= 1 && imaxCol + ishiftC <= 144)
2890             {
2891               for(Int_t it = 0; it < fTimeMax; it++){
2892                 sum[it] += phvalue[imaxRow + ishiftR-1][imaxCol + ishiftC-1][it];
2893               }
2894             }
2895         } // col shift
2896     }// row shift
2897
2898   Int_t nbcl = 0;
2899   Double_t sumcharge = 0.0;
2900   for(Int_t it = 0; it < fTimeMax; it++){
2901     sumcharge += sum[it];
2902     if(sum[it] > fThresholdClustersDAQ)  nbcl++;
2903   }
2904
2905
2906   /////////////////////////////////////////////////////////
2907   // Debug
2908   ////////////////////////////////////////////////////////
2909   if(fDebugLevel > 0){
2910     if ( !fDebugStreamer ) {
2911       //debug stream
2912       TDirectory *backup = gDirectory;
2913       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2914       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2915     }
2916
2917     Double_t amph0 = sum[0];
2918     Double_t amphlast = sum[fTimeMax-1];
2919     Double_t rms      = TMath::RMS(fTimeMax,sum);
2920     Int_t    goodtracklet = (Int_t) fGoodTracklet;
2921     for(Int_t it = 0; it < fTimeMax; it++){
2922       Double_t clustera = sum[it];
2923
2924     (* fDebugStreamer) << "FillDAQa"<<
2925       "ampTotal="<<sumcharge<<
2926       "row="<<imaxRow<<
2927       "col="<<imaxCol<<
2928       "detector="<<idect<<
2929       "amph0="<<amph0<<
2930       "amphlast="<<amphlast<<
2931       "goodtracklet="<<goodtracklet<<
2932       "clustera="<<clustera<<
2933       "it="<<it<<
2934       "rms="<<rms<<
2935       "nbcl="<<nbcl<<
2936       "\n";
2937     }
2938   }
2939
2940   ////////////////////////////////////////////////////////
2941   // fill
2942   ///////////////////////////////////////////////////////
2943   //printf("fNumberClusters %d, fNumberClustersf %d\n",fNumberClusters,fNumberClustersf);
2944   if(sum[0] > 100.0) used = 1;
2945   if(nbcl < fNumberClusters) used = 1;
2946   if(nbcl > fNumberClustersf) used = 1;
2947
2948   //if(fDetectorPreviousTrack == 15){
2949   //  printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2950   //}
2951   //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2952   if(used == 0){
2953     for(Int_t it = 0; it < fTimeMax; it++){
2954       if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2955       else{
2956         if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2957       }
2958       //if(fFillWithZero) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2959       //else{
2960       // if(sum[it] > 0.0) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2961       //}
2962     }
2963
2964
2965     //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2966     used = 2;
2967     //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2968
2969   }
2970
2971   return used;
2972
2973 }
2974 //____________Online trackling in AliTRDtrigger________________________________
2975 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2976 {
2977   //
2978   // For the DAQ
2979   // Fill a simple average pulse height
2980   //
2981
2982
2983   ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2984
2985
2986   return kTRUE;
2987
2988 }
2989 //____________Write_____________________________________________________
2990 //_____________________________________________________________________
2991 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2992 {
2993   //
2994   //  Write infos to file
2995   //
2996
2997   //For debugging
2998   if ( fDebugStreamer ) {
2999     delete fDebugStreamer;
3000     fDebugStreamer = 0x0;
3001   }
3002
3003   AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
3004                ,fNumberTrack
3005                ,fNumberUsedCh[0]
3006                ,fNumberUsedCh[1]
3007                ,fNumberUsedPh[0]
3008                ,fNumberUsedPh[1]));
3009
3010   TDirectory *backup = gDirectory;
3011   TString option;
3012
3013   if ( append )
3014     option = "update";
3015   else
3016     option = "recreate";
3017
3018   TFile f(filename,option.Data());
3019
3020   TStopwatch stopwatch;
3021   stopwatch.Start();
3022   if(fVector2d) {
3023     f.WriteTObject(fCalibraVector);
3024   }
3025
3026   if (fCH2dOn ) {
3027     if (fHisto2d) {
3028       f.WriteTObject(fCH2d);
3029     }
3030   }
3031   if (fPH2dOn ) {
3032     if (fHisto2d) {
3033       f.WriteTObject(fPH2d);
3034     }
3035   }
3036   if (fPRF2dOn) {
3037     if (fHisto2d) {
3038         f.WriteTObject(fPRF2d);
3039     }
3040   }
3041   if(fLinearFitterOn){
3042     if(fLinearFitterDebugOn) AnalyseLinearFitter();
3043     f.WriteTObject(fLinearVdriftFit);
3044   }
3045
3046   f.Close();
3047
3048   if ( backup ) backup->cd();
3049
3050   AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
3051                ,stopwatch.RealTime(),stopwatch.CpuTime()));
3052 }
3053 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
3054 // Stats stuff
3055 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
3056 //___________________________________________probe the histos__________________________________________________
3057 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
3058 {
3059   //
3060   // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
3061   // debug mode with 2 for TH2I and 3 for TProfile2D
3062   // It gives a pointer to a Double_t[7] with the info following...
3063   // [0] : number of calibration groups with entries
3064   // [1] : minimal number of entries found
3065   // [2] : calibration group number of the min
3066   // [3] : maximal number of entries found
3067   // [4] : calibration group number of the max
3068   // [5] : mean number of entries found
3069   // [6] : mean relative error
3070   //
3071
3072   Double_t *info = new Double_t[7];
3073
3074   // Number of Xbins (detectors or groups of pads)
3075   Int_t    nbins   = h->GetNbinsY(); //number of calibration groups
3076   Int_t    nxbins  = h->GetNbinsX(); //number of bins per histo
3077
3078   // Initialise
3079   Double_t nbwe = 0; //number of calibration groups with entries
3080   Double_t minentries = 0; //minimal number of entries found
3081   Double_t maxentries = 0; //maximal number of entries found
3082   Double_t placemin = 0; //calibration group number of the min
3083   Double_t placemax = -1; //calibration group number of the max
3084   Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
3085   Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
3086
3087   Double_t counter = 0;
3088
3089   //Debug
3090   TH1F *nbEntries = 0x0;//distribution of the number of entries
3091   TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
3092   TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
3093
3094   // Beginning of the loop over the calibration groups
3095   for (Int_t idect = 0; idect < nbins; idect++) {
3096
3097     TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
3098     projch->SetDirectory(0);
3099
3100     // Number of entries for this calibration group
3101     Double_t nentries = 0.0;
3102     if((i%2) == 0){
3103       for (Int_t k = 0; k < nxbins; k++) {
3104         nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
3105       }
3106     }
3107     else{
3108       for (Int_t k = 0; k < nxbins; k++) {
3109         nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
3110         if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
3111           meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
3112           counter++;
3113         }
3114       }
3115     }
3116
3117     //Debug
3118     if(i > 1){
3119       if(nentries > 0){
3120         if(!((Bool_t)nbEntries)) nbEntries = new TH1F("Number of entries","Number of entries",100,(Int_t)nentries/2,nentries*2);
3121         nbEntries->SetDirectory(0);
3122         nbEntries->Fill(nentries);
3123         if(!((Bool_t)nbEntriesPerGroup)) nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group",nbins,0,nbins);
3124         nbEntriesPerGroup->SetDirectory(0);
3125         nbEntriesPerGroup->Fill(idect+0.5,nentries);
3126         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));
3127         nbEntriesPerSp->SetDirectory(0);
3128         nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
3129       }
3130     }
3131
3132     //min amd max
3133     if(nentries > maxentries){
3134       maxentries = nentries;
3135       placemax = idect;
3136     }
3137     if(idect == 0) {
3138       minentries = nentries;
3139     }
3140     if(nentries < minentries){
3141       minentries = nentries;
3142       placemin = idect;
3143     }
3144     //nbwe
3145     if(nentries > 0) {
3146       nbwe++;
3147       meanstats += nentries;
3148     }
3149   }//calibration groups loop
3150
3151   if(nbwe > 0) meanstats /= nbwe;
3152   if(counter > 0) meanrelativerror /= counter;
3153
3154   AliInfo(Form("There are %f calibration groups with entries",nbwe));
3155   AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
3156   AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
3157   AliInfo(Form("The mean number of entries is %f",meanstats));
3158   if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
3159
3160   info[0] = nbwe;
3161   info[1] = minentries;
3162   info[2] = placemin;
3163   info[3] = maxentries;
3164   info[4] = placemax;
3165   info[5] = meanstats;
3166   info[6] = meanrelativerror;
3167
3168   if(nbEntries && nbEntriesPerSp && nbEntriesPerGroup){
3169     gStyle->SetPalette(1);
3170     gStyle->SetOptStat(1111);
3171     gStyle->SetPadBorderMode(0);
3172     gStyle->SetCanvasColor(10);
3173     gStyle->SetPadLeftMargin(0.13);
3174     gStyle->SetPadRightMargin(0.01);
3175     TCanvas *stat = new TCanvas("stat","",50,50,600,800);
3176     stat->Divide(2,1);
3177     stat->cd(1);
3178     nbEntries->Draw("");
3179     stat->cd(2);
3180     nbEntriesPerSp->SetStats(0);
3181     nbEntriesPerSp->Draw("");
3182     TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
3183     stat1->cd();
3184     nbEntriesPerGroup->SetStats(0);
3185     nbEntriesPerGroup->Draw("");
3186   }
3187
3188   return info;
3189
3190 }
3191 //____________________________________________________________________________
3192 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
3193 {
3194   //
3195   // Return a Int_t[4] with:
3196   // 0 Mean number of entries
3197   // 1 median of number of entries
3198   // 2 rms of number of entries
3199   // 3 number of group with entries
3200   //
3201
3202   Double_t *stat = new Double_t[4];
3203   stat[3]             = 0.0;
3204
3205   Int_t    nbofgroups = CalculateTotalNumberOfBins(0);
3206
3207   Double_t *weight = new Double_t[nbofgroups];
3208   Double_t *nonul = new Double_t[nbofgroups];
3209
3210   for(Int_t k = 0; k < nbofgroups; k++){
3211     if(fEntriesCH[k] > 0) {
3212       weight[k] = 1.0;
3213       nonul[(Int_t)stat[3]] = fEntriesCH[k];
3214       stat[3]++;
3215     }
3216     else weight[k] = 0.0;
3217   }
3218   stat[0]          = TMath::Mean(nbofgroups,fEntriesCH,weight);
3219   stat[1]          = TMath::Median(nbofgroups,fEntriesCH,weight);
3220   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
3221
3222   delete [] weight;
3223   delete [] nonul;
3224
3225   return stat;
3226
3227 }
3228 //____________________________________________________________________________
3229 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
3230 {
3231   //
3232   // Return a Int_t[4] with:
3233   // 0 Mean number of entries
3234   // 1 median of number of entries
3235   // 2 rms of number of entries
3236   // 3 number of group with entries
3237   //
3238
3239   Double_t *stat      = new Double_t[4];
3240   stat[3]             = 0.0;
3241
3242   Int_t    nbofgroups = 540;
3243   Double_t *weight    = new Double_t[nbofgroups];
3244   Int_t    *nonul     = new Int_t[nbofgroups];
3245
3246   for(Int_t k = 0; k < nbofgroups; k++){
3247     if(fEntriesLinearFitter[k] > 0) {
3248       weight[k] = 1.0;
3249       nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
3250       stat[3]++;
3251     }
3252     else weight[k] = 0.0;
3253   }
3254   stat[0]          = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
3255   stat[1]          = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
3256   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
3257
3258   delete [] weight;
3259   delete [] nonul;
3260
3261   return stat;
3262
3263 }
3264 //////////////////////////////////////////////////////////////////////////////////////
3265 // Create Histos
3266 //////////////////////////////////////////////////////////////////////////////////////
3267 //_____________________________________________________________________________
3268 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
3269 {
3270   //
3271   // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
3272   // If fNgroupprf is zero then no binning in tnp
3273   //
3274
3275   TString name("Nz");
3276   name += fCalibraMode->GetNz(2);
3277   name += "Nrphi";
3278   name += fCalibraMode->GetNrphi(2);
3279   name += "Ngp";
3280   name += fNgroupprf;
3281
3282   if(fNgroupprf != 0){
3283
3284     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3285                             ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
3286     fPRF2d->SetYTitle("Det/pad groups");
3287     fPRF2d->SetXTitle("Position x/W [pad width units]");
3288     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3289     fPRF2d->SetStats(0);
3290   }
3291   else{
3292     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
3293                             ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
3294     fPRF2d->SetYTitle("Det/pad groups");
3295     fPRF2d->SetXTitle("Position x/W [pad width units]");
3296     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
3297     fPRF2d->SetStats(0);
3298   }
3299
3300 }
3301
3302 //_____________________________________________________________________________
3303 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
3304 {
3305   //
3306   // Create the 2D histos
3307   //
3308
3309   TString name("Ver");
3310   name += fVersionVdriftUsed;
3311   name += "Subver";
3312   name += fSubVersionVdriftUsed;
3313   name += "Nz";
3314   name += fCalibraMode->GetNz(1);
3315   name += "Nrphi";
3316   name += fCalibraMode->GetNrphi(1);
3317
3318   fPH2d = new TProfile2D("PH2d",(const Char_t *) name
3319                          ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
3320                          ,nn,0,nn);
3321   fPH2d->SetYTitle("Det/pad groups");
3322   fPH2d->SetXTitle("time [#mus]");
3323   fPH2d->SetZTitle("<PH> [a.u.]");
3324   fPH2d->SetStats(0);
3325
3326 }
3327 //_____________________________________________________________________________
3328 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
3329 {
3330   //
3331   // Create the 2D histos
3332   //
3333
3334   TString name("Ver");
3335   name += fVersionGainUsed;
3336   name += "Subver";
3337   name += fSubVersionGainUsed;
3338   name += "Nz";
3339   name += fCalibraMode->GetNz(0);
3340   name += "Nrphi";
3341   name += fCalibraMode->GetNrphi(0);
3342
3343   fCH2d = new TH2I("CH2d",(const Char_t *) name
3344                    ,fNumberBinCharge,0,300,nn,0,nn);
3345   fCH2d->SetYTitle("Det/pad groups");
3346   fCH2d->SetXTitle("charge deposit [a.u]");
3347   fCH2d->SetZTitle("counts");
3348   fCH2d->SetStats(0);
3349   fCH2d->Sumw2();
3350
3351 }
3352 //////////////////////////////////////////////////////////////////////////////////
3353 // Set relative scale
3354 /////////////////////////////////////////////////////////////////////////////////
3355 //_____________________________________________________________________________
3356 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
3357 {
3358   //
3359   // Set the factor that will divide the deposited charge
3360   // to fit in the histo range [0,300]
3361   //
3362
3363   if (RelativeScale > 0.0) {
3364     fRelativeScale = RelativeScale;
3365   }
3366   else {
3367     AliInfo("RelativeScale must be strict positif!");
3368   }
3369
3370 }
3371 //////////////////////////////////////////////////////////////////////////////////
3372 // Quick way to fill a histo
3373 //////////////////////////////////////////////////////////////////////////////////
3374 //_____________________________________________________________________
3375 void  AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3376 {
3377   //
3378   // FillCH2d: Marian style
3379   //
3380
3381   //skip simply the value out of range
3382   if((y>=300.0) || (y<0.0)) return;
3383
3384   //Calcul the y place
3385   Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3386   Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3387
3388   //Fill
3389   fCH2d->GetArray()[place]++;
3390
3391 }
3392
3393 //////////////////////////////////////////////////////////////////////////////////
3394 // Geometrical functions
3395 ///////////////////////////////////////////////////////////////////////////////////
3396 //_____________________________________________________________________________
3397 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
3398 {
3399   //
3400   // Reconstruct the layer number from the detector number
3401   //
3402
3403   return ((Int_t) (d % 6));
3404
3405 }
3406
3407 //_____________________________________________________________________________
3408 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
3409 {
3410   //
3411   // Reconstruct the stack number from the detector number
3412   //
3413   const Int_t kNlayer = 6;
3414
3415   return ((Int_t) (d % 30) / kNlayer);
3416
3417 }
3418
3419 //_____________________________________________________________________________
3420 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3421 {
3422   //
3423   // Reconstruct the sector number from the detector number
3424   //
3425   Int_t fg = 30;
3426
3427   return ((Int_t) (d / fg));
3428
3429 }
3430 ///////////////////////////////////////////////////////////////////////////////////
3431 // Getter functions for DAQ of the CH2d and the PH2d
3432 //////////////////////////////////////////////////////////////////////////////////
3433 //_____________________________________________________________________
3434 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3435 {
3436     //
3437     // return pointer to fPH2d TProfile2D
3438     // create a new TProfile2D if it doesn't exist allready
3439     //
3440     if ( fPH2d )
3441         return fPH2d;
3442
3443     // Some parameters
3444     fTimeMax = nbtimebin;
3445     fSf      = samplefrequency;
3446
3447     CreatePH2d(540);
3448
3449     return fPH2d;
3450 }
3451 //_____________________________________________________________________
3452 TH2I* AliTRDCalibraFillHisto::GetCH2d()
3453 {
3454     //
3455     // return pointer to fCH2d TH2I
3456     // create a new TH2I if it doesn't exist allready
3457     //
3458     if ( fCH2d )
3459         return fCH2d;
3460
3461     CreateCH2d(540);
3462
3463     return fCH2d;
3464 }
3465 ////////////////////////////////////////////////////////////////////////////////////////////
3466 // Drift velocity calibration
3467 ///////////////////////////////////////////////////////////////////////////////////////////
3468 //_____________________________________________________________________
3469 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3470 {
3471     //
3472     // return pointer to TLinearFitter Calibration
3473     // if force is true create a new TLinearFitter if it doesn't exist allready
3474     //
3475
3476   if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3477     return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3478   }
3479
3480   // if we are forced and TLinearFitter doesn't yet exist create it
3481
3482   // new TLinearFitter
3483   TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3484   fLinearFitterArray.AddAt(linearfitter,detector);
3485   return linearfitter;
3486 }
3487
3488 //____________________________________________________________________________
3489 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3490 {
3491   //
3492   // Analyse array of linear fitter because can not be written
3493   // Store two arrays: one with the param the other one with the error param + number of entries
3494   //
3495
3496   for(Int_t k = 0; k < 540; k++){
3497     TLinearFitter *linearfitter = GetLinearFitter(k);
3498     if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3499       TVectorD  *par  = new TVectorD(2);
3500       TVectorD   pare = TVectorD(2);
3501       TVectorD  *parE = new TVectorD(3);
3502       linearfitter->Eval();
3503       linearfitter->GetParameters(*par);
3504       linearfitter->GetErrors(pare);
3505       Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3506       (*parE)[0] = pare[0]*ppointError;
3507       (*parE)[1] = pare[1]*ppointError;
3508       (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3509       ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3510       ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);
3511     }
3512   }
3513 }
3514
3515
3516