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