TRD/AliTRDCalibraFillHisto.cxx

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