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