TRD/AliTRDCalibraFillHisto.cxx

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