TRD/AliTRDCalibTask.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /////////////////////////////////////////////////////////////////////////////////
  18 //
  19 // AliTRDCalibTask
  20 //
  21 // Offline TRD calibration task
  22 //
  23 // Author:
  24 //   R. Bailhache (R.Bailhache@gsi.de)
  25 //
  26 //////////////////////////////////////////////////////////////////////////////////////
  27
  28 #include <iostream>
  29 using namespace std;
  30 #include "Riostream.h"
  31 #include "TChain.h"
  32 #include "TTree.h"
  33 #include "TFile.h"
  34 #include "TProfile2D.h"
  35 #include "TH2I.h"
  36 #include "TH1F.h"
  37 #include "TH2F.h"
  38 #include "TH2S.h"
  39 #include "TList.h"
  40 #include "TMath.h"
  41 #include "TObject.h"
  42 #include "TObjArray.h"
  43 #include "TString.h"
  44 #include "TCanvas.h"
  45 #include "TLegend.h"
  46 #include "TStyle.h"
  47 #include "TLine.h"
  48 #include "TIterator.h"
  49 #include "TLinearFitter.h"
  50 #include "TVectorD.h"
  51
  52 #include "AliAnalysisTask.h"
  53 #include "AliAnalysisManager.h"
  54
  55 #include "AliExternalTrackParam.h"
  56 #include "AliESDVertex.h"
  57 #include "AliESDEvent.h"
  58 #include "AliESDfriend.h"
  59 #include "AliCentrality.h"
  60 #include "AliESDInputHandler.h"
  61 #include "AliESDtrack.h"
  62 #include "AliESDfriendTrack.h"
  63 #include "AliTRDtrackV1.h"
  64 #include "AliTRDseedV1.h"
  65 #include "AliTRDcluster.h"
  66 #include "AliTRDgeometry.h"
  67 #include "AliESDtrackCuts.h"
  68 #include "AliESDVertex.h"
  69 #include "AliTRDCalDet.h"
  70
  71 #include "AliTRDCalibraVector.h"
  72 #include "AliTRDCalibraFillHisto.h"
  73 #include "AliTRDCalibraVdriftLinearFit.h"
  74 #include "AliTRDCalibraExbAltFit.h"
  75
  76 #include "AliTRDcalibDB.h"
  77 #include "AliCDBId.h"
  78 #include "AliLog.h"
  79
  80 #include "AliTRDCalibChamberStatus.h"
  81 #include "AliTRDdEdxBaseUtils.h"
  82 #include "AliTRDdEdxCalibHistArray.h"
  83 #include "AliTRDdEdxCalibUtils.h"
  84
  85 #include "AliTRDCalibTask.h"
  86
  87
  88 ClassImp(AliTRDCalibTask)
  89
  90 //________________________________________________________________________
  91 AliTRDCalibTask::AliTRDCalibTask(const char *name)
  92 : AliAnalysisTaskSE(name), fESD(0),
  93   fkEsdTrack(0),
  94   fFriendTrack(0),
  95   fCalibObject(0),
  96   fTrdTrack(0),
  97   fCl(0),
  98   fListHist(0),
  99   fTRDCalibraFillHisto(0),
 100   fTRDChamberStatus(0),
 101   fNEvents(0),
 102   fNEventsInput(0),
 103   fNEventsTrigger(0),
 104   fNbTRDTrack(0),
 105   fNbTRDTrackOffline(0),
 106   fNbTRDTrackStandalone(0),
 107   fNbTPCTRDtrack(0),
 108   fNbGoodTracks(0),
 109   fNbTimeBin(0),
 110   fNbTimeBinOffline(0),
 111   fNbTimeBinStandalone(0),
 112   fNbClusters(0),
 113   fNbClustersOffline(0),
 114   fNbClustersStandalone(0),
 115   fNbTracklets(0),
 116   fNbTrackletsOffline(0),
 117   fNbTrackletsStandalone(0),
 118   fAbsoluteGain(0),
 119   fTOFbc(0),
 120   fCH2dSum(0),
 121   fPH2dSum(0),
 122   fCH2dSM(0),
 123   fPH2dSM(0),
 124   fCH2dTest(0),
 125   fPH2dTest(0),
 126   fLinearVdriftTest(0),
 127   fOnInstance(kTRUE),
 128   fHisto2d(kTRUE),
 129   fVector2d(kFALSE),
 130   fVdriftLinear(kTRUE),
 131   fExbAlt(kFALSE),
 132   fDebugLevelTRDCalibraFillHisto(0),
 133   fNbTimeBins(0),
 134   fNumberBinCharge(100),
 135   fRangeCharge(150),
 136   fVdBindx(32),
 137   fVdBindy(70),
 138   fVdRangex(0.8),
 139   fVdRangey(1.4),
 140   fSelectTrigger(kTRUE),
 141   fSelectedTrigger(new TObjArray()),
 142   fRejected(kTRUE),
 143   fEsdTrackCuts(0),
 144   fRequirePrimaryVertex(kFALSE),
 145   fVtxTPC(kFALSE),
 146   fVtxSPD(kFALSE),
 147   fMinNbContributors(0),
 148   fRangePrimaryVertexZ(9999999.0),
 149   fRejectPileUpWithSPD(kFALSE),
 150   fRejectPileUpWithTOF(kFALSE),
 151   fRejectPileUpWithTOFOrITS(kFALSE),
 152   fMinNbTracks(9),
 153   fMaxNbTracks(999999999),
 154   fCutWithVdriftCalib(kFALSE),
 155   fMinNbTRDtracklets(0),
 156   fMinTRDMomentum(0),
 157   fScaleGainWithTPCSignal(kFALSE),
 158   fLow(0),
 159   fHigh(30),
 160   fFillZero(kFALSE),
 161   fNormalizeNbOfCluster(kFALSE),
 162   fRelativeScale(0.0),
 163   fMaxCluster(100.0),
 164   fNbMaxCluster(2),
 165   fOfflineTracks(kFALSE),
 166   fStandaloneTracks(kFALSE),
 167   fFirstRunGain(-1),
 168   fVersionGainUsed(-1),
 169   fSubVersionGainUsed(-1),
 170   fFirstRunGainLocal(-1),
 171   fVersionGainLocalUsed(-1),
 172   fSubVersionGainLocalUsed(-1),
 173   fFirstRunVdrift(-1),
 174   fVersionVdriftUsed(-1),
 175   fSubVersionVdriftUsed(-1),
 176   fFirstRunExB(-1),
 177   fVersionExBUsed(-1),
 178   fSubVersionExBUsed(-1),
 179   fCalDetGain(0x0),
 180   fMaxEvent(0),
 181   fCounter(0),
 182   fPHQon(kTRUE)
 183 {
 184   //
 185   // Default constructor
 186   //
 187
 188   fNz[0] = 0;
 189   fNz[1] = 0;
 190   fNz[2] = 0;
 191
 192   fNrphi[0] = 0;
 193   fNrphi[1] = 0;
 194   fNrphi[2] = 0;
 195
 196
 197   // Define input and output slots here
 198   // Input slot #0 works with a TChain
 199   DefineInput(0, TChain::Class());
 200
 201   // Output slot #0 writes into a TList container
 202   DefineOutput(1, TList::Class());
 203
 204
 205 }
 206 //____________________________________________________________________________________
 207 AliTRDCalibTask::~AliTRDCalibTask()
 208 {
 209   //
 210   // AliTRDCalibTask destructor
 211   //
 212
 213   // Pointeur
 214   if(fNEvents) delete fNEvents;
 215   if(fNEventsInput) delete fNEventsInput;
 216   if(fNEventsTrigger) delete fNEventsTrigger;
 217   if(fNbTRDTrack) delete fNbTRDTrack;
 218   if(fNbTRDTrackOffline) delete fNbTRDTrackOffline;
 219   if(fNbTRDTrackStandalone) delete fNbTRDTrackStandalone;
 220   if(fNbTPCTRDtrack) delete fNbTPCTRDtrack;
 221   if(fNbGoodTracks) delete fNbGoodTracks;
 222   if(fNbTimeBin) delete fNbTimeBin;
 223   if(fNbTimeBinOffline) delete fNbTimeBinOffline;
 224   if(fNbTimeBinStandalone) delete fNbTimeBinStandalone;
 225   if(fNbClusters) delete fNbClusters;
 226   if(fNbClustersOffline) delete fNbClustersOffline;
 227   if(fNbClustersStandalone) delete fNbClustersStandalone;
 228   if(fNbTracklets) delete fNbTracklets;
 229   if(fNbTrackletsOffline) delete fNbTrackletsOffline;
 230   if(fNbTrackletsStandalone) delete fNbTrackletsStandalone;
 231   if(fAbsoluteGain) delete fAbsoluteGain;
 232   if(fTOFbc) delete fTOFbc;
 233   if(fCH2dSum) delete fCH2dSum;
 234   if(fPH2dSum) delete fPH2dSum;
 235   if(fCH2dSM) delete fCH2dSM;
 236   if(fPH2dSM) delete fPH2dSM;
 237   if(fCH2dTest) delete fCH2dTest;
 238   if(fPH2dTest) delete fPH2dTest;
 239   if(fLinearVdriftTest) delete fLinearVdriftTest;
 240   if(IsPHQon()){
 241     AliTRDdEdxCalibUtils::DeleteHistArray();
 242   }
 243
 244   if(fCalDetGain) delete fCalDetGain;
 245
 246   if(fSelectedTrigger) {
 247     fSelectedTrigger->Delete();
 248     delete fSelectedTrigger;
 249   }
 250   if(fEsdTrackCuts) {
 251     delete fEsdTrackCuts;
 252   }
 253
 254   if(fTRDChamberStatus) delete fTRDChamberStatus;
 255
 256 }
 257
 258 //________________________________________________________________________
 259 void AliTRDCalibTask::UserCreateOutputObjects()
 260 {
 261   //
 262   // Create the histos
 263   //
 264   //cout << "AliTRDCalibTask::CreateOutputObjects() IN" << endl;
 265
 266   // Number of time bins
 267   if(fNbTimeBins==0) {
 268     AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 269     fNbTimeBins = cal->GetNumberOfTimeBinsDCS();
 270     if(fNbTimeBins <= 0){
 271       AliWarning(Form("No of TimeBins from DB [%d] use default [30]", fNbTimeBins));
 272       fNbTimeBins = 30;
 273     }
 274   }
 275
 276   // output list
 277   fListHist = new TList();
 278   fListHist->SetOwner();
 279
 280   // init chamber status
 281   fTRDChamberStatus = new AliTRDCalibChamberStatus();
 282   fTRDChamberStatus->Init();
 283   fListHist->Add(fTRDChamberStatus->GetSparseI());
 284
 285   // instance calibration
 286   fTRDCalibraFillHisto = AliTRDCalibraFillHisto::Instance();
 287   if(fOnInstance) {
 288     fTRDCalibraFillHisto->SetNumberBinCharge(fNumberBinCharge); // set number of bin of the charge distribution
 289     fTRDCalibraFillHisto->SetCutWithVdriftCalib(fCutWithVdriftCalib); // cut vdrift
 290     fTRDCalibraFillHisto->SetMinNbTRDtracklets(fMinNbTRDtracklets); // min number of TRD tracklets
 291     fTRDCalibraFillHisto->SetMinTRDMomentum(fMinTRDMomentum); // min TRD momentum
 292     fTRDCalibraFillHisto->SetHisto2d(fHisto2d); // choose to use histograms
 293     fTRDCalibraFillHisto->SetVector2d(fVector2d); // choose to use vectors
 294     fTRDCalibraFillHisto->SetCH2dOn();  // choose to calibrate the gain
 295     fTRDCalibraFillHisto->SetPH2dOn();  // choose to calibrate the drift velocity
 296     fTRDCalibraFillHisto->SetPRF2dOn(); // choose to look at the PRF
 297     fTRDCalibraFillHisto->SetLinearFitterOn(fVdriftLinear); // Other possibility vdrift VDRIFT
 298     fTRDCalibraFillHisto->SetLinearFitterDebugOn(fVdriftLinear); // Other possibility vdrift
 299     fTRDCalibraFillHisto->SetExbAltFitOn(fExbAlt); // Alternative method for exb
 300     fTRDCalibraFillHisto->SetScaleWithTPCSignal(fScaleGainWithTPCSignal); // Scale Gain with TPC signal
 301     for(Int_t k = 0; k < 3; k++){
 302       if(((fNz[k] != 10) && (fNrphi[k] != 10)) && ((fNz[k] != 100) && (fNrphi[k] != 100))) {
 303         fTRDCalibraFillHisto->SetNz(k,fNz[k]);                                    // Mode calibration
 304         fTRDCalibraFillHisto->SetNrphi(k,fNrphi[k]);                             // Mode calibration
 305       }
 306       else {
 307         if((fNz[k] == 100) && (fNrphi[k] == 100))  {
 308           if(fVector2d) AliInfo("The mode all together is not supported by the vector method");
 309           fTRDCalibraFillHisto->SetAllTogether(k);
 310         }
 311         if((fNz[k] == 10) && (fNrphi[k] == 10))  {
 312           if(fVector2d) AliInfo("The mode per supermodule is not supported by the vector method");
 313           fTRDCalibraFillHisto->SetPerSuperModule(k);
 314         }
 315       }
 316     }
 317     // Variables for how to fill
 318     fTRDCalibraFillHisto->SetFillWithZero(fFillZero);
 319     fTRDCalibraFillHisto->SetNormalizeNbOfCluster(fNormalizeNbOfCluster);
 320     fTRDCalibraFillHisto->SetMaxCluster(fMaxCluster);
 321     fTRDCalibraFillHisto->SetNbMaxCluster(fNbMaxCluster);
 322
 323     // Init with 30 timebins
 324     fTRDCalibraFillHisto->Init2Dhistos(fNbTimeBins); // initialise the histos
 325     fTRDCalibraFillHisto->SetNumberClusters(fLow); // At least 11 clusters
 326     fTRDCalibraFillHisto->SetNumberClustersf(fHigh); // At least 11 clusters
 327
 328     // For testing only
 329     fTRDCalibraFillHisto->SetDebugLevel(fDebugLevelTRDCalibraFillHisto); //debug stuff
 330
 331     if(fHisto2d) {
 332       fListHist->Add(fTRDCalibraFillHisto->GetCH2d());
 333       fListHist->Add(fTRDCalibraFillHisto->GetPH2d());
 334       fListHist->Add(fTRDCalibraFillHisto->GetPRF2d());
 335     }
 336     if(fVdriftLinear) {
 337       AliTRDCalibraVdriftLinearFit *fvdl = fTRDCalibraFillHisto->GetVdriftLinearFit();
 338       if(fvdl){
 339         fvdl->SetNbBindx(fVdBindx);
 340         fvdl->SetNbBindy(fVdBindy);
 341         fvdl->SetRangedx(fVdRangex);
 342         fvdl->SetRangedy(fVdRangey);
 343       }
 344       fListHist->Add((TObject *)fTRDCalibraFillHisto->GetVdriftLinearFit());
 345     }
 346     if(fVector2d) fListHist->Add((TObject *) fTRDCalibraFillHisto->GetCalibraVector()); //calibra vector
 347     if(fExbAlt) fListHist->Add((TObject *)fTRDCalibraFillHisto->GetExbAltFit());
 348   }
 349   fRelativeScale = fTRDCalibraFillHisto->GetRelativeScale(); // Get the relative scale for the gain
 350
 351   fNEvents = new TH1I(Form("NEvents_%s",(const char*)fName),"NEvents", 2, 0, 2);
 352   fListHist->Add(fNEvents);
 353   fNEventsTrigger = new TH1I(Form("NEventsTrigger_%s",(const char*)fName),"NEventsTrigger", 2, 0, 2);
 354   fListHist->Add(fNEventsTrigger);
 355   fNEventsInput = new TH1I(Form("NEventsInput_%s",(const char*)fName),"NEventsInput", 2, 0, 2);
 356   fListHist->Add(fNEventsInput);
 357
 358   // absolute gain calibration even without AliESDfriend
 359   Int_t nBinsPt = 25;
 360   Double_t minPt = 0.001;
 361   Double_t maxPt = 10.0;
 362
 363   Double_t *binLimLogPt = new Double_t[nBinsPt+1];
 364   Double_t *binLimPt    = new Double_t[nBinsPt+1];
 365   for(Int_t i=0; i<=nBinsPt; i++) binLimLogPt[i]=(Double_t)TMath::Log10(minPt) + (TMath::Log10(maxPt)-TMath::Log10(minPt))/nBinsPt*(Double_t)i ;
 366   for(Int_t i=0; i<=nBinsPt; i++) binLimPt[i]=(Double_t)TMath::Power(10,binLimLogPt[i]);
 367
 368   fAbsoluteGain = new TH2F(Form("AbsoluteGain_%s",(const char*)fName),"AbsoluteGain", 200, 0.0, 700.0, nBinsPt, binLimPt);
 369   fAbsoluteGain->SetYTitle("Momentum at TRD");
 370   fAbsoluteGain->SetXTitle("charge deposit [a.u]");
 371   fAbsoluteGain->SetZTitle("counts");
 372   fAbsoluteGain->SetStats(0);
 373   fAbsoluteGain->Sumw2();
 374   fListHist->Add(fAbsoluteGain);
 375
 376   fTOFbc = new TH2F(Form("TOFbc_%s",(const char*)fName),"TOFbc", 200, -150., 50., nBinsPt, binLimPt);
 377   fTOFbc->SetYTitle("Momentum");
 378   fTOFbc->SetXTitle("TOF bc");
 379   fTOFbc->SetZTitle("counts");
 380   fTOFbc->SetStats(0);
 381   fTOFbc->Sumw2();
 382   fListHist->Add(fTOFbc);
 383
 384   if(IsPHQon()){
 385     printf("\n        AliTRDCalibTask PHQ is on!!     \n\n");
 386     AliTRDdEdxBaseUtils::PrintControl();
 387     AliTRDdEdxCalibUtils::IniHistArray(fListHist, kTRUE);
 388   }
 389   else{
 390     printf("\n        AliTRDCalibTask PHQ is off!!     \n\n");
 391   }
 392
 393   /////////////////////////////////////////
 394   // First debug level
 395   ///////////////////////////////////////
 396   if(fDebug > 0) {
 397
 398     fLinearVdriftTest = new TH2S(Form("LFDV0testversion_%s",(const char*)fName),"LFDV0testversion",36,-0.9,0.9,48,-1.2,1.2);
 399     fLinearVdriftTest->SetXTitle("tan(phi_{track})");
 400     fLinearVdriftTest->SetYTitle("dy/dt");
 401     fLinearVdriftTest->SetZTitle("Number of tracklets");
 402     fLinearVdriftTest->SetStats(0);
 403     fLinearVdriftTest->SetDirectory(0);
 404
 405     // Standart with AliESDfriend
 406     fPH2dTest = new TProfile2D(Form("PH2dTest_%s",(const char*)fName),"Nz0Nrphi0"
 407                             ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)
 408                            ,540,0,540);
 409     fPH2dTest->SetYTitle("Det/pad groups");
 410     fPH2dTest->SetXTitle("time [#mus]");
 411     fPH2dTest->SetZTitle("<PH> [a.u.]");
 412     fPH2dTest->SetStats(0);
 413     //
 414     fCH2dTest = new TH2I(Form("CH2dTest_%s",(const char*)fName),"Nz0Nrphi0",50,0,300,540,0,540);
 415     fCH2dTest->SetYTitle("Det/pad groups");
 416     fCH2dTest->SetXTitle("charge deposit [a.u]");
 417     fCH2dTest->SetZTitle("counts");
 418     fCH2dTest->SetStats(0);
 419     fCH2dTest->Sumw2();
 420
 421     //
 422     fPH2dSM = new TProfile2D(Form("PH2dSM_%s",(const char*)fName),"Nz10Nrphi10"
 423                             ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)
 424                            ,18,0,18);
 425     fPH2dSM->SetYTitle("Det/pad groups");
 426     fPH2dSM->SetXTitle("time [#mus]");
 427     fPH2dSM->SetZTitle("<PH> [a.u.]");
 428     fPH2dSM->SetStats(0);
 429     //
 430     fCH2dSM = new TH2I(Form("CH2dSM_%s",(const char*)fName),"Nz10Nrphi10",50,0,300,18,0,18);
 431     fCH2dSM->SetYTitle("Det/pad groups");
 432     fCH2dSM->SetXTitle("charge deposit [a.u]");
 433     fCH2dSM->SetZTitle("counts");
 434     fCH2dSM->SetStats(0);
 435     fCH2dSM->Sumw2();
 436     //
 437     fPH2dSum = new TProfile2D(Form("PH2dSum_%s",(const char*)fName),"Nz100Nrphi100"
 438                             ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)
 439                             ,1,0,1);
 440     fPH2dSum->SetYTitle("Det/pad groups");
 441     fPH2dSum->SetXTitle("time [#mus]");
 442     fPH2dSum->SetZTitle("<PH> [a.u.]");
 443     fPH2dSum->SetStats(0);
 444     //
 445     fCH2dSum = new TH2I(Form("CH2dSum_%s",(const char*)fName),"Nz100Nrphi100",50,0,300,1,0,1);
 446     fCH2dSum->SetYTitle("Det/pad groups");
 447     fCH2dSum->SetXTitle("charge deposit [a.u]");
 448     fCH2dSum->SetZTitle("counts");
 449     fCH2dSum->SetStats(0);
 450     fCH2dSum->Sumw2();
 451
 452
 453     // Add them
 454     fListHist->Add(fLinearVdriftTest);
 455     fListHist->Add(fPH2dTest);
 456     fListHist->Add(fCH2dTest);
 457     fListHist->Add(fPH2dSM);
 458     fListHist->Add(fCH2dSM);
 459     fListHist->Add(fPH2dSum);
 460     fListHist->Add(fCH2dSum);
 461
 462   }
 463
 464   /////////////////////////////////////////
 465   // Second debug level
 466   ///////////////////////////////////////
 467   if(fDebug > 1) {
 468
 469     fNbGoodTracks = new TH2F(Form("NbGoodTracks_%s",(const char*)fName),"NbGoodTracks",500,0.0,2500.0,200,0.0,100.0);
 470     fNbGoodTracks->SetXTitle("Nb of good tracks");
 471     fNbGoodTracks->SetYTitle("Centrality");
 472     fNbGoodTracks->SetStats(0);
 473
 474     fNbTRDTrack = new TH1F(Form("TRDTrack_%s",(const char*)fName),"TRDTrack",50,0,50);
 475     fNbTRDTrack->Sumw2();
 476     fNbTRDTrackOffline = new TH1F(Form("TRDTrackOffline_%s",(const char*)fName),"TRDTrackOffline",50,0,50);
 477     fNbTRDTrackOffline->Sumw2();
 478     fNbTRDTrackStandalone = new TH1F(Form("TRDTrackStandalone_%s",(const char*)fName),"TRDTrackStandalone",50,0,50);
 479     fNbTRDTrackStandalone->Sumw2();
 480     fNbTPCTRDtrack = new TH2F(Form("NbTPCTRDtrack_%s",(const char*)fName),"NbTPCTRDtrack",100,0,100,100,0,100);
 481     fNbTPCTRDtrack->Sumw2();
 482     //
 483     fNbTimeBin = new TH1F(Form("NbTimeBin_%s",(const char*)fName),"NbTimeBin",35,0,35);
 484     fNbTimeBin->Sumw2();
 485     fNbTimeBinOffline = new TH1F(Form("NbTimeBinOffline_%s",(const char*)fName),"NbTimeBinOffline",35,0,35);
 486     fNbTimeBinOffline->Sumw2();
 487     fNbTimeBinStandalone = new TH1F(Form("NbTimeBinStandalone_%s",(const char*)fName),"NbTimeBinStandalone",35,0,35);
 488     fNbTimeBinStandalone->Sumw2();
 489     //
 490     fNbClusters = new TH1F(Form("NbClusters_%s",(const char*)fName),"",35,0,35);
 491     fNbClusters->Sumw2();
 492     fNbClustersOffline = new TH1F(Form("NbClustersOffline_%s",(const char*)fName),"",35,0,35);
 493     fNbClustersOffline->Sumw2();
 494     fNbClustersStandalone = new TH1F(Form("NbClustersStandalone_%s",(const char*)fName),"",35,0,35);
 495     fNbClustersStandalone->Sumw2();
 496     //
 497     fNbTracklets = new TH1F(Form("NbTracklets_%s",(const char*)fName),"NbTracklets",540,0.,540.);
 498     fNbTracklets->Sumw2();
 499     fNbTrackletsOffline = new TH1F(Form("NbTrackletsOffline_%s",(const char*)fName),"NbTrackletsOffline",540,0.,540.);
 500     fNbTrackletsOffline->Sumw2();
 501     fNbTrackletsStandalone = new TH1F(Form("NbTrackletsStandalone_%s",(const char*)fName),"NbTrackletsStandalone",540,0.,540.);
 502     fNbTrackletsStandalone->Sumw2();
 503
 504     fListHist->Add(fNbGoodTracks);
 505
 506     fListHist->Add(fNbTRDTrack);
 507     fListHist->Add(fNbTRDTrackOffline);
 508     fListHist->Add(fNbTRDTrackStandalone);
 509     fListHist->Add(fNbTPCTRDtrack);
 510
 511     fListHist->Add(fNbTimeBin);
 512     fListHist->Add(fNbTimeBinOffline);
 513     fListHist->Add(fNbTimeBinStandalone);
 514     fListHist->Add(fNbClusters);
 515     fListHist->Add(fNbClustersOffline);
 516     fListHist->Add(fNbClustersStandalone);
 517     fListHist->Add(fNbTracklets);
 518     fListHist->Add(fNbTrackletsOffline);
 519     fListHist->Add(fNbTrackletsStandalone);
 520
 521   }
 522
 523   delete [] binLimLogPt;
 524   delete [] binLimPt;
 525
 526   PostData(1,fListHist);
 527
 528   //cout << "AliTRDCalibTask::UserCreateOutputObjects() OUT" << endl;
 529
 530 }
 531
 532 //________________________________________________________________________
 533 void AliTRDCalibTask::UserExec(Option_t *)
 534 {
 535   //
 536   // Filling of the histos
 537   //
 538   //cout << "AliTRDCalibTask::Exec() IN" << endl;
 539
 540   // Init Versions and subversions used
 541   if((fFirstRunGain==-1) || (fVersionGainUsed==-1) || (fSubVersionGainUsed==-1) || (fFirstRunGainLocal==-1) || (fVersionGainLocalUsed==-1) || (fSubVersionGainLocalUsed==-1) || (fFirstRunVdrift==-1) || (fVersionVdriftUsed==-1) || (fSubVersionVdriftUsed==-1)) {
 542     if(!SetVersionSubversion()) {
 543       PostData(1, fListHist);
 544       return;
 545     }
 546   }
 547   if(fCounter==0) {
 548     if(fOnInstance) {
 549       fTRDCalibraFillHisto->SetFirstRunGain(fFirstRunGain); // Gain Used
 550       fTRDCalibraFillHisto->SetVersionGainUsed(fVersionGainUsed); // Gain Used
 551       fTRDCalibraFillHisto->SetSubVersionGainUsed(fSubVersionGainUsed); // Gain Used
 552       fTRDCalibraFillHisto->SetFirstRunGainLocal(fFirstRunGainLocal); // Gain Used
 553       fTRDCalibraFillHisto->SetVersionGainLocalUsed(fVersionGainLocalUsed); // Gain Used
 554       fTRDCalibraFillHisto->SetSubVersionGainLocalUsed(fSubVersionGainLocalUsed); // Gain Used
 555       fTRDCalibraFillHisto->SetFirstRunVdrift(fFirstRunVdrift); // Vdrift Used
 556       fTRDCalibraFillHisto->SetVersionVdriftUsed(fVersionVdriftUsed); // Vdrift Used
 557       fTRDCalibraFillHisto->SetSubVersionVdriftUsed(fSubVersionVdriftUsed); // Vdrift Used
 558       if((fFirstRunExB != -1) && (fVersionExBUsed != -1) && (fSubVersionExBUsed != -1)){
 559         fTRDCalibraFillHisto->SetFirstRunExB(fFirstRunExB); // ExB Used
 560         fTRDCalibraFillHisto->SetVersionExBUsed(fVersionExBUsed); // ExB Used
 561         fTRDCalibraFillHisto->SetSubVersionExBUsed(fSubVersionExBUsed); // ExB Used
 562       }
 563       fTRDCalibraFillHisto->InitCalDet();
 564     }
 565     if(fDebug > 1){
 566       // title CH2dTest
 567       TString name("Ver");
 568       name += fVersionGainUsed;
 569       name += "Subver";
 570       name += fSubVersionGainUsed;
 571       name += "FirstRun";
 572       name += fFirstRunGain;
 573       name += "Nz0Nrphi0";
 574       fCH2dTest->SetTitle(name);
 575       // title PH2dTest
 576       TString namee("Ver");
 577       namee += fVersionVdriftUsed;
 578       namee += "Subver";
 579       namee += fSubVersionVdriftUsed;
 580       namee += "FirstRun";
 581       namee += fFirstRunVdrift;
 582       namee += "Nz0Nrphi0";
 583       fPH2dTest->SetTitle(namee);
 584     }
 585   }
 586
 587   //  AliLog::SetGlobalLogLevel(AliLog::kError);
 588   //  cout << "AliTRDCalibTask::Exec() 1" << endl;
 589   fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
 590   if(!fESD){
 591     AliError("ESD Event missing");
 592     PostData(1, fListHist);
 593     return;
 594   }
 595
 596   const char* type = fESD->GetBeamType();
 597
 598
 599   //printf("Counter %d\n",fCounter);
 600
 601   fCounter++;
 602   fNEventsInput->Fill(1);
 603
 604   //cout << "maxEvent = " << fMaxEvent << endl;
 605   //if(fCounter%100==0) cout << "fCounter = " << fCounter << endl;
 606   if((fMaxEvent != 0) && (fMaxEvent < fCounter)) {
 607     PostData(1, fListHist);
 608     return;
 609   }
 610   //if(fCounter%100==0) cout << "fCounter1 = " << fCounter << endl;
 611   //cout << "event = " << fCounter << endl;
 612
 613   //printf("Counter %d\n",fCounter);
 614
 615   ///////////////////
 616   // Check trigger
 617   ///////////////////
 618   Bool_t pass = kTRUE;
 619
 620   if (fSelectTrigger) {
 621
 622     //printf("Will check the triggers\n");
 623
 624     Int_t numberOfTriggerSelected = fSelectedTrigger->GetEntriesFast();
 625     //printf("numberofTriggerSelected %d\n",numberOfTriggerSelected);
 626     if(fRejected) {
 627       pass = kTRUE;
 628       for(Int_t k = 0; k < numberOfTriggerSelected; k++){
 629         const TObjString *const obString=(TObjString*)fSelectedTrigger->At(k);
 630         const TString tString=obString->GetString();
 631         if(fESD->IsTriggerClassFired((const char*)tString)) {
 632           pass = kFALSE;
 633         }
 634       }
 635     }
 636     else {
 637       pass = kFALSE;
 638       for(Int_t k = 0; k < numberOfTriggerSelected; k++){
 639         const TObjString *const obString=(TObjString*)fSelectedTrigger->At(k);
 640         const TString tString=obString->GetString();
 641         if(fESD->IsTriggerClassFired((const char*)tString)) {
 642           pass = kTRUE;
 643         }
 644       }
 645     }
 646     if(!pass) {
 647       PostData(1, fListHist);
 648       return;
 649     }
 650
 651   }
 652
 653   TString classfired = fESD->GetFiredTriggerClasses();
 654   Bool_t isWU =  classfired.Contains("WU");
 655   if(isWU) fNEventsTrigger->Fill(1.5);
 656   else fNEventsTrigger->Fill(0.5);
 657
 658   //printf("Class Fired %s\n",(const char*)fESD->GetFiredTriggerClasses());
 659   //printf("Trigger passed\n");
 660
 661   ///////////////////////////////
 662   // Require a primary vertex
 663   //////////////////////////////
 664   if(fRequirePrimaryVertex) {
 665     const AliESDVertex* vtxESD = 0x0;
 666     if      (fVtxTPC) vtxESD = fESD->GetPrimaryVertexTPC() ;
 667     else if (fVtxSPD) vtxESD = fESD->GetPrimaryVertexSPD() ;
 668     else              vtxESD = fESD->GetPrimaryVertexTracks() ;
 669     if(!vtxESD){
 670       PostData(1, fListHist);
 671       return;
 672     }
 673     Int_t nCtrb = vtxESD->GetNContributors();
 674     if(nCtrb < fMinNbContributors) {
 675       PostData(1, fListHist);
 676       return;
 677     }
 678     Double_t zPosition = vtxESD->GetZ();
 679     if(TMath::Abs(zPosition) > fRangePrimaryVertexZ) {
 680       PostData(1, fListHist);
 681       return;
 682     }
 683
 684   }
 685
 686   //printf("Primary vertex passed\n");
 687
 688   //////////////////////////////////
 689   // Reject pile-up with SPD
 690   //////////////////////////////////
 691
 692   if(fRejectPileUpWithSPD) {
 693     if(fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5)){
 694       //printf("test\n");
 695       PostData(1, fListHist);
 696       return;
 697     }
 698   }
 699
 700   //////////////////////////////////////
 701   // Requirement on number of good tracks
 702   //////////////////////////////////////
 703   Int_t nGoodParticles = 0;
 704   Double_t nbTracks = fESD->GetNumberOfTracks();
 705   for(Int_t itrack = 0; itrack < nbTracks; itrack++) {
 706     if(ParticleGood(itrack)) nGoodParticles++;
 707   }
 708   if(fDebug > 1)  {
 709     // Centrality
 710     AliCentrality *esdCentrality = fESD->GetCentrality();
 711     Float_t centrality = esdCentrality->GetCentralityPercentile("V0M");
 712     //Float_t centralityb = esdCentrality->GetCentralityPercentile("CL1");
 713     fNbGoodTracks->Fill(nGoodParticles,centrality);
 714     //printf("centrality %f, centralityb %f\n",centrality,centralityb);
 715   }
 716
 717   //printf("Beam type %s\n",(const char*)type);
 718   if (!strstr(type,"p-p")) {
 719     //if (strstr(type,"A-A")) {
 720     //printf("Will check the number of good tracks %f %f %f\n",nGoodParticles,fMinNbTracks,fMaxNbTracks);
 721     if((nGoodParticles < fMinNbTracks) || (nGoodParticles > fMaxNbTracks)) {
 722       PostData(1, fListHist);
 723       return;
 724     }
 725   }
 726
 727   fNEvents->Fill(1);
 728
 729   // In total
 730   Int_t nbTrdTracks = 0;
 731   // standalone
 732   Int_t nbTrdTracksStandalone = 0;
 733   // offline
 734   Int_t nbTrdTracksOffline = 0;
 735   // TPC
 736   Int_t nbtrackTPC = 0;
 737
 738
 739
 740   if (nbTracks <= 0.0) {
 741
 742     if(fDebug > 1) {
 743       fNbTRDTrack->Fill(nbTrdTracks);
 744       fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
 745       fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
 746     }
 747     PostData(1, fListHist);
 748     return;
 749   }
 750
 751
 752   fESDfriend = dynamic_cast<AliESDfriend*> (fESD->FindListObject("AliESDfriend"));
 753   if(!fESDfriend){
 754     AliError("fESDfriend not available");
 755     PostData(1, fListHist);
 756     return;
 757   }
 758
 759   if(fESDfriend->TestSkipBit()) {
 760     PostData(1, fListHist);
 761     return;
 762   }
 763
 764   //printf("has friends\n");
 765
 766   /////////////////////////////////////
 767   // Loop on AliESDtrack
 768   ////////////////////////////////////
 769   //printf("Nb of tracks %f\n",nbTracks);
 770   Int_t nbTracksfriends = fESDfriend->GetNumberOfTracks();
 771   for(int itrk=0; itrk < nbTracksfriends; ++itrk){
 772
 773     // Get ESD track
 774     fkEsdTrack = fESD->GetTrack(itrk);
 775     if(!fkEsdTrack) continue;
 776     ULong_t status = fkEsdTrack->GetStatus();
 777     if(status&(AliESDtrack::kTPCout)) ++nbtrackTPC;
 778
 779     // Fix suggested by Alex
 780     fFriendTrack = fESDfriend->GetTrack(itrk);
 781     //printf("itrk %d\n",itrk);
 782     //fFriendTrack = (fESDfriend->GetNumberOfTracks()>itrk)?fESDfriend->GetTrack(itrk):NULL;
 783     if(!fFriendTrack)  {
 784       //printf("No friend track %d\n",itrk);
 785       continue;
 786     }
 787
 788     // Other cuts
 789     fTrdTrack = 0x0;
 790     Bool_t good = kTRUE;
 791     Bool_t standalonetrack = kFALSE;
 792     Bool_t offlinetrack = kFALSE;
 793     //ULong_t status = fkEsdTrack->GetStatus();
 794
 795     //////////////////////////////////////
 796     // Loop on calibration objects
 797     //////////////////////////////////////
 798     Int_t icalib=0;
 799     Int_t nTRDtrackV1=0;
 800     while((fCalibObject = (TObject *)(fFriendTrack->GetCalibObject(icalib++)))){
 801       //printf("Name %s\n",fCalibObject->IsA()->GetName());
 802       if(strcmp(fCalibObject->IsA()->GetName(), "AliTRDtrackV1") != 0) continue;
 803       //printf("Find the calibration object\n");
 804       ++nTRDtrackV1;
 805
 806       if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {
 807         standalonetrack = kTRUE;
 808       }
 809       if((status&(AliESDtrack::kTRDin))) {
 810         offlinetrack = kTRUE;
 811       }
 812       if(fOfflineTracks){
 813         if(!offlinetrack){
 814           good = kFALSE;
 815         }
 816       }
 817       else if(fStandaloneTracks){
 818         if(!standalonetrack){
 819           good = kFALSE;
 820         }
 821       }
 822
 823       fTrdTrack = (AliTRDtrackV1 *)fCalibObject;
 824       // process chamberstatus
 825       fTRDChamberStatus->ProcessTrack(fTrdTrack);
 826     }
 827
 828     // Quality cuts on the AliESDtrack
 829     if((fEsdTrackCuts) && (!fEsdTrackCuts->IsSelected((AliVParticle *)fkEsdTrack))) {
 830       //printf("Not a good track\n");
 831       continue;
 832     }
 833
 834     Int_t nbcrossing = fkEsdTrack->GetTOFBunchCrossing();
 835     fTOFbc->Fill(nbcrossing,fkEsdTrack->Pt());
 836     // TOF pile-up rejection is asked
 837     if(fRejectPileUpWithTOF) {
 838       if(TMath::Abs(nbcrossing)>0.5) continue;
 839     }
 840
 841     // ITS or TOF
 842     if(fRejectPileUpWithTOFOrITS) {
 843       ULong_t statusits = fkEsdTrack->GetStatus();
 844       UChar_t itsPixel = fkEsdTrack->GetITSClusterMap();
 845       Bool_t itskany = kFALSE;
 846       if(((statusits & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) && ((TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1)))) itskany = kTRUE;
 847       if(!(itskany || (TMath::Abs(nbcrossing)<0.5))) continue;
 848     }
 849
 850     // First Absolute gain calibration
 851     Int_t trdNTracklets = (Int_t) fkEsdTrack->GetTRDntracklets();
 852     Int_t trdNTrackletsPID = (Int_t) fkEsdTrack->GetTRDntrackletsPID();
 853     //printf("Number of trd tracklets %d and PID trd tracklets %d\n",trdNTracklets,trdNTrackletsPID);
 854     if((trdNTracklets > 0) && (trdNTrackletsPID > 0)) {
 855       for(Int_t iPlane = 0; iPlane < 6; ++iPlane){
 856         //Double_t slide = fkEsdTrack->GetTRDslice(iPlane);
 857         //printf("Number of slide %d\n",fkEsdTrack->GetNumberOfTRDslices());
 858         //Double_t momentum = fkEsdTrack->GetTRDmomentum(iPlane);
 859         //printf("momentum %f, slide %f\n",momentum,slide);
 860         if(fkEsdTrack->GetTRDslice(iPlane) > 0.0)
 861           fAbsoluteGain->Fill(fkEsdTrack->GetTRDslice(iPlane)*8.0/100.0,
 862                               fkEsdTrack->GetTRDmomentum(iPlane));
 863       }
 864     }
 865
 866
 867     if(!fTrdTrack) continue;
 868
 869     if(good && fOnInstance) {
 870       //cout << "good" << endl;
 871       fTRDCalibraFillHisto->UpdateHistogramsV1(fTrdTrack,fkEsdTrack);
 872       //printf("Fill fTRDCalibraFillHisto\n");
 873     }
 874
 875     if(IsPHQon()){
 876       const Double_t mag = AliTRDdEdxBaseUtils::IsExBOn() ? fESD->GetMagneticField() : -1;
 877       const Int_t charge = AliTRDdEdxBaseUtils::IsExBOn() ? fkEsdTrack->Charge() : -1;
 878       const Double_t toTPCscale = AliTRDdEdxCalibUtils::GetCalibTPCscale(fkEsdTrack->GetTPCncls(), fkEsdTrack->GetTPCsignal());
 879       if(toTPCscale>0){
 880         AliTRDdEdxCalibUtils::FillHist(fTrdTrack, 0, mag, charge, toTPCscale);
 881       }
 882     }
 883
 884     //////////////////////////////////
 885     // Debug
 886     ////////////////////////////////
 887
 888     if(fDebug > 0) {
 889
 890       //printf("Enter debug\n");
 891
 892       Int_t nbtracklets = 0;
 893
 894       // Check some stuff
 895       // Coverity
 896       //Bool_t standalonetracklet = kFALSE;
 897       const AliTRDseedV1 *tracklet = 0x0;
 898       //////////////////////////////////////
 899       // Loop tracklets
 900       /////////////////////////////////////
 901       Int_t nbclusters=0;
 902       Double_t phtb[AliTRDseedV1::kNtb];
 903       memset(phtb, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
 904       Double_t sum = 0.0;
 905       Float_t normalisation = 1.13;
 906       Int_t detector = 0;
 907       Int_t sector = 0;
 908       for(Int_t itr = 0; itr < 6; ++itr){
 909
 910         if(!(tracklet = fTrdTrack->GetTracklet(itr))) continue;
 911         if(!tracklet->IsOK()) continue;
 912         ++nbtracklets;
 913         // Coverity
 914         //standalonetracklet = kFALSE;
 915         //if(tracklet->IsStandAlone()) standalonetracklet = kTRUE;
 916
 917         nbclusters = 0;
 918         memset(phtb, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
 919         sum = 0.0;
 920         detector = 0;
 921         sector = 0;
 922         //Int_t crossrow = 0;
 923
 924           // Check no shared clusters
 925           //for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
 926           //  if((fcl = tracklet->GetClusters(icc)))  crossrow = 1;
 927           // }
 928
 929           // Loop on clusters
 930         Int_t time = 0;
 931         Float_t ch = 0;
 932         Float_t qcl = 0;
 933         for(int ic=0; ic<AliTRDseedV1::kNtb; ++ic){
 934
 935           if(!(fCl = tracklet->GetClusters(ic))) continue;
 936           ++nbclusters;
 937           time = fCl->GetPadTime();
 938           //ch =  tracklet->GetdQdl(ic);
 939           ch =  tracklet->GetQperTB(ic);
 940           qcl = TMath::Abs(fCl->GetQ());
 941           detector = fCl->GetDetector();
 942           // Add the charge if shared cluster
 943           if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
 944             if((fCl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) {
 945               qcl += TMath::Abs(fCl->GetQ());
 946               //printf("Add the cluster charge\n");
 947             }
 948           }
 949           if((time>-1) && (time<fNbTimeBins)) phtb[time]=qcl;
 950           if((fCalDetGain) && (fCalDetGain->GetValue(detector) > 0.0)) sum += ch*fCalDetGain->GetValue(detector)/normalisation;
 951           else sum += ch/normalisation;
 952
 953           if(fDebug > 1) {
 954             fNbTimeBin->Fill(time);
 955             if(tracklet->IsStandAlone()) fNbTimeBinStandalone->Fill(time);
 956             else fNbTimeBinOffline->Fill(time);
 957           }
 958         }
 959         sector = AliTRDgeometry::GetSector(detector);
 960
 961         if(fDebug > 1) {
 962           fNbTracklets->Fill(detector);
 963           if(tracklet->IsStandAlone()) fNbTrackletsStandalone->Fill(detector);
 964           else fNbTrackletsOffline->Fill(detector);
 965
 966           fNbClusters->Fill(nbclusters);
 967           if(tracklet->IsStandAlone())  fNbClustersStandalone->Fill(nbclusters);
 968           else  fNbClustersOffline->Fill(nbclusters);
 969         }
 970
 971         if((nbclusters > fLow) && (nbclusters < fHigh)){
 972           if(fRelativeScale > 0.0) sum = sum/fRelativeScale;
 973           fCH2dTest->Fill(sum,detector+0.5);
 974           fCH2dSM->Fill(sum,sector+0.5);
 975           fCH2dSum->Fill(sum,0.5);
 976           Bool_t checknoise = kTRUE;
 977           if(fMaxCluster > 0) {
 978             if(phtb[0] > fMaxCluster) checknoise = kFALSE;
 979             if(fNbTimeBins > fNbMaxCluster) {
 980               for(Int_t k = (fNbTimeBins-fNbMaxCluster); k < fNbTimeBins; k++){
 981                 if(phtb[k] > fMaxCluster) checknoise = kFALSE;
 982               }
 983             }
 984           }
 985           if(checknoise) {
 986             for(int ic=0; ic<fNbTimeBins; ic++){
 987               if(fFillZero) {
 988                 fPH2dTest->Fill((Double_t)(ic/10.0),detector+0.5,(Double_t)phtb[ic]);
 989                 fPH2dSum->Fill((Double_t)(ic/10.0),0.5,(Double_t)phtb[ic]);
 990                 fPH2dSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
 991               }
 992               else {
 993                 if(phtb[ic] > 0.0) {
 994                   fPH2dTest->Fill((Double_t)(ic/10.0),detector+0.5,(Double_t)phtb[ic]);
 995                   fPH2dSum->Fill((Double_t)(ic/10.0),0.0,(Double_t)phtb[ic]);
 996                   fPH2dSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
 997                 }
 998               }
 999             }
1000           }
1001         }
1002         if(detector == 0) FindP1TrackPHtrackletV1Test(tracklet,nbclusters);
1003
1004       } // loop on tracklets
1005
1006
1007     } // debug
1008
1009     if(nTRDtrackV1 > 0) {
1010       ++nbTrdTracks;
1011       if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {
1012         ++nbTrdTracksStandalone;
1013       }
1014       if((status&(AliESDtrack::kTRDin))) {
1015         ++nbTrdTracksOffline;
1016       }
1017     }
1018     //delete fFriendTrack;
1019   } // loop ESD track
1020
1021   if(fDebug > 1) {
1022     fNbTRDTrack->Fill(nbTrdTracks);
1023     fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
1024     fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
1025     fNbTPCTRDtrack->Fill(nbTrdTracks,nbtrackTPC);
1026   }
1027
1028   // Post output data
1029   PostData(1, fListHist);
1030   //cout << "AliTRDCalibTask::Exec() OUT" << endl;
1031 }
1032
1033 //________________________________________________________________________
1034 void AliTRDCalibTask::Terminate(Option_t *)
1035 {
1036   //
1037   // Terminate
1038   //
1039
1040   if(fTRDCalibraFillHisto) fTRDCalibraFillHisto->DestroyDebugStreamer();
1041
1042
1043 }
1044 //_______________________________________________________
1045 Bool_t AliTRDCalibTask::Load(const Char_t *filename)
1046 {
1047   //
1048   // Generic container loader
1049   //
1050
1051   if(!TFile::Open(filename)){
1052     //AliWarning(Form("Couldn't open file %s.", filename));
1053     return kFALSE;
1054   }
1055   TList *o = 0x0;
1056   if(!(o = (TList*)gFile->Get(GetName()))){
1057     //AliWarning("Missing histogram container.");
1058     return kFALSE;
1059   }
1060   fListHist = (TList*)o->Clone(GetName());
1061   gFile->Close();
1062   return kTRUE;
1063 }
1064 //_______________________________________________________
1065 Bool_t AliTRDCalibTask::Load(TList *lister)
1066 {
1067   //
1068   // Generic container loader
1069   //
1070
1071   fListHist = (TList*)lister->Clone(GetName());
1072   return kTRUE;
1073 }
1074 //_______________________________________________________________________________________
1075 void  AliTRDCalibTask::AddTask(const AliTRDCalibTask * calibTask) {
1076
1077   //
1078   // Add stats
1079   //
1080
1081   TList *listcalibTask = calibTask->GetList();
1082   if(!listcalibTask) return;
1083
1084   THnSparseI *histoEntries = (THnSparseI *) listcalibTask->FindObject("NumberOfEntries");
1085
1086   TH1I *nEvents  = (TH1I *) listcalibTask->FindObject(Form("NEvents_%s",(const char*)calibTask->GetName()));
1087   TH1I *nEventsInput  = (TH1I *) listcalibTask->FindObject(Form("NEventsInput_%s",(const char*)calibTask->GetName()));
1088   TH2F *absoluteGain  = (TH2F *) listcalibTask->FindObject(Form("AbsoluteGain_%s",(const char*)calibTask->GetName()));
1089
1090   TH1F *trdTrack = (TH1F *) listcalibTask->FindObject(Form("TRDTrack_%s",(const char*)calibTask->GetName()));
1091   TH1F *trdTrackOffline = (TH1F *) listcalibTask->FindObject(Form("TRDTrackOffline_%s",(const char*)calibTask->GetName()));
1092   TH1F *trdTrackStandalone = (TH1F *) listcalibTask->FindObject(Form("TRDTrackStandalone_%s",(const char*)calibTask->GetName()));
1093
1094   TH2F *tpctrdTrack = (TH2F *) listcalibTask->FindObject(Form("NbTPCTRDtrack_%s",(const char*)calibTask->GetName()));
1095
1096   TH1F *nbTimeBin = (TH1F *) listcalibTask->FindObject(Form("NbTimeBin_%s",(const char*)calibTask->GetName()));
1097   TH1F *nbTimeBinOffline = (TH1F *) listcalibTask->FindObject(Form("NbTimeBinOffline_%s",(const char*)calibTask->GetName()));
1098   TH1F *nbTimeBinStandalone = (TH1F *) listcalibTask->FindObject(Form("NbTimeBinStandalone_%s",(const char*)calibTask->GetName()));
1099
1100   TH1F *nbClusters = (TH1F *) listcalibTask->FindObject(Form("NbClusters_%s",(const char*)calibTask->GetName()));
1101   TH1F *nbClustersOffline = (TH1F *) listcalibTask->FindObject(Form("NbClustersOffline_%s",(const char*)calibTask->GetName()));
1102   TH1F *nbClustersStandalone = (TH1F *) listcalibTask->FindObject(Form("NbClustersStandalone_%s",(const char*)calibTask->GetName()));
1103
1104   TH1F *nbTracklets = (TH1F *) listcalibTask->FindObject(Form("NbTracklets_%s",(const char*)calibTask->GetName()));
1105   TH1F *nbTrackletsOffline = (TH1F *) listcalibTask->FindObject(Form("NbTrackletsOffline_%s",(const char*)calibTask->GetName()));
1106   TH1F *nbTrackletsStandalone = (TH1F *) listcalibTask->FindObject(Form("NbTrackletsStandalone_%s",(const char*)calibTask->GetName()));
1107
1108   TH2I *ch2d = (TH2I *) listcalibTask->FindObject("CH2d");
1109   TProfile2D *ph2d = (TProfile2D *) listcalibTask->FindObject("PH2d");
1110   TProfile2D *prf2d = (TProfile2D *) listcalibTask->FindObject("PRF2d");
1111
1112   TH2I *ch2dSum = (TH2I *) listcalibTask->FindObject(Form("CH2dSum_%s",(const char*)calibTask->GetName()));
1113   TProfile2D *ph2dSum = (TProfile2D *) listcalibTask->FindObject(Form("PH2dSum_%s",(const char*)calibTask->GetName()));
1114
1115   TH2I *ch2dSM = (TH2I *) listcalibTask->FindObject(Form("CH2dSM_%s",(const char*)calibTask->GetName()));
1116   TProfile2D *ph2dSM = (TProfile2D *) listcalibTask->FindObject(Form("PH2dSM_%s",(const char*)calibTask->GetName()));
1117
1118   AliTRDCalibraVdriftLinearFit *linearfit = (AliTRDCalibraVdriftLinearFit *) listcalibTask->FindObject("AliTRDCalibraVdriftLinearFit");
1119   AliTRDCalibraExbAltFit *exbaltfit = (AliTRDCalibraExbAltFit *) listcalibTask->FindObject("AliTRDCalibraExbAltFit");
1120   AliTRDCalibraVector *calibraVector = (AliTRDCalibraVector *) listcalibTask->FindObject("AliTRDCalibraVector");
1121
1122   //
1123
1124   THnSparseI *inhistoEntries = (THnSparseI *) fListHist->FindObject("NumberOfEntries");
1125
1126   TH1I *inEventsInput  = (TH1I *) fListHist->FindObject(Form("NEventsInput_%s",(const char*)fName));
1127   TH1I *inEvents  = (TH1I *) fListHist->FindObject(Form("NEvents_%s",(const char*)fName));
1128   TH2F *iabsoluteGain  = (TH2F *) fListHist->FindObject(Form("AbsoluteGain_%s",(const char*)fName));
1129
1130   TH1F *itrdTrack = (TH1F *) fListHist->FindObject(Form("TRDTrack_%s",(const char*)fName));
1131   TH1F *itrdTrackOffline = (TH1F *) fListHist->FindObject(Form("TRDTrackOffline_%s",(const char*)fName));
1132   TH1F *itrdTrackStandalone = (TH1F *) fListHist->FindObject(Form("TRDTrackStandalone_%s",(const char*)fName));
1133
1134   TH2F *itpctrdTrack = (TH2F *) fListHist->FindObject(Form("NbTPCTRDtrack_%s",(const char*)fName));
1135
1136   TH1F *inbTimeBin = (TH1F *) fListHist->FindObject(Form("NbTimeBin_%s",(const char*)fName));
1137   TH1F *inbTimeBinOffline = (TH1F *) fListHist->FindObject(Form("NbTimeBinOffline_%s",(const char*)fName));
1138   TH1F *inbTimeBinStandalone = (TH1F *) fListHist->FindObject(Form("NbTimeBinStandalone_%s",(const char*)fName));
1139
1140   TH1F *inbClusters = (TH1F *) fListHist->FindObject(Form("NbClusters_%s",(const char*)fName));
1141   TH1F *inbClustersOffline = (TH1F *) fListHist->FindObject(Form("NbClustersOffline_%s",(const char*)fName));
1142   TH1F *inbClustersStandalone = (TH1F *) fListHist->FindObject(Form("NbClustersStandalone_%s",(const char*)fName));
1143
1144   TH1F *inbTracklets = (TH1F *) fListHist->FindObject(Form("NbTracklets_%s",(const char*)fName));
1145   TH1F *inbTrackletsOffline = (TH1F *) fListHist->FindObject(Form("NbTrackletsOffline_%s",(const char*)fName));
1146   TH1F *inbTrackletsStandalone = (TH1F *) fListHist->FindObject(Form("NbTrackletsStandalone_%s",(const char*)fName));
1147
1148   TH2I *ich2d = (TH2I *) fListHist->FindObject("CH2d");
1149   TProfile2D *iph2d = (TProfile2D *) fListHist->FindObject("PH2d");
1150   TProfile2D *iprf2d = (TProfile2D *) fListHist->FindObject("PRF2d");
1151
1152   TH2I *ich2dSum = (TH2I *) fListHist->FindObject(Form("CH2dSum_%s",(const char*)fName));
1153   TProfile2D *iph2dSum = (TProfile2D *) fListHist->FindObject(Form("PH2dSum_%s",(const char*)fName));
1154
1155   TH2I *ich2dSM = (TH2I *) fListHist->FindObject(Form("CH2dSM_%s",(const char*)fName));
1156   TProfile2D *iph2dSM = (TProfile2D *) fListHist->FindObject(Form("PH2dSM_%s",(const char*)fName));
1157
1158   AliTRDCalibraVdriftLinearFit *ilinearfit = (AliTRDCalibraVdriftLinearFit *) fListHist->FindObject("AliTRDCalibraVdriftLinearFit");
1159   AliTRDCalibraExbAltFit *iexbaltfit = (AliTRDCalibraExbAltFit *) fListHist->FindObject("AliTRDCalibraExbAltFit");
1160   AliTRDCalibraVector *icalibraVector = (AliTRDCalibraVector *) fListHist->FindObject("AliTRDCalibraVector");
1161
1162
1163   // Add
1164
1165   if(histoEntries) {
1166     if(inhistoEntries) {
1167       inhistoEntries->Add(histoEntries);
1168       //printf("Add Events\n");
1169     }
1170     else {
1171       //printf("Create new Events\n");
1172       inhistoEntries = (THnSparseI *) histoEntries->Clone();
1173       fListHist->Add(inhistoEntries);
1174     }
1175   }
1176
1177   if(nEventsInput) {
1178     if(inEventsInput) {
1179       inEventsInput->Add(nEventsInput);
1180       //printf("Add Events\n");
1181     }
1182     else {
1183       //printf("Create new Events\n");
1184       inEventsInput = new TH1I(*nEventsInput);
1185       fListHist->Add(inEventsInput);
1186     }
1187   }
1188
1189   if(nEvents) {
1190     if(inEvents) {
1191       inEvents->Add(nEvents);
1192       //printf("Add Events\n");
1193     }
1194     else {
1195       //printf("Create new Events\n");
1196       inEvents = new TH1I(*nEvents);
1197       fListHist->Add(inEvents);
1198     }
1199   }
1200
1201   if(absoluteGain) {
1202     if(iabsoluteGain) iabsoluteGain->Add(absoluteGain);
1203     else {
1204       iabsoluteGain = new TH2F(*absoluteGain);
1205       fListHist->Add(iabsoluteGain);
1206     }
1207   }
1208
1209   if(trdTrack) {
1210     if(itrdTrack) itrdTrack->Add(trdTrack);
1211     else {
1212      itrdTrack = new TH1F(*trdTrack);
1213      fListHist->Add(itrdTrack);
1214     }
1215   }
1216
1217   if(trdTrackOffline) {
1218     if(itrdTrackOffline) itrdTrackOffline->Add(trdTrackOffline);
1219     else {
1220       itrdTrackOffline = new TH1F(*trdTrackOffline);
1221       fListHist->Add(itrdTrackOffline);
1222     }
1223   }
1224
1225   if(trdTrackStandalone) {
1226     if(itrdTrackStandalone) itrdTrackStandalone->Add(trdTrackStandalone);
1227     else {
1228       itrdTrackStandalone = new TH1F(*trdTrackStandalone);
1229       fListHist->Add(itrdTrackStandalone);
1230     }
1231   }
1232
1233   if(tpctrdTrack) {
1234     if(itpctrdTrack) itpctrdTrack->Add(tpctrdTrack);
1235     else {
1236       itpctrdTrack = new TH2F(*tpctrdTrack);
1237       fListHist->Add(itpctrdTrack);
1238     }
1239   }
1240
1241   if(nbTimeBin) {
1242     if(inbTimeBin) inbTimeBin->Add(nbTimeBin);
1243     else {
1244       inbTimeBin = new TH1F(*inbTimeBin);
1245       fListHist->Add(inbTimeBin);
1246     }
1247   }
1248
1249   if(nbTimeBinOffline) {
1250     if(inbTimeBinOffline) inbTimeBinOffline->Add(nbTimeBinOffline);
1251     else {
1252       inbTimeBinOffline = new TH1F(*nbTimeBinOffline);
1253       fListHist->Add(inbTimeBinOffline);
1254     }
1255   }
1256
1257   if(nbTimeBinStandalone) {
1258     if(inbTimeBinStandalone) inbTimeBinStandalone->Add(nbTimeBinStandalone);
1259     else {
1260       inbTimeBinStandalone = new TH1F(*nbTimeBinStandalone);
1261       fListHist->Add(inbTimeBinStandalone);
1262     }
1263   }
1264
1265   if(nbClusters) {
1266     if(inbClusters) inbClusters->Add(nbClusters);
1267     else {
1268       inbClusters = new TH1F(*nbClusters);
1269       fListHist->Add(inbClusters);
1270     }
1271   }
1272
1273   if(nbClustersOffline) {
1274     if(inbClustersOffline) inbClustersOffline->Add(nbClustersOffline);
1275     else {
1276       inbClustersOffline = new TH1F(*nbClustersOffline);
1277       fListHist->Add(inbClustersOffline);
1278     }
1279   }
1280
1281   if(nbClustersStandalone) {
1282     if(inbClustersStandalone) inbClustersStandalone->Add(nbClustersStandalone);
1283     else {
1284       inbClustersStandalone = new TH1F(*nbClustersStandalone);
1285       fListHist->Add(inbClustersStandalone);
1286     }
1287   }
1288
1289   if(nbTracklets) {
1290     if(inbTracklets) inbTracklets->Add(nbTracklets);
1291     else {
1292       inbTracklets = new TH1F(*nbTracklets);
1293       fListHist->Add(inbTracklets);
1294     }
1295   }
1296
1297   if(nbTrackletsOffline) {
1298     if(inbTrackletsOffline) inbTrackletsOffline->Add(nbTrackletsOffline);
1299     else {
1300       inbTrackletsOffline = new TH1F(*nbTrackletsOffline);
1301       fListHist->Add(inbTrackletsOffline);
1302     }
1303   }
1304
1305   if(nbTrackletsStandalone) {
1306     if(inbTrackletsStandalone) inbTrackletsStandalone->Add(nbTrackletsStandalone);
1307     else {
1308       inbTrackletsStandalone = new TH1F(*nbTrackletsStandalone);
1309       fListHist->Add(inbTrackletsStandalone);
1310     }
1311   }
1312
1313   if(ch2d) {
1314     if(ich2d) ich2d->Add(ch2d);
1315     else {
1316       ich2d = new TH2I(*ch2d);
1317       fListHist->Add(ich2d);
1318     }
1319   }
1320
1321   if(ph2d) {
1322     if(iph2d) iph2d->Add(ph2d);
1323     else {
1324       iph2d = new TProfile2D(*ph2d);
1325       fListHist->Add(iph2d);
1326     }
1327   }
1328
1329   if(prf2d) {
1330     if(iprf2d) iprf2d->Add(prf2d);
1331     else {
1332       iprf2d = new TProfile2D(*prf2d);
1333       fListHist->Add(iprf2d);
1334     }
1335   }
1336
1337   if(ch2dSum) {
1338     if(ich2dSum) ich2dSum->Add(ch2dSum);
1339     else {
1340       ich2dSum = new TH2I(*ch2dSum);
1341       fListHist->Add(ich2dSum);
1342     }
1343   }
1344
1345   if(ph2dSum) {
1346     if(iph2dSum) iph2dSum->Add(ph2dSum);
1347     else {
1348       iph2dSum = new TProfile2D(*ph2dSum);
1349       fListHist->Add(iph2dSum);
1350     }
1351   }
1352
1353   if(ch2dSM) {
1354     if(ich2dSM) ich2dSM->Add(ch2dSM);
1355     else {
1356       ich2dSM = new TH2I(*ch2dSM);
1357       fListHist->Add(ich2dSM);
1358     }
1359   }
1360
1361   if(ph2dSM) {
1362     if(iph2dSM) iph2dSM->Add(ph2dSM);
1363     else {
1364       iph2dSM = new TProfile2D(*ph2dSM);
1365       fListHist->Add(iph2dSM);
1366     }
1367   }
1368
1369   if(linearfit) {
1370     if(ilinearfit) ilinearfit->Add(linearfit);
1371     else {
1372       ilinearfit = new AliTRDCalibraVdriftLinearFit(*linearfit);
1373       fListHist->Add(ilinearfit);
1374     }
1375   }
1376
1377   if(exbaltfit) {
1378     if(iexbaltfit) iexbaltfit->Add(exbaltfit);
1379     else {
1380       iexbaltfit = new AliTRDCalibraExbAltFit(*exbaltfit);
1381       fListHist->Add(iexbaltfit);
1382     }
1383   }
1384
1385   if(calibraVector) {
1386     if(icalibraVector) icalibraVector->Add(calibraVector);
1387     else {
1388       icalibraVector = new AliTRDCalibraVector(*calibraVector);
1389       fListHist->Add(icalibraVector);
1390     }
1391   }
1392
1393 }
1394 //________________________________________________________________________________
1395 Long64_t AliTRDCalibTask::Merge(TCollection *li) {
1396
1397   //
1398   // merge component
1399   //
1400
1401   TIterator* iter = li->MakeIterator();
1402   AliTRDCalibTask* cal = 0;
1403
1404   while ((cal = (AliTRDCalibTask*)iter->Next())) {
1405     if (!cal->InheritsFrom(AliTRDCalibTask::Class())) {
1406       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
1407       return -1;
1408     }
1409
1410     // add histograms here...
1411     this->AddTask(cal);
1412
1413   }
1414
1415   return 0;
1416
1417 }
1418 //_____________________________________________________
1419 Bool_t AliTRDCalibTask::SetVersionSubversion(){
1420   //
1421   // Load Chamber Gain factors into the Tender supply
1422   //
1423
1424   //printf("SetVersionSubversion\n");
1425
1426   //find previous entry from the UserInfo
1427   TTree *tree=((TChain*)GetInputData(0))->GetTree();
1428   if (!tree) {
1429     AliError("Tree not found in ESDhandler");
1430     return kFALSE;
1431   }
1432
1433   TList *userInfo=(TList*)tree->GetUserInfo();
1434   if (!userInfo) {
1435     AliError("No UserInfo found in tree");
1436     return kFALSE;
1437   }
1438
1439   TList *cdbList=(TList*)userInfo->FindObject("cdbList");
1440   if (!cdbList) {
1441     AliError("No cdbList found in UserInfo");
1442     if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print();
1443     return kFALSE;
1444   }
1445
1446   TIter nextCDB(cdbList);
1447   TObjString *os=0x0;
1448   while ( (os=(TObjString*)nextCDB()) ){
1449     if(os->GetString().Contains("TRD/Calib/ChamberGainFactor")){
1450       // Get Old gain calibration
1451       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1452       fFirstRunGain = id->GetFirstRun();
1453       fVersionGainUsed = id->GetVersion();
1454       fSubVersionGainUsed = id->GetSubVersion();
1455     } else if(os->GetString().Contains("TRD/Calib/ChamberVdrift")){
1456       // Get Old drift velocity calibration
1457       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1458       fFirstRunVdrift = id->GetFirstRun();
1459       fVersionVdriftUsed = id->GetVersion();
1460       fSubVersionVdriftUsed = id->GetSubVersion();
1461     } else if(os->GetString().Contains("TRD/Calib/LocalGainFactor")){
1462       // Get Old drift velocity calibration
1463       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1464       fFirstRunGainLocal = id->GetFirstRun();
1465       fVersionGainLocalUsed = id->GetVersion();
1466       fSubVersionGainLocalUsed = id->GetSubVersion();
1467     } else if((os->GetString().Contains("TRD/Calib/ChamberExB")) && (!os->GetString().Contains("TRD/Calib/ChamberExBAlt"))){
1468       // Get Old drift velocity calibration
1469       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1470       fFirstRunExB = id->GetFirstRun();
1471       fVersionExBUsed = id->GetVersion();
1472       fSubVersionExBUsed = id->GetSubVersion();
1473       //printf("Version %d and subversion %d\n",fVersionExBUsed,fSubVersionExBUsed);
1474     }
1475   }
1476
1477   //printf("VersionGain %d, SubversionGain %d, VersionLocalGain %d, Subversionlocalgain %d, Versionvdrift %d, Subversionvdrift %d\n",fVersionGainUsed,fSubVersionGainUsed,fVersionGainLocalUsed,fSubVersionGainLocalUsed,fVersionVdriftUsed,fSubVersionVdriftUsed);
1478
1479   // Check
1480   if((fFirstRunGain < 0)            ||
1481      (fFirstRunGainLocal < 0)       ||
1482      (fFirstRunVdrift < 0)          ||
1483      (fVersionGainUsed < 0)         ||
1484      (fVersionGainLocalUsed < 0)    ||
1485      (fSubVersionGainUsed < 0)      ||
1486      (fSubVersionGainLocalUsed < 0) ||
1487      (fVersionVdriftUsed < 0)       ||
1488      (fSubVersionVdriftUsed < 0)) {
1489     AliError("No recent calibration found");
1490     return kFALSE;
1491   }
1492   else return kTRUE;
1493
1494 }
1495 //_________________________________________________________________________________________________________________________
1496 Bool_t AliTRDCalibTask::ParticleGood(int i) const {
1497
1498   //
1499   // Definition of good tracks
1500   //
1501
1502
1503   AliESDtrack *track = fESD->GetTrack(i);
1504   if (!track->IsOn(AliESDtrack::kTPCrefit)) return 0;        // TPC refit
1505   if (track->GetTPCNcls() < 90) return 0;                    // number of TPC clusters
1506   if (fabs(track->Eta())>0.8) return 0;                         // fiducial pseudorapidity
1507   Float_t r,z;
1508   track->GetImpactParametersTPC(r,z);
1509   if (fabs(z)>2.0) return 0;                          // impact parameter in z
1510   if (fabs(r)>2.0) return 0;                          // impact parameter in xy
1511   if (r==0) return 0;
1512   return 1;
1513
1514
1515 }
1516 //______________________________________________________________________________________________________________________
1517 Bool_t AliTRDCalibTask::FindP1TrackPHtrackletV1Test(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1518 {
1519   //
1520   // Drift velocity calibration:
1521   // Fit the clusters with a straight line
1522   // From the slope find the drift velocity
1523   //
1524
1525   ////////////////////////////////////////////////
1526   //Number of points: if less than 3 return kFALSE
1527   /////////////////////////////////////////////////
1528   if(nbclusters <= 2) return kFALSE;
1529
1530   ////////////
1531   //Variables
1532   ////////////
1533   // results of the linear fit
1534   Double_t dydt                       = 0.0;                                // dydt tracklet after straight line fit
1535   // Coverity
1536   //Double_t errorpar                   = 0.0;                                // error after straight line fit on dy/dt
1537   Double_t pointError                 = 0.0;                                // error after straight line fit
1538   // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
1539   Int_t    crossrow                   = 0;                                  // if it crosses a pad row
1540   Int_t    rowp                       = -1;                                 // if it crosses a pad row
1541   Float_t  tnt                        = tracklet->GetTilt();                // tan tiltingangle
1542   TLinearFitter linearFitterTracklet(2,"pol1");
1543   linearFitterTracklet.StoreData(kTRUE);
1544
1545
1546   ///////////////////////////////////////////
1547   // Take the parameters of the track
1548   //////////////////////////////////////////
1549   // take now the snp, tnp and tgl from the track
1550   Double_t snp = tracklet->GetSnp();             // sin dy/dx at the end of the chamber
1551   Double_t tnp = 0.0;                            // dy/dx at the end of the chamber
1552   if( TMath::Abs(snp) <  1.){
1553     tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1554   }
1555   Double_t tgl  = tracklet->GetTgl();           // dz/dl
1556   Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp);   // dz/dx calculated from dz/dl
1557   // at the entrance
1558   //Double_t tnp = tracklet->GetYref(1);      // dy/dx at the entrance of the chamber
1559   //Double_t tgl = tracklet->GetZref(1);      // dz/dl at the entrance of the chamber
1560   //Double_t dzdx = tgl;                      //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1561   // at the end with correction due to linear fit
1562   //Double_t tnp = tracklet->GetYfit(1);      // dy/dx at the end of the chamber after fit correction
1563   //Double_t tgl = tracklet->GetZfit(1);      // dz/dl at the end of the chamber after fit correction
1564
1565
1566   ////////////////////////////
1567   // loop over the clusters
1568   ////////////////////////////
1569   Int_t  nbli = 0;
1570   AliTRDcluster *cl                   = 0x0;
1571   //////////////////////////////
1572   // Check no shared clusters
1573   //////////////////////////////
1574   for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
1575     cl = tracklet->GetClusters(icc);
1576     if(cl)  crossrow = 1;
1577   }
1578   //////////////////////////////////
1579   // Loop clusters
1580   //////////////////////////////////
1581   for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1582     if(!(cl = tracklet->GetClusters(ic))) continue;
1583     //if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
1584
1585     Double_t ycluster                 = cl->GetY();
1586     Int_t time                        = cl->GetPadTime();
1587     Double_t timeis                   = time/10.0;
1588     //See if cross two pad rows
1589     Int_t    row                      = cl->GetPadRow();
1590     if(rowp==-1) rowp                 = row;
1591     if(row != rowp) crossrow          = 1;
1592
1593     linearFitterTracklet.AddPoint(&timeis,ycluster,1);
1594     nbli++;
1595
1596
1597   }
1598
1599   ////////////////////////////////////
1600   // Do the straight line fit now
1601   ///////////////////////////////////
1602   if(nbli <= 2){
1603     linearFitterTracklet.ClearPoints();
1604     return kFALSE;
1605   }
1606   TVectorD pars;
1607   linearFitterTracklet.Eval();
1608   linearFitterTracklet.GetParameters(pars);
1609   pointError  =  TMath::Sqrt(linearFitterTracklet.GetChisquare()/(nbli-2));
1610   // Coverity
1611   //errorpar    =  linearFitterTracklet.GetParError(1)*pointError;
1612   dydt        = pars[1];
1613   //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",linearFitterTracklet->GetChisquare(),nbli,pointError,linearFitterTracklet->GetParError(1),errorpar);
1614   linearFitterTracklet.ClearPoints();
1615
1616   /////////////////////////
1617   // Cuts quality
1618   ////////////////////////
1619
1620   if(nbclusters < fLow) return kFALSE;
1621   if(nbclusters > fHigh) return kFALSE;
1622   if(pointError >= 0.3) return kFALSE;
1623   if(crossrow == 1) return kTRUE;
1624
1625   ///////////////////////
1626   // Fill
1627   //////////////////////
1628
1629   if(fDebug > 0){
1630     //Add to the linear fitter of the detector
1631     if( TMath::Abs(snp) <  1.){
1632       Double_t x = tnp-dzdx*tnt;
1633       //if(!fLinearVdriftTest) printf("Not there\n");
1634       Double_t nbentries = fLinearVdriftTest->GetEntries();
1635       if(nbentries < (5.0*32767)) fLinearVdriftTest->Fill(x,dydt);
1636     }
1637   }
1638
1639   return kTRUE;
1640 }
1641