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