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