TRD/AliTRDCalibraFillHisto.cxx

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