TRD/AliTRDCalibTask.cxx

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