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