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   Float_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]+=((Float_t)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(TMath::Abs(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 = 0.0, rms = 0.0;
 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     if(signalSum > 0.0) mean/=signalSum;
 903
 904     for(Int_t ip = 0; ip<AliTRDfeeParam::GetNcol(); ip++)
 905       rms+=sigArr[ip]*(ip-mean)*(ip-mean);
 906
 907     if(signalSum > 0.0) rms = TMath::Sqrt(TMath::Abs(rms/signalSum));
 908
 909     rms -= TMath::Abs(dz*tilt);
 910     dydx -= dzdx*tilt;
 911   }
 912
 913   ////////////////////////////////
 914   // Debug stuff
 915   ///////////////////////////////
 916
 917
 918   //if(fDebugLevel > 0){
 919     if ( !fDebugStreamer ) {
 920       //debug stream
 921       TDirectory *backup = gDirectory;
 922       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
 923       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 924     }
 925
 926     float xcoord = tnp-dzdx*tnt;
 927     float pt = tracklet->GetPt();
 928     Int_t layer = GetLayer(fDetectorPreviousTrack);
 929
 930     (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
 931       //"snpright="<<snpright<<
 932       "nbli="<<nbli<<
 933       "nbclusters="<<nbclusters<<
 934       "detector="<<fDetectorPreviousTrack<<
 935       "layer="<<layer<<
 936       "snp="<<snp<<
 937       "tnp="<<tnp<<
 938       "tgl="<<tgl<<
 939       "tnt="<<tnt<<
 940       "dydt="<<dydt<<
 941       "dzdx="<<dzdx<<
 942       "crossrow="<<crossrow<<
 943       "errorpar="<<errorpar<<
 944       "pointError="<<pointError<<
 945       "xcoord="<<xcoord<<
 946       "pt="<<pt<<
 947       "rms="<<rms<<
 948       "dydx="<<dydx<<
 949       "dz="<<dz<<
 950       "tilt="<<tilt<<
 951       "ncl="<<ncl<<
 952       "\n";
 953
 954     //}
 955
 956   /////////////////////////
 957   // Cuts quality
 958   ////////////////////////
 959
 960   if(nbclusters < fNumberClusters) return kFALSE;
 961   if(nbclusters > fNumberClustersf) return kFALSE;
 962   if(pointError >= 0.3) return kFALSE;
 963   if(crossrow == 1) return kTRUE;
 964
 965   ///////////////////////
 966   // Fill
 967   //////////////////////
 968
 969   if(fLinearFitterOn){
 970     //Add to the linear fitter of the detector
 971     if( TMath::Abs(snp) <  1.){
 972       Double_t x = tnp-dzdx*tnt;
 973       if(fLinearFitterDebugOn) {
 974         (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
 975         fEntriesLinearFitter[fDetectorPreviousTrack]++;
 976       }
 977       fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
 978     }
 979   }
 980   if(fExbAltFitOn){
 981     fExbAltFit->Update(fDetectorPreviousTrack,dydx,rms);
 982   }
 983
 984   return kTRUE;
 985 }
 986 //____________Offine tracking in the AliTRDtracker_____________________________
 987 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
 988 {
 989   //
 990   // PRF width calibration
 991   // Assume a Gaussian shape: determinate the position of the three pad clusters
 992   // Fit with a straight line
 993   // Take the fitted values for all the clusters (3 or 2 pad clusters)
 994   // Fill the PRF as function of angle of the track
 995   //
 996   //
 997
 998   //printf("begin\n");
 999   ///////////////////////////////////////////
1000   // Take the parameters of the track
1001   //////////////////////////////////////////
1002   // take now the snp, tnp and tgl from the track
1003   Double_t snp = tracklet->GetSnp();             // sin dy/dx at the end of the chamber
1004   Double_t tnp = 0.0;                            // dy/dx at the end of the chamber
1005   if( TMath::Abs(snp) <  1.){
1006     tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1007   }
1008   Double_t tgl  = tracklet->GetTgl();           // dz/dl
1009   Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp);   // dz/dx calculated from dz/dl
1010   // at the entrance
1011   //Double_t tnp = tracklet->GetYref(1);      // dy/dx at the entrance of the chamber
1012   //Double_t tgl = tracklet->GetZref(1);      // dz/dl at the entrance of the chamber
1013   //Double_t dzdx = tgl;                      //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1014   // at the end with correction due to linear fit
1015   //Double_t tnp = tracklet->GetYfit(1);      // dy/dx at the end of the chamber after fit correction
1016   //Double_t tgl = tracklet->GetZfit(1);      // dz/dl at the end of the chamber after fit correction
1017
1018   ///////////////////////////////
1019   // Calculate tnp group shift
1020   ///////////////////////////////
1021   Bool_t echec   = kFALSE;
1022   Double_t shift = 0.0;
1023   //Calculate the shift in x coresponding to this tnp
1024   if(fNgroupprf != 0.0){
1025     shift      = -3.0*(fNgroupprf-1)-1.5;
1026     Double_t limithigh  = -0.2*(fNgroupprf-1);
1027     if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1028     else{
1029       while(tnp > limithigh){
1030         limithigh += 0.2;
1031         shift += 3.0;
1032       }
1033     }
1034   }
1035   // do nothing if out of tnp range
1036   //printf("echec %d\n",(Int_t)echec);
1037   if(echec) return kFALSE;
1038
1039   ///////////////////////
1040   // Variables
1041   //////////////////////
1042
1043   Int_t nb3pc    = 0;              // number of three pads clusters used for fit
1044   // to see the difference between the fit and the 3 pad clusters position
1045   Double_t padPositions[AliTRDseedV1::kNtb];
1046   memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
1047   fLinearFitterTracklet->ClearPoints();
1048
1049   //printf("loop clusters \n");
1050   ////////////////////////////
1051   // loop over the clusters
1052   ////////////////////////////
1053   AliTRDcluster *cl                   = 0x0;
1054   for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1055     // reject shared clusters on pad row
1056     if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
1057       cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1058       if(cl) continue;
1059     }
1060     cl = tracklet->GetClusters(ic);
1061     if(!cl) continue;
1062     Double_t     time  = cl->GetPadTime();
1063     if((time<=7) || (time>=21)) continue;
1064     Short_t  *signals  = cl->GetSignals();
1065     Float_t xcenter    = 0.0;
1066     Bool_t  echec1      = kTRUE;
1067
1068     /////////////////////////////////////////////////////////////
1069     // Center 3 balanced: position with the center of the pad
1070     /////////////////////////////////////////////////////////////
1071     if ((((Float_t) signals[3]) > 0.0) &&
1072         (((Float_t) signals[2]) > 0.0) &&
1073         (((Float_t) signals[4]) > 0.0)) {
1074       echec1 = kFALSE;
1075       // Security if the denomiateur is 0
1076       if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1077            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1078         xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1079           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1080                         / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1081       }
1082       else {
1083         echec1 = kTRUE;
1084       }
1085     }
1086     if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1087     if(echec1) continue;
1088
1089     ////////////////////////////////////////////////////////
1090     //if no echec1: calculate with the position of the pad
1091     // Position of the cluster
1092     // fill the linear fitter
1093     ///////////////////////////////////////////////////////
1094     Double_t       padPosition = xcenter +  cl->GetPadCol();
1095     padPositions[ic]            = padPosition;
1096     nb3pc++;
1097     fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1098
1099
1100   }//clusters loop
1101
1102   //printf("Fin loop clusters \n");
1103   //////////////////////////////
1104   // fit with a straight line
1105   /////////////////////////////
1106   if(nb3pc < 3){
1107     fLinearFitterTracklet->ClearPoints();
1108     return kFALSE;
1109   }
1110   fLinearFitterTracklet->Eval();
1111   TVectorD line(2);
1112   fLinearFitterTracklet->GetParameters(line);
1113   Float_t  pointError  = -1.0;
1114   if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1115   pointError  =  TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1116   }
1117   fLinearFitterTracklet->ClearPoints();
1118
1119   //printf("PRF second loop \n");
1120   ////////////////////////////////////////////////
1121   // Fill the PRF: Second loop over clusters
1122   //////////////////////////////////////////////
1123   for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1124     // reject shared clusters on pad row
1125     cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb);
1126     if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl)) continue;
1127     //
1128     cl = tracklet->GetClusters(ic);
1129     if(!cl) continue;
1130
1131     Short_t  *signals      = cl->GetSignals();              // signal
1132     Double_t     time      = cl->GetPadTime();         // time bin
1133     Float_t padPosTracklet = line[0]+line[1]*time;          // reconstruct position from fit
1134     Float_t padPos         = cl->GetPadCol();               // middle pad
1135     Double_t dpad          = padPosTracklet - padPos;       // reconstruct position relative to middle pad from fit
1136     Float_t ycenter        = 0.0;                           // relative center charge
1137     Float_t ymin           = 0.0;                           // relative left charge
1138     Float_t ymax           = 0.0;                           // relative right charge
1139
1140     ////////////////////////////////////////////////////////////////
1141     // Calculate ycenter, ymin and ymax even for two pad clusters
1142     ////////////////////////////////////////////////////////////////
1143     if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1144        ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1145       Float_t sum     = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1146       if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1147       if(sum > 0.0) ymin    = ((Float_t) signals[2])/ sum;
1148       if(sum > 0.0) ymax    = ((Float_t) signals[4])/ sum;
1149     }
1150
1151     /////////////////////////
1152     // Calibration group
1153     ////////////////////////
1154     Int_t     rowcl        = cl->GetPadRow();                           // row of cluster
1155     Int_t     colcl        = cl->GetPadCol();                           // col of cluster
1156     Int_t     grouplocal   = CalculateCalibrationGroup(2,rowcl,colcl);  // calcul the corresponding group
1157     Int_t     caligroup    = fCalibraMode->GetXbins(2)+ grouplocal;     // calcul the corresponding group
1158     Float_t   xcl          = cl->GetY();                                // y cluster
1159     Float_t   qcl          = cl->GetQ();                                // charge cluster
1160     Int_t     layer        = GetLayer(fDetectorPreviousTrack);          // layer
1161     Int_t     stack        = GetStack(fDetectorPreviousTrack);          // stack
1162     Double_t  xdiff        = dpad;                                      // reconstructed position constant
1163     Double_t  x            = dpad;                                      // reconstructed position moved
1164     Float_t   ep           = pointError;                                // error of fit
1165     Float_t   signal1      = (Float_t)signals[1];                       // signal at the border
1166     Float_t   signal3      = (Float_t)signals[3];                       // signal
1167     Float_t   signal2      = (Float_t)signals[2];                       // signal
1168     Float_t   signal4      = (Float_t)signals[4];                       // signal
1169     Float_t   signal5      = (Float_t)signals[5];                       // signal at the border
1170
1171
1172
1173     /////////////////////
1174     // Debug stuff
1175     ////////////////////
1176
1177     if(fDebugLevel > 0){
1178       if ( !fDebugStreamer ) {
1179         //debug stream
1180         TDirectory *backup = gDirectory;
1181         fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1182         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1183       }
1184
1185       x = xdiff;
1186       Int_t type=0;
1187       Float_t y = ycenter;
1188       (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1189         "caligroup="<<caligroup<<
1190         "detector="<<fDetectorPreviousTrack<<
1191         "layer="<<layer<<
1192         "stack="<<stack<<
1193         "npoints="<<nbclusters<<
1194         "Np="<<nb3pc<<
1195         "ep="<<ep<<
1196         "type="<<type<<
1197         "snp="<<snp<<
1198         "tnp="<<tnp<<
1199         "tgl="<<tgl<<
1200         "dzdx="<<dzdx<<
1201         "padPos="<<padPos<<
1202         "padPosition="<<padPositions[ic]<<
1203         "padPosTracklet="<<padPosTracklet<<
1204         "x="<<x<<
1205         "y="<<y<<
1206         "xcl="<<xcl<<
1207         "qcl="<<qcl<<
1208         "signal1="<<signal1<<
1209         "signal2="<<signal2<<
1210         "signal3="<<signal3<<
1211         "signal4="<<signal4<<
1212         "signal5="<<signal5<<
1213         "time="<<time<<
1214         "\n";
1215       x=-(xdiff+1);
1216       y = ymin;
1217       type=-1;
1218       (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1219         "caligroup="<<caligroup<<
1220         "detector="<<fDetectorPreviousTrack<<
1221         "layer="<<layer<<
1222         "stack="<<stack<<
1223         "npoints="<<nbclusters<<
1224         "Np="<<nb3pc<<
1225         "ep="<<ep<<
1226         "type="<<type<<
1227         "snp="<<snp<<
1228         "tnp="<<tnp<<
1229         "tgl="<<tgl<<
1230         "dzdx="<<dzdx<<
1231         "padPos="<<padPos<<
1232         "padPosition="<<padPositions[ic]<<
1233         "padPosTracklet="<<padPosTracklet<<
1234         "x="<<x<<
1235         "y="<<y<<
1236         "xcl="<<xcl<<
1237         "qcl="<<qcl<<
1238         "signal1="<<signal1<<
1239         "signal2="<<signal2<<
1240         "signal3="<<signal3<<
1241         "signal4="<<signal4<<
1242         "signal5="<<signal5<<
1243         "time="<<time<<
1244         "\n";
1245       x=1-xdiff;
1246       y = ymax;
1247       type=1;
1248       (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1249         "caligroup="<<caligroup<<
1250         "detector="<<fDetectorPreviousTrack<<
1251         "layer="<<layer<<
1252         "stack="<<stack<<
1253         "npoints="<<nbclusters<<
1254         "Np="<<nb3pc<<
1255         "ep="<<ep<<
1256         "type="<<type<<
1257         "snp="<<snp<<
1258         "tnp="<<tnp<<
1259         "tgl="<<tgl<<
1260         "dzdx="<<dzdx<<
1261         "padPos="<<padPos<<
1262         "padPosition="<<padPositions[ic]<<
1263         "padPosTracklet="<<padPosTracklet<<
1264         "x="<<x<<
1265         "y="<<y<<
1266         "xcl="<<xcl<<
1267         "qcl="<<qcl<<
1268         "signal1="<<signal1<<
1269         "signal2="<<signal2<<
1270         "signal3="<<signal3<<
1271         "signal4="<<signal4<<
1272         "signal5="<<signal5<<
1273         "time="<<time<<
1274         "\n";
1275
1276     }
1277
1278     /////////////////////
1279     // Cuts quality
1280     /////////////////////
1281     if(nbclusters < fNumberClusters) continue;
1282     if(nbclusters > fNumberClustersf) continue;
1283     if(nb3pc <= 5) continue;
1284     if((time >= 21) || (time < 7)) continue;
1285     if(TMath::Abs(qcl) < 80) continue;
1286     if( TMath::Abs(snp) >  1.) continue;
1287
1288
1289     ////////////////////////
1290     // Fill the histos
1291     ///////////////////////
1292     if (fHisto2d) {
1293       if(TMath::Abs(dpad) < 1.5) {
1294         fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1295         fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1296         //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1297       }
1298       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1299         fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1300         fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1301       }
1302       if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1303         fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1304         fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1305       }
1306     }
1307     // vector method
1308     if (fVector2d) {
1309       if(TMath::Abs(dpad) < 1.5) {
1310         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1311         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1312       }
1313       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1314         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1315         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1316       }
1317       if((ymax > 0.0)  && (TMath::Abs(dpad-1.0) < 1.5)) {
1318         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1319         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1320       }
1321     }
1322   } // second loop over clusters
1323
1324
1325   return kTRUE;
1326 }
1327 ///////////////////////////////////////////////////////////////////////////////////////
1328 // Pad row col stuff: see if masked or not
1329 ///////////////////////////////////////////////////////////////////////////////////////
1330 //_____________________________________________________________________________
1331 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(const AliTRDcluster *cl)
1332 {
1333   //
1334   // See if we are not near a masked pad
1335   //
1336
1337   if(cl->IsMasked()) fGoodTracklet = kFALSE;
1338
1339
1340 }
1341 //_____________________________________________________________________________
1342 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(const Int_t detector,const Int_t row,const Int_t col)
1343 {
1344   //
1345   // See if we are not near a masked pad
1346   //
1347
1348   if (!IsPadOn(detector, col, row)) {
1349     fGoodTracklet = kFALSE;
1350   }
1351
1352   if (col > 0) {
1353     if (!IsPadOn(detector, col-1, row)) {
1354       fGoodTracklet = kFALSE;
1355     }
1356   }
1357
1358   if (col < 143) {
1359     if (!IsPadOn(detector, col+1, row)) {
1360       fGoodTracklet = kFALSE;
1361     }
1362   }
1363
1364 }
1365 //_____________________________________________________________________________
1366 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1367 {
1368   //
1369   // Look in the choosen database if the pad is On.
1370   // If no the track will be "not good"
1371   //
1372
1373   // Get the parameter object
1374   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1375   if (!cal) {
1376     AliInfo("Could not get calibDB");
1377     return kFALSE;
1378   }
1379
1380   if (!cal->IsChamberGood(detector)     ||
1381        cal->IsChamberNoData(detector)        ||
1382        cal->IsPadMasked(detector,col,row)) {
1383     return kFALSE;
1384   }
1385   else {
1386     return kTRUE;
1387   }
1388
1389 }
1390 ///////////////////////////////////////////////////////////////////////////////////////
1391 // Calibration groups: calculate the number of groups, localise...
1392 ////////////////////////////////////////////////////////////////////////////////////////
1393 //_____________________________________________________________________________
1394 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1395 {
1396   //
1397   // Calculate the calibration group number for i
1398   //
1399
1400   // Row of the cluster and position in the pad groups
1401   Int_t posr = 0;
1402   if (fCalibraMode->GetNnZ(i) != 0) {
1403     posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1404   }
1405
1406
1407   // Col of the cluster and position in the pad groups
1408   Int_t posc = 0;
1409   if (fCalibraMode->GetNnRphi(i) != 0) {
1410     posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1411   }
1412
1413   return posc*fCalibraMode->GetNfragZ(i)+posr;
1414
1415 }
1416 //____________________________________________________________________________________
1417 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1418 {
1419   //
1420   // Calculate the total number of calibration groups
1421   //
1422
1423   Int_t ntotal = 0;
1424
1425   // All together
1426   if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1427     ntotal = 1;
1428     AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1429     return ntotal;
1430   }
1431
1432   // Per Supermodule
1433   if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1434     ntotal = 18;
1435     AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1436     return ntotal;
1437   }
1438
1439   // More
1440   fCalibraMode->ModePadCalibration(2,i);
1441   fCalibraMode->ModePadFragmentation(0,2,0,i);
1442   fCalibraMode->SetDetChamb2(i);
1443   ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1444   fCalibraMode->ModePadCalibration(0,i);
1445   fCalibraMode->ModePadFragmentation(0,0,0,i);
1446   fCalibraMode->SetDetChamb0(i);
1447   ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1448   AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1449   return ntotal;
1450
1451 }
1452 //_____________________________________________________________________________
1453 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1454 {
1455   //
1456   // Set the mode of calibration group in the z direction for the parameter i
1457   //
1458
1459   if ((Nz >= 0) &&
1460       (Nz <  5)) {
1461     fCalibraMode->SetNz(i, Nz);
1462   }
1463   else {
1464     AliInfo("You have to choose between 0 and 4");
1465   }
1466
1467 }
1468 //_____________________________________________________________________________
1469 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1470 {
1471   //
1472   // Set the mode of calibration group in the rphi direction for the parameter i
1473   //
1474
1475   if ((Nrphi >= 0) &&
1476       (Nrphi <  7)) {
1477     fCalibraMode->SetNrphi(i ,Nrphi);
1478   }
1479   else {
1480     AliInfo("You have to choose between 0 and 6");
1481   }
1482
1483 }
1484
1485 //_____________________________________________________________________________
1486 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1487 {
1488   //
1489   // Set the mode of calibration group all together
1490   //
1491   if(fVector2d == kTRUE) {
1492     AliInfo("Can not work with the vector method");
1493     return;
1494   }
1495   fCalibraMode->SetAllTogether(i);
1496
1497 }
1498
1499 //_____________________________________________________________________________
1500 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1501 {
1502   //
1503   // Set the mode of calibration group per supermodule
1504   //
1505   if(fVector2d == kTRUE) {
1506     AliInfo("Can not work with the vector method");
1507     return;
1508   }
1509   fCalibraMode->SetPerSuperModule(i);
1510
1511 }
1512
1513 //____________Set the pad calibration variables for the detector_______________
1514 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1515 {
1516   //
1517   // For the detector calcul the first Xbins and set the number of row
1518   // and col pads per calibration groups, the number of calibration
1519   // groups in the detector.
1520   //
1521
1522   // first Xbins of the detector
1523   if (fCH2dOn) {
1524     fCalibraMode->CalculXBins(detector,0);
1525   }
1526   if (fPH2dOn) {
1527     fCalibraMode->CalculXBins(detector,1);
1528   }
1529   if (fPRF2dOn) {
1530     fCalibraMode->CalculXBins(detector,2);
1531   }
1532
1533   // fragmentation of idect
1534   for (Int_t i = 0; i < 3; i++) {
1535     fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1536     fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1537                        , (Int_t) GetStack(detector)
1538                        , (Int_t) GetSector(detector),i);
1539   }
1540
1541   return kTRUE;
1542
1543 }
1544 //_____________________________________________________________________________
1545 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1546 {
1547   //
1548   // Should be between 0 and 6
1549   //
1550
1551   if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1552     AliInfo("The number of groups must be between 0 and 6!");
1553   }
1554   else {
1555     fNgroupprf = numberGroupsPRF;
1556   }
1557
1558 }
1559 ///////////////////////////////////////////////////////////////////////////////////////////
1560 // Per tracklet: store or reset the info, fill the histos with the info
1561 //////////////////////////////////////////////////////////////////////////////////////////
1562 //_____________________________________________________________________________
1563 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)
1564 {
1565   //
1566   // Store the infos in fAmpTotal, fPHPlace and fPHValue
1567   // Correct from the gain correction before
1568   // cls is shared cluster if any
1569   //
1570
1571   //printf("StoreInfoCHPHtrack\n");
1572
1573   // time bin of the cluster not corrected
1574   Int_t    time     = cl->GetPadTime();
1575   Float_t  charge   = TMath::Abs(cl->GetQ());
1576   if(cls) {
1577     charge += TMath::Abs(cls->GetQ());
1578     //printf("AliTRDCalibraFillHisto::Add the cluster charge");
1579   }
1580
1581   //printf("Store::time %d, amplitude %f\n",time,dqdl);
1582
1583   //Correct for the gain coefficient used in the database for reconstruction
1584   Float_t correctthegain = 1.0;
1585   if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1586   else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1587   Float_t correction    = 1.0;
1588   Float_t normalisation = 1.13; //org: 6.67; 1st: 1.056; 2nd: 1.13;
1589   // we divide with gain in AliTRDclusterizer::Transform...
1590   if( correctthegain > 0 ) normalisation /= correctthegain;
1591
1592
1593   // take dd/dl corrected from the angle of the track
1594   correction = dqdl / (normalisation);
1595
1596
1597   // Fill the fAmpTotal with the charge
1598   if (fCH2dOn) {
1599     if((!fLimitChargeIntegration) || (cl->IsInChamber())) {
1600       //printf("Store::group %d, amplitude %f\n",group[0],correction);
1601       fAmpTotal[(Int_t) group[0]] += correction;
1602     }
1603   }
1604
1605   // Fill the fPHPlace and value
1606   if (fPH2dOn) {
1607     fPHPlace[time] = group[1];
1608     fPHValue[time] = charge;
1609   }
1610
1611 }
1612 //____________Offine tracking in the AliTRDtracker_____________________________
1613 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1614 {
1615   //
1616   // Reset values per tracklet
1617   //
1618
1619   //Reset good tracklet
1620   fGoodTracklet = kTRUE;
1621
1622   // Reset the fPHValue
1623   if (fPH2dOn) {
1624     //Reset the fPHValue and fPHPlace
1625     for (Int_t k = 0; k < fTimeMax; k++) {
1626       fPHValue[k] = 0.0;
1627       fPHPlace[k] = -1;
1628     }
1629   }
1630
1631   // Reset the fAmpTotal where we put value
1632   if (fCH2dOn) {
1633     for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1634       fAmpTotal[k] = 0.0;
1635     }
1636   }
1637 }
1638 //____________Offine tracking in the AliTRDtracker_____________________________
1639 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1640 {
1641   //
1642   // For the offline tracking or mcm tracklets
1643   // This function will be called in the functions UpdateHistogram...
1644   // to fill the info of a track for the relativ gain calibration
1645   //
1646
1647   Int_t nb            =  0;   // Nombre de zones traversees
1648   Int_t fd            = -1;   // Premiere zone non nulle
1649   Float_t totalcharge = 0.0;  // Total charge for the supermodule histo
1650
1651   //printf("CH2d nbclusters %d, fNumberClusters %d, fNumberClustersf %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1652
1653   if(nbclusters < fNumberClusters) return;
1654   if(nbclusters > fNumberClustersf) return;
1655
1656
1657   // Normalize with the number of clusters
1658   Double_t normalizeCst = fRelativeScale;
1659   if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1660
1661   //printf("Number of groups in one detector %d\n",fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0));
1662
1663   // See if the track goes through different zones
1664   for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1665     //printf("fAmpTotal %f for %d\n",fAmpTotal[k],k);
1666     if (fAmpTotal[k] > 0.0) {
1667       totalcharge += fAmpTotal[k];
1668       nb++;
1669       if (nb == 1) {
1670         fd = k;
1671       }
1672     }
1673   }
1674
1675   //printf("CH2d: nb %d, fd %d, calibration group %d, amplitude %f, detector %d\n",nb,fd,fCalibraMode->GetXbins(0),fAmpTotal[fd]/normalizeCst,fDetectorPreviousTrack);
1676
1677   switch (nb)
1678     {
1679     case 1:
1680       fNumberUsedCh[0]++;
1681       fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1682       if (fHisto2d) {
1683         FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1684         //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1685       }
1686       if (fVector2d) {
1687         fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1688       }
1689       break;
1690     case 2:
1691       if ((fAmpTotal[fd]   > 0.0) &&
1692           (fAmpTotal[fd+1] > 0.0)) {
1693         // One of the two very big
1694         if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1695           if (fHisto2d) {
1696             FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1697             //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1698           }
1699           if (fVector2d) {
1700             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1701           }
1702           fNumberUsedCh[1]++;
1703           fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1704         }
1705         if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1706           if (fHisto2d) {
1707             FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
1708             //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
1709           }
1710           if (fVector2d) {
1711             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
1712           }
1713           fNumberUsedCh[1]++;
1714           fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1715         }
1716       }
1717       if (fCalibraMode->GetNfragZ(0) > 1) {
1718         if (fAmpTotal[fd] > 0.0) {
1719           if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1720             if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1721               // One of the two very big
1722               if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1723                 if (fHisto2d) {
1724                   FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
1725                   //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
1726                 }
1727                 if (fVector2d) {
1728                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
1729                 }
1730                 fNumberUsedCh[1]++;
1731                 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1732               }
1733               if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1734                 if (fHisto2d) {
1735                   FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1736                   //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1737                 }
1738                 fNumberUsedCh[1]++;
1739                 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1740                 if (fVector2d) {
1741                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
1742                 }
1743               }
1744             }
1745           }
1746         }
1747       }
1748       break;
1749     default: break;
1750   }
1751 }
1752 //____________Offine tracking in the AliTRDtracker_____________________________
1753 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1754 {
1755   //
1756   // For the offline tracking or mcm tracklets
1757   // This function will be called in the functions UpdateHistogram...
1758   // to fill the info of a track for the drift velocity  calibration
1759   //
1760
1761   Int_t nb  =  1; // Nombre de zones traversees 1, 2 ou plus de 3
1762   Int_t fd1 = -1; // Premiere zone non nulle
1763   Int_t fd2 = -1; // Deuxieme zone non nulle
1764   Int_t k1  = -1; // Debut de la premiere zone
1765   Int_t k2  = -1; // Debut de la seconde zone
1766   Int_t nbclusters = 0; // number of clusters
1767
1768
1769
1770   // See if the track goes through different zones
1771   for (Int_t k = 0; k < fTimeMax; k++) {
1772     if (fPHValue[k] > 0.0) {
1773       nbclusters++;
1774       if (fd1 == -1) {
1775         fd1 = fPHPlace[k];
1776         k1  = k;
1777       }
1778       if (fPHPlace[k] != fd1) {
1779         if (fd2 == -1) {
1780           k2  = k;
1781           fd2 = fPHPlace[k];
1782           nb  = 2;
1783         }
1784         if (fPHPlace[k] != fd2) {
1785           nb = 3;
1786         }
1787       }
1788     }
1789   }
1790
1791   // See if noise before and after
1792   if(fMaxCluster > 0) {
1793     if(fPHValue[0] > fMaxCluster) return;
1794     if(fTimeMax > fNbMaxCluster) {
1795       for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
1796         if(fPHValue[k] > fMaxCluster) return;
1797       }
1798     }
1799   }
1800
1801   //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
1802
1803   if(nbclusters < fNumberClusters) return;
1804   if(nbclusters > fNumberClustersf) return;
1805
1806   switch(nb)
1807     {
1808     case 1:
1809       fNumberUsedPh[0]++;
1810       for (Int_t i = 0; i < fTimeMax; i++) {
1811         if (fHisto2d) {
1812           if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1813           else {
1814             if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1815               }
1816           //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
1817         }
1818         if (fVector2d) {
1819           if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1820           else {
1821             if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1822           }
1823         }
1824       }
1825       break;
1826     case 2:
1827       if ((fd1 == fd2+1) ||
1828           (fd2 == fd1+1)) {
1829         // One of the two fast all the think
1830         if (k2 > (k1+fDifference)) {
1831           //we choose to fill the fd1 with all the values
1832           fNumberUsedPh[1]++;
1833           for (Int_t i = 0; i < fTimeMax; i++) {
1834             if (fHisto2d) {
1835               if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1836               else {
1837                 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1838                   }
1839             }
1840             if (fVector2d) {
1841               if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1842               else {
1843                 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1844                   }
1845             }
1846           }
1847         }
1848         if ((k2+fDifference) < fTimeMax) {
1849           //we choose to fill the fd2 with all the values
1850           fNumberUsedPh[1]++;
1851           for (Int_t i = 0; i < fTimeMax; i++) {
1852             if (fHisto2d) {
1853               if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1854               else {
1855                 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1856               }
1857             }
1858           if (fVector2d) {
1859             if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1860             else {
1861               if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1862             }
1863           }
1864           }
1865         }
1866       }
1867       // Two zones voisines sinon rien!
1868       if (fCalibraMode->GetNfragZ(1) > 1) {
1869         // Case 2
1870         if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1871           if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1872             // One of the two fast all the think
1873             if (k2 > (k1+fDifference)) {
1874               //we choose to fill the fd1 with all the values
1875               fNumberUsedPh[1]++;
1876               for (Int_t i = 0; i < fTimeMax; i++) {
1877                 if (fHisto2d) {
1878                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1879                   else {
1880                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1881                   }
1882                 }
1883                 if (fVector2d) {
1884                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1885                   else {
1886                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1887                   }
1888                 }
1889               }
1890             }
1891             if ((k2+fDifference) < fTimeMax) {
1892               //we choose to fill the fd2 with all the values
1893               fNumberUsedPh[1]++;
1894               for (Int_t i = 0; i < fTimeMax; i++) {
1895                 if (fHisto2d) {
1896                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1897                   else {
1898                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1899                   }
1900                 }
1901                 if (fVector2d) {
1902                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1903                   else {
1904                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1905                   }
1906                 }
1907               }
1908             }
1909           }
1910         }
1911         // Two zones voisines sinon rien!
1912         // Case 3
1913         if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
1914           if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
1915             // One of the two fast all the think
1916             if (k2 > (k1 + fDifference)) {
1917               //we choose to fill the fd1 with all the values
1918               fNumberUsedPh[1]++;
1919               for (Int_t i = 0; i < fTimeMax; i++) {
1920                 if (fHisto2d) {
1921                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1922                   else {
1923                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1924                   }
1925                 }
1926                 if (fVector2d) {
1927                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1928                   else {
1929                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1930                   }
1931                 }
1932               }
1933             }
1934             if ((k2+fDifference) < fTimeMax) {
1935               //we choose to fill the fd2 with all the values
1936               fNumberUsedPh[1]++;
1937               for (Int_t i = 0; i < fTimeMax; i++) {
1938                 if (fHisto2d) {
1939                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1940                   else {
1941                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1942                   }
1943                 }
1944                 if (fVector2d) {
1945                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1946                   else {
1947                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1948                   }
1949                 }
1950               }
1951             }
1952           }
1953         }
1954       }
1955       break;
1956     default: break;
1957     }
1958 }
1959 //////////////////////////////////////////////////////////////////////////////////////////
1960 // DAQ process functions
1961 /////////////////////////////////////////////////////////////////////////////////////////
1962 //_____________________________________________________________________
1963 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader)
1964  { //main
1965   //
1966   // Event Processing loop - AliTRDrawStream
1967   //
1968   // 0 timebin problem
1969   // 1 no input
1970   // 2 input
1971   // Same algorithm as TestBeam but different reader
1972   //
1973
1974   AliTRDrawStream *rawStream = new AliTRDrawStream(rawReader);
1975
1976   AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(kTRUE);
1977   digitsManager->CreateArrays();
1978
1979   Int_t withInput = 1;
1980
1981   Double_t phvalue[16][144][36];
1982   for(Int_t k = 0; k < 36; k++){
1983     for(Int_t j = 0; j < 16; j++){
1984       for(Int_t c = 0; c < 144; c++){
1985         phvalue[j][c][k] = 0.0;
1986       }
1987     }
1988   }
1989
1990   fDetectorPreviousTrack = -1;
1991   fMCMPrevious           = -1;
1992   fROBPrevious           = -1;
1993
1994   Int_t nbtimebin = 0;
1995   Int_t baseline  = 10;
1996
1997
1998     fTimeMax = 0;
1999
2000     Int_t det    = 0;
2001     while ((det = rawStream->NextChamber(digitsManager, NULL, NULL)) >= 0) { //idetector
2002
2003       if (digitsManager->GetIndexes(det)->HasEntry()) {//QA
2004         //      printf("there is ADC data on this chamber!\n");
2005
2006         AliTRDarrayADC *digits = (AliTRDarrayADC *) digitsManager->GetDigits(det); //mod
2007         if (digits->HasData()) { //array
2008
2009           AliTRDSignalIndex   *indexes = digitsManager->GetIndexes(det);
2010           if (indexes->IsAllocated() == kFALSE) {
2011             AliError("Indexes do not exist!");
2012             break;
2013           }
2014           Int_t iRow  = 0;
2015           Int_t iCol  = 0;
2016           indexes->ResetCounters();
2017
2018           while (indexes->NextRCIndex(iRow, iCol)) { //column,row
2019             //printf(" det %d \t row %d \t col %d \t digit\n",det,iRow,iCol);
2020             //while (rawStream->Next()) {
2021
2022             Int_t idetector = det;                                             //  current detector
2023             //Int_t imcm      = rawStream->GetMCM();                            //  current MCM
2024             //Int_t irob      = rawStream->GetROB();                            //  current ROB
2025
2026
2027             if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)) {
2028               // Fill
2029               withInput = TMath::Max(FillDAQ(phvalue),withInput);
2030
2031               // reset
2032               for(Int_t k = 0; k < 36; k++){
2033                 for(Int_t j = 0; j < 16; j++){
2034                   for(Int_t c = 0; c < 144; c++){
2035                     phvalue[j][c][k] = 0.0;
2036                   }
2037                 }
2038               }
2039             }
2040
2041             fDetectorPreviousTrack = idetector;
2042             //fMCMPrevious           = imcm;
2043             //fROBPrevious           = irob;
2044
2045             //    nbtimebin              = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
2046             AliTRDdigitsParam *digitParam = (AliTRDdigitsParam *)digitsManager->GetDigitsParam();
2047             nbtimebin              = digitParam->GetNTimeBins(det);              //  number of time bins read from data
2048             baseline               = digitParam->GetADCbaseline(det);            //  baseline
2049
2050             if(nbtimebin == 0) return 0;
2051             if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2052             fTimeMax          = nbtimebin;
2053
2054             fNumberClustersf    = fTimeMax;
2055             fNumberClusters     = (Int_t)(fNumberClustersProcent*fTimeMax);
2056
2057
2058             for(Int_t itime = 0; itime < nbtimebin; itime++) {
2059               //            phvalue[row][col][itime] = signal[itime]-baseline;
2060               phvalue[iRow][iCol][itime] = (Short_t)(digits->GetData(iRow,iCol,itime) - baseline);
2061               /*if(phvalue[iRow][iCol][itime] >= 20) {
2062                  printf("----------> phvalue[%d][%d][%d] %d  baseline %d \n",
2063                        iRow,
2064                        iCol,
2065                        itime,
2066                        (Short_t)(digits->GetData(iRow,iCol,itime)),
2067                        baseline);
2068                        }*/
2069             }
2070
2071           }//column,row
2072
2073           // fill the last one
2074           if(fDetectorPreviousTrack != -1){
2075
2076             // Fill
2077             withInput = TMath::Max(FillDAQ(phvalue),withInput);
2078             //      printf("\n ---> withinput %d\n\n",withInput);
2079             // reset
2080             for(Int_t k = 0; k < 36; k++){
2081               for(Int_t j = 0; j < 16; j++){
2082                 for(Int_t c = 0; c < 144; c++){
2083                   phvalue[j][c][k] = 0.0;
2084                 }
2085               }
2086             }
2087           }
2088
2089         }//array
2090       }//QA
2091       digitsManager->ClearArrays(det);
2092     }//idetector
2093     delete digitsManager;
2094
2095     delete rawStream;
2096     return withInput;
2097  }//main
2098 //_____________________________________________________________________
2099 //////////////////////////////////////////////////////////////////////////////
2100 // Routine inside the DAQ process
2101 /////////////////////////////////////////////////////////////////////////////
2102 //_______________________________________________________________________
2103 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2104
2105   //
2106   // Look for the maximum by collapsing over the time
2107   // Sum over four pad col and two pad row
2108   //
2109
2110   Int_t used = 0;
2111
2112
2113   Int_t idect = fDetectorPreviousTrack;
2114   //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2115   Double_t sum[36];
2116   for(Int_t tb = 0; tb < 36; tb++){
2117     sum[tb] = 0.0;
2118   }
2119
2120   //fGoodTracklet = kTRUE;
2121   //fDetectorPreviousTrack = 0;
2122
2123
2124   ///////////////////////////
2125   // look for maximum
2126   /////////////////////////
2127
2128   Int_t imaxRow = 0;
2129   Int_t imaxCol = 0;
2130   Double_t integralMax = -1;
2131
2132   for (Int_t ir = 1; ir <= 15; ir++)
2133     {
2134       for (Int_t ic = 2; ic <= 142; ic++)
2135         {
2136           Double_t integral = 0;
2137           for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2138             {
2139               for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2140                 {
2141                   if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2142                       ic + ishiftC >= 1 && ic + ishiftC <= 144)
2143                     {
2144
2145                       for(Int_t tb = 0; tb< fTimeMax; tb++){
2146                         integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2147                       }// addtb
2148                     } //addsignal
2149                 } //shiftC
2150             } // shiftR
2151           if (integralMax < integral)
2152             {
2153               imaxRow = ir;
2154               imaxCol = ic;
2155               integralMax = integral;
2156
2157             } // check max integral
2158         } //ic
2159     } // ir
2160
2161   //  printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2162   //if((imaxRow == 0) || (imaxRow >= 15) || (imaxCol <= 3) || (imaxCol >= 140)) {
2163   //  used=1;
2164   //  return used;
2165   // }
2166
2167   if(((imaxRow + fNumberRowDAQ) > 16) || (imaxRow == 0) || ((imaxCol - fNumberColDAQ) <= 1) || ((imaxCol + fNumberColDAQ) >= 144)) {
2168     used=1;
2169     return used;
2170   }
2171   //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2172   //if(!fGoodTracklet) used = 1;;
2173
2174   //  /////////////////////////////////////////////////////
2175   // sum ober 2 row and 4 pad cols for each time bins
2176   //  ////////////////////////////////////////////////////
2177
2178
2179
2180   for (Int_t ishiftR = 0; ishiftR < fNumberRowDAQ; ishiftR++)
2181     {
2182       for (Int_t ishiftC = -fNumberColDAQ; ishiftC < fNumberColDAQ; ishiftC++)
2183         {
2184           if (imaxRow + ishiftR >= 1 && imaxRow + ishiftR <= 16 &&
2185               imaxCol + ishiftC >= 1 && imaxCol + ishiftC <= 144)
2186             {
2187               for(Int_t it = 0; it < fTimeMax; it++){
2188                 sum[it] += phvalue[imaxRow + ishiftR-1][imaxCol + ishiftC-1][it];
2189               }
2190             }
2191         } // col shift
2192     }// row shift
2193
2194   Int_t nbcl = 0;
2195   Double_t sumcharge = 0.0;
2196   for(Int_t it = 0; it < fTimeMax; it++){
2197     sumcharge += sum[it];
2198     if(sum[it] > fThresholdClustersDAQ)  nbcl++;
2199   }
2200
2201
2202   /////////////////////////////////////////////////////////
2203   // Debug
2204   ////////////////////////////////////////////////////////
2205   if(fDebugLevel > 0){
2206     if ( !fDebugStreamer ) {
2207       //debug stream
2208       TDirectory *backup = gDirectory;
2209       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2210       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2211     }
2212
2213     Double_t amph0 = sum[0];
2214     Double_t amphlast = sum[fTimeMax-1];
2215     Double_t rms      = TMath::RMS(fTimeMax,sum);
2216     Int_t    goodtracklet = (Int_t) fGoodTracklet;
2217     for(Int_t it = 0; it < fTimeMax; it++){
2218       Double_t clustera = sum[it];
2219
2220     (* fDebugStreamer) << "FillDAQa"<<
2221       "ampTotal="<<sumcharge<<
2222       "row="<<imaxRow<<
2223       "col="<<imaxCol<<
2224       "detector="<<idect<<
2225       "amph0="<<amph0<<
2226       "amphlast="<<amphlast<<
2227       "goodtracklet="<<goodtracklet<<
2228       "clustera="<<clustera<<
2229       "it="<<it<<
2230       "rms="<<rms<<
2231       "nbcl="<<nbcl<<
2232       "\n";
2233     }
2234   }
2235
2236   ////////////////////////////////////////////////////////
2237   // fill
2238   ///////////////////////////////////////////////////////
2239   //printf("fNumberClusters %d, fNumberClustersf %d\n",fNumberClusters,fNumberClustersf);
2240   if(sum[0] > 100.0) used = 1;
2241   if(nbcl < fNumberClusters) used = 1;
2242   if(nbcl > fNumberClustersf) used = 1;
2243
2244   //if(fDetectorPreviousTrack == 15){
2245   //  printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2246   //}
2247   //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2248   if(used == 0){
2249     for(Int_t it = 0; it < fTimeMax; it++){
2250       if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2251       else{
2252         if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2253       }
2254       //if(fFillWithZero) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2255       //else{
2256       // if(sum[it] > 0.0) UpdateDAQ(0,0,0,it,sum[it],fTimeMax);
2257       //}
2258     }
2259
2260
2261     //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2262     used = 2;
2263     //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2264
2265   }
2266
2267   return used;
2268
2269 }
2270 //____________Online trackling in AliTRDtrigger________________________________
2271 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2272 {
2273   //
2274   // For the DAQ
2275   // Fill a simple average pulse height
2276   //
2277
2278
2279   ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2280
2281
2282   return kTRUE;
2283
2284 }
2285 //____________Write_____________________________________________________
2286 //_____________________________________________________________________
2287 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2288 {
2289   //
2290   //  Write infos to file
2291   //
2292
2293   //For debugging
2294   if ( fDebugStreamer ) {
2295     delete fDebugStreamer;
2296     fDebugStreamer = 0x0;
2297   }
2298
2299   AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2300                ,fNumberTrack
2301                ,fNumberUsedCh[0]
2302                ,fNumberUsedCh[1]
2303                ,fNumberUsedPh[0]
2304                ,fNumberUsedPh[1]));
2305
2306   TDirectory *backup = gDirectory;
2307   TString option;
2308
2309   if ( append )
2310     option = "update";
2311   else
2312     option = "recreate";
2313
2314   TFile f(filename,option.Data());
2315
2316   TStopwatch stopwatch;
2317   stopwatch.Start();
2318   if(fVector2d) {
2319     f.WriteTObject(fCalibraVector);
2320   }
2321
2322   if (fCH2dOn ) {
2323     if (fHisto2d) {
2324       f.WriteTObject(fCH2d);
2325     }
2326   }
2327   if (fPH2dOn ) {
2328     if (fHisto2d) {
2329       f.WriteTObject(fPH2d);
2330     }
2331   }
2332   if (fPRF2dOn) {
2333     if (fHisto2d) {
2334         f.WriteTObject(fPRF2d);
2335     }
2336   }
2337   if(fLinearFitterOn){
2338     if(fLinearFitterDebugOn) AnalyseLinearFitter();
2339     f.WriteTObject(fLinearVdriftFit);
2340   }
2341
2342   f.Close();
2343
2344   if ( backup ) backup->cd();
2345
2346   AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2347                ,stopwatch.RealTime(),stopwatch.CpuTime()));
2348 }
2349 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2350 // Stats stuff
2351 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2352 //___________________________________________probe the histos__________________________________________________
2353 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2354 {
2355   //
2356   // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2357   // debug mode with 2 for TH2I and 3 for TProfile2D
2358   // It gives a pointer to a Double_t[7] with the info following...
2359   // [0] : number of calibration groups with entries
2360   // [1] : minimal number of entries found
2361   // [2] : calibration group number of the min
2362   // [3] : maximal number of entries found
2363   // [4] : calibration group number of the max
2364   // [5] : mean number of entries found
2365   // [6] : mean relative error
2366   //
2367
2368   Double_t *info = new Double_t[7];
2369
2370   // Number of Xbins (detectors or groups of pads)
2371   Int_t    nbins   = h->GetNbinsY(); //number of calibration groups
2372   Int_t    nxbins  = h->GetNbinsX(); //number of bins per histo
2373
2374   // Initialise
2375   Double_t nbwe = 0; //number of calibration groups with entries
2376   Double_t minentries = 0; //minimal number of entries found
2377   Double_t maxentries = 0; //maximal number of entries found
2378   Double_t placemin = 0; //calibration group number of the min
2379   Double_t placemax = -1; //calibration group number of the max
2380   Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2381   Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2382
2383   Double_t counter = 0;
2384
2385   //Debug
2386   TH1F *nbEntries = 0x0;//distribution of the number of entries
2387   TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2388   TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2389
2390   // Beginning of the loop over the calibration groups
2391   for (Int_t idect = 0; idect < nbins; idect++) {
2392
2393     TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2394     projch->SetDirectory(0);
2395
2396     // Number of entries for this calibration group
2397     Double_t nentries = 0.0;
2398     if((i%2) == 0){
2399       for (Int_t k = 0; k < nxbins; k++) {
2400         nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2401       }
2402     }
2403     else{
2404       for (Int_t k = 0; k < nxbins; k++) {
2405         nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2406         if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2407           meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
2408           counter++;
2409         }
2410       }
2411     }
2412
2413     //Debug
2414     if(i > 1){
2415       if(nentries > 0){
2416         if(!((Bool_t)nbEntries)) nbEntries = new TH1F("Number of entries","Number of entries",100,(Int_t)nentries/2,nentries*2);
2417         nbEntries->SetDirectory(0);
2418         nbEntries->Fill(nentries);
2419         if(!((Bool_t)nbEntriesPerGroup)) nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group",nbins,0,nbins);
2420         nbEntriesPerGroup->SetDirectory(0);
2421         nbEntriesPerGroup->Fill(idect+0.5,nentries);
2422         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));
2423         nbEntriesPerSp->SetDirectory(0);
2424         nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
2425       }
2426     }
2427
2428     //min amd max
2429     if(nentries > maxentries){
2430       maxentries = nentries;
2431       placemax = idect;
2432     }
2433     if(idect == 0) {
2434       minentries = nentries;
2435     }
2436     if(nentries < minentries){
2437       minentries = nentries;
2438       placemin = idect;
2439     }
2440     //nbwe
2441     if(nentries > 0) {
2442       nbwe++;
2443       meanstats += nentries;
2444     }
2445   }//calibration groups loop
2446
2447   if(nbwe > 0) meanstats /= nbwe;
2448   if(counter > 0) meanrelativerror /= counter;
2449
2450   AliInfo(Form("There are %f calibration groups with entries",nbwe));
2451   AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
2452   AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
2453   AliInfo(Form("The mean number of entries is %f",meanstats));
2454   if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
2455
2456   info[0] = nbwe;
2457   info[1] = minentries;
2458   info[2] = placemin;
2459   info[3] = maxentries;
2460   info[4] = placemax;
2461   info[5] = meanstats;
2462   info[6] = meanrelativerror;
2463
2464   if(nbEntries && nbEntriesPerSp && nbEntriesPerGroup){
2465     gStyle->SetPalette(1);
2466     gStyle->SetOptStat(1111);
2467     gStyle->SetPadBorderMode(0);
2468     gStyle->SetCanvasColor(10);
2469     gStyle->SetPadLeftMargin(0.13);
2470     gStyle->SetPadRightMargin(0.01);
2471     TCanvas *stat = new TCanvas("stat","",50,50,600,800);
2472     stat->Divide(2,1);
2473     stat->cd(1);
2474     nbEntries->Draw("");
2475     stat->cd(2);
2476     nbEntriesPerSp->SetStats(0);
2477     nbEntriesPerSp->Draw("");
2478     TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
2479     stat1->cd();
2480     nbEntriesPerGroup->SetStats(0);
2481     nbEntriesPerGroup->Draw("");
2482   }
2483
2484   return info;
2485
2486 }
2487 //____________________________________________________________________________
2488 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
2489 {
2490   //
2491   // Return a Int_t[4] with:
2492   // 0 Mean number of entries
2493   // 1 median of number of entries
2494   // 2 rms of number of entries
2495   // 3 number of group with entries
2496   //
2497
2498   Double_t *stat = new Double_t[4];
2499   stat[3]             = 0.0;
2500
2501   Int_t    nbofgroups = CalculateTotalNumberOfBins(0);
2502
2503   Double_t *weight = new Double_t[nbofgroups];
2504   Double_t *nonul = new Double_t[nbofgroups];
2505
2506   for(Int_t k = 0; k < nbofgroups; k++){
2507     if(fEntriesCH[k] > 0) {
2508       weight[k] = 1.0;
2509       nonul[(Int_t)stat[3]] = fEntriesCH[k];
2510       stat[3]++;
2511     }
2512     else weight[k] = 0.0;
2513   }
2514   stat[0]          = TMath::Mean(nbofgroups,fEntriesCH,weight);
2515   stat[1]          = TMath::Median(nbofgroups,fEntriesCH,weight);
2516   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
2517
2518   delete [] weight;
2519   delete [] nonul;
2520
2521   return stat;
2522
2523 }
2524 //____________________________________________________________________________
2525 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
2526 {
2527   //
2528   // Return a Int_t[4] with:
2529   // 0 Mean number of entries
2530   // 1 median of number of entries
2531   // 2 rms of number of entries
2532   // 3 number of group with entries
2533   //
2534
2535   Double_t *stat      = new Double_t[4];
2536   stat[3]             = 0.0;
2537
2538   Int_t    nbofgroups = 540;
2539   Double_t *weight    = new Double_t[nbofgroups];
2540   Int_t    *nonul     = new Int_t[nbofgroups];
2541
2542   for(Int_t k = 0; k < nbofgroups; k++){
2543     if(fEntriesLinearFitter[k] > 0) {
2544       weight[k] = 1.0;
2545       nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
2546       stat[3]++;
2547     }
2548     else weight[k] = 0.0;
2549   }
2550   stat[0]          = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
2551   stat[1]          = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
2552   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
2553
2554   delete [] weight;
2555   delete [] nonul;
2556
2557   return stat;
2558
2559 }
2560 //////////////////////////////////////////////////////////////////////////////////////
2561 // Create Histos
2562 //////////////////////////////////////////////////////////////////////////////////////
2563 //_____________________________________________________________________________
2564 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
2565 {
2566   //
2567   // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
2568   // If fNgroupprf is zero then no binning in tnp
2569   //
2570
2571   TString name("Nz");
2572   name += fCalibraMode->GetNz(2);
2573   name += "Nrphi";
2574   name += fCalibraMode->GetNrphi(2);
2575   name += "Ngp";
2576   name += fNgroupprf;
2577
2578   if(fNgroupprf != 0){
2579
2580     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2581                             ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
2582     fPRF2d->SetYTitle("Det/pad groups");
2583     fPRF2d->SetXTitle("Position x/W [pad width units]");
2584     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2585     fPRF2d->SetStats(0);
2586   }
2587   else{
2588     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2589                             ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
2590     fPRF2d->SetYTitle("Det/pad groups");
2591     fPRF2d->SetXTitle("Position x/W [pad width units]");
2592     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2593     fPRF2d->SetStats(0);
2594   }
2595
2596 }
2597
2598 //_____________________________________________________________________________
2599 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
2600 {
2601   //
2602   // Create the 2D histos
2603   //
2604
2605   TString name("Ver");
2606   name += fVersionVdriftUsed;
2607   name += "Subver";
2608   name += fSubVersionVdriftUsed;
2609   name += "FirstRun";
2610   name += fFirstRunVdrift;
2611   name += "Nz";
2612   name += fCalibraMode->GetNz(1);
2613   name += "Nrphi";
2614   name += fCalibraMode->GetNrphi(1);
2615
2616   fPH2d = new TProfile2D("PH2d",(const Char_t *) name
2617                          ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
2618                          ,nn,0,nn);
2619   fPH2d->SetYTitle("Det/pad groups");
2620   fPH2d->SetXTitle("time [#mus]");
2621   fPH2d->SetZTitle("<PH> [a.u.]");
2622   fPH2d->SetStats(0);
2623
2624 }
2625 //_____________________________________________________________________________
2626 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
2627 {
2628   //
2629   // Create the 2D histos
2630   //
2631
2632   TString name("Ver");
2633   name += fVersionGainUsed;
2634   name += "Subver";
2635   name += fSubVersionGainUsed;
2636   name += "FirstRun";
2637   name += fFirstRunGain;
2638   name += "Nz";
2639   name += fCalibraMode->GetNz(0);
2640   name += "Nrphi";
2641   name += fCalibraMode->GetNrphi(0);
2642
2643   fCH2d = new TH2I("CH2d",(const Char_t *) name
2644                    ,fNumberBinCharge,0,300,nn,0,nn);
2645   fCH2d->SetYTitle("Det/pad groups");
2646   fCH2d->SetXTitle("charge deposit [a.u]");
2647   fCH2d->SetZTitle("counts");
2648   fCH2d->SetStats(0);
2649   fCH2d->Sumw2();
2650
2651 }
2652 //////////////////////////////////////////////////////////////////////////////////
2653 // Set relative scale
2654 /////////////////////////////////////////////////////////////////////////////////
2655 //_____________________________________________________________________________
2656 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
2657 {
2658   //
2659   // Set the factor that will divide the deposited charge
2660   // to fit in the histo range [0,300]
2661   //
2662
2663   if (RelativeScale > 0.0) {
2664     fRelativeScale = RelativeScale;
2665   }
2666   else {
2667     AliInfo("RelativeScale must be strict positif!");
2668   }
2669
2670 }
2671 //////////////////////////////////////////////////////////////////////////////////
2672 // Quick way to fill a histo
2673 //////////////////////////////////////////////////////////////////////////////////
2674 //_____________________________________________________________________
2675 void  AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
2676 {
2677   //
2678   // FillCH2d: Marian style
2679   //
2680
2681   //skip simply the value out of range
2682   if((y>=300.0) || (y<0.0)) return;
2683
2684   //Calcul the y place
2685   Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
2686   Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
2687
2688   //Fill
2689   fCH2d->GetArray()[place]++;
2690
2691 }
2692
2693 //////////////////////////////////////////////////////////////////////////////////
2694 // Geometrical functions
2695 ///////////////////////////////////////////////////////////////////////////////////
2696 //_____________________________________________________________________________
2697 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
2698 {
2699   //
2700   // Reconstruct the layer number from the detector number
2701   //
2702
2703   return ((Int_t) (d % 6));
2704
2705 }
2706
2707 //_____________________________________________________________________________
2708 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
2709 {
2710   //
2711   // Reconstruct the stack number from the detector number
2712   //
2713   const Int_t kNlayer = 6;
2714
2715   return ((Int_t) (d % 30) / kNlayer);
2716
2717 }
2718
2719 //_____________________________________________________________________________
2720 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2721 {
2722   //
2723   // Reconstruct the sector number from the detector number
2724   //
2725   Int_t fg = 30;
2726
2727   return ((Int_t) (d / fg));
2728
2729 }
2730 ///////////////////////////////////////////////////////////////////////////////////
2731 // Getter functions for DAQ of the CH2d and the PH2d
2732 //////////////////////////////////////////////////////////////////////////////////
2733 //_____________________________________________________________________
2734 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
2735 {
2736     //
2737     // return pointer to fPH2d TProfile2D
2738     // create a new TProfile2D if it doesn't exist allready
2739     //
2740     if ( fPH2d )
2741         return fPH2d;
2742
2743     // Some parameters
2744     fTimeMax = nbtimebin;
2745     fSf      = samplefrequency;
2746
2747     CreatePH2d(540);
2748
2749     return fPH2d;
2750 }
2751 //_____________________________________________________________________
2752 TH2I* AliTRDCalibraFillHisto::GetCH2d()
2753 {
2754     //
2755     // return pointer to fCH2d TH2I
2756     // create a new TH2I if it doesn't exist allready
2757     //
2758     if ( fCH2d )
2759         return fCH2d;
2760
2761     CreateCH2d(540);
2762
2763     return fCH2d;
2764 }
2765 ////////////////////////////////////////////////////////////////////////////////////////////
2766 // Drift velocity calibration
2767 ///////////////////////////////////////////////////////////////////////////////////////////
2768 //_____________________________________________________________________
2769 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
2770 {
2771     //
2772     // return pointer to TLinearFitter Calibration
2773     // if force is true create a new TLinearFitter if it doesn't exist allready
2774     //
2775
2776   if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
2777     return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
2778   }
2779
2780   // if we are forced and TLinearFitter doesn't yet exist create it
2781
2782   // new TLinearFitter
2783   TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
2784   fLinearFitterArray.AddAt(linearfitter,detector);
2785   return linearfitter;
2786 }
2787
2788 //____________________________________________________________________________
2789 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
2790 {
2791   //
2792   // Analyse array of linear fitter because can not be written
2793   // Store two arrays: one with the param the other one with the error param + number of entries
2794   //
2795
2796   for(Int_t k = 0; k < 540; k++){
2797     TLinearFitter *linearfitter = GetLinearFitter(k);
2798     if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
2799       TVectorD  *par  = new TVectorD(2);
2800       TVectorD   pare = TVectorD(2);
2801       TVectorD  *parE = new TVectorD(3);
2802       if((linearfitter->EvalRobust(0.8)==0)) {
2803         //linearfitter->Eval();
2804         linearfitter->GetParameters(*par);
2805         //linearfitter->GetErrors(pare);
2806         //Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
2807         //(*parE)[0] = pare[0]*ppointError;
2808         //(*parE)[1] = pare[1]*ppointError;
2809
2810         (*parE)[0] = 0.0;
2811         (*parE)[1] = 0.0;
2812         (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
2813         ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
2814         ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);
2815       }
2816     }
2817   }
2818 }
2819
2820
2821