TRD/AliTRDCalibraFillHisto.cxx

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