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