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