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