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