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