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