TRD/AliTRDCalibraFillHisto.cxx

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