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