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