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   //printf("Beam type %s\n",(const char*)type);
 642   if (!strstr(type,"p-p")) {
 643     //if (strstr(type,"A-A")) {
 644     //printf("Will check the number of good tracks %f %f %f\n",nGoodParticles,fMinNbTracks,fMaxNbTracks);
 645     if((nGoodParticles < fMinNbTracks) || (nGoodParticles > fMaxNbTracks)) {
 646       PostData(1, fListHist);
 647       return;
 648     }
 649   }
 650
 651   fNEvents->Fill(1);
 652
 653   // In total
 654   Int_t nbTrdTracks = 0;
 655   // standalone
 656   Int_t nbTrdTracksStandalone = 0;
 657   // offline
 658   Int_t nbTrdTracksOffline = 0;
 659   // TPC
 660   Int_t nbtrackTPC = 0;
 661
 662
 663
 664   if (nbTracks <= 0.0) {
 665
 666     if(fDebug > 1) {
 667       fNbTRDTrack->Fill(nbTrdTracks);
 668       fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
 669       fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
 670     }
 671     PostData(1, fListHist);
 672     return;
 673   }
 674
 675
 676   fESDfriend = dynamic_cast<AliESDfriend*> (fESD->FindListObject("AliESDfriend"));
 677   if(!fESDfriend){
 678     AliError("fESDfriend not available");
 679     PostData(1, fListHist);
 680     return;
 681   }
 682
 683   if(fESDfriend->TestSkipBit()) {
 684     PostData(1, fListHist);
 685     return;
 686   }
 687
 688   //printf("has friends\n");
 689
 690   /////////////////////////////////////
 691   // Loop on AliESDtrack
 692   ////////////////////////////////////
 693   //printf("Nb of tracks %f\n",nbTracks);
 694   for(int itrk=0; itrk < nbTracks; ++itrk){
 695
 696     // Get ESD track
 697     fkEsdTrack = fESD->GetTrack(itrk);
 698     if(!fkEsdTrack) continue;
 699     ULong_t status = fkEsdTrack->GetStatus();
 700     if(status&(AliESDtrack::kTPCout)) ++nbtrackTPC;
 701
 702     fFriendTrack = fESDfriend->GetTrack(itrk);
 703     if(!fFriendTrack)  {
 704       //printf("No friend track %d\n",itrk);
 705       continue;
 706     }
 707
 708     // Other cuts
 709     fTrdTrack = 0x0;
 710     Bool_t good = kTRUE;
 711     Bool_t standalonetrack = kFALSE;
 712     Bool_t offlinetrack = kFALSE;
 713     //ULong_t status = fkEsdTrack->GetStatus();
 714
 715     //////////////////////////////////////
 716     // Loop on calibration objects
 717     //////////////////////////////////////
 718     Int_t icalib=0;
 719     Int_t nTRDtrackV1=0;
 720     while((fCalibObject = (TObject *)(fFriendTrack->GetCalibObject(icalib++)))){
 721       //printf("Name %s\n",fCalibObject->IsA()->GetName());
 722       if(strcmp(fCalibObject->IsA()->GetName(), "AliTRDtrackV1") != 0) continue;
 723       //printf("Find the calibration object\n");
 724       ++nTRDtrackV1;
 725
 726       if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {
 727         standalonetrack = kTRUE;
 728       }
 729       if((status&(AliESDtrack::kTRDin))) {
 730         offlinetrack = kTRUE;
 731       }
 732       if(fOfflineTracks){
 733         if(!offlinetrack){
 734           good = kFALSE;
 735         }
 736       }
 737       else if(fStandaloneTracks){
 738         if(!standalonetrack){
 739           good = kFALSE;
 740         }
 741       }
 742
 743       fTrdTrack = (AliTRDtrackV1 *)fCalibObject;
 744       // process chamberstatus
 745       fTRDChamberStatus->ProcessTrack(fTrdTrack);
 746     }
 747
 748     // Quality cuts on the AliESDtrack
 749     if((fEsdTrackCuts) && (!fEsdTrackCuts->IsSelected((AliVParticle *)fkEsdTrack))) {
 750       //printf("Not a good track\n");
 751       continue;
 752     }
 753
 754     // First Absolute gain calibration
 755     Int_t trdNTracklets = (Int_t) fkEsdTrack->GetTRDntracklets();
 756     Int_t trdNTrackletsPID = (Int_t) fkEsdTrack->GetTRDntrackletsPID();
 757     //printf("Number of trd tracklets %d and PID trd tracklets %d\n",trdNTracklets,trdNTrackletsPID);
 758     if((trdNTracklets > 0) && (trdNTrackletsPID > 0)) {
 759       for(Int_t iPlane = 0; iPlane < 6; ++iPlane){
 760         //Double_t slide = fkEsdTrack->GetTRDslice(iPlane);
 761         //printf("Number of slide %d\n",fkEsdTrack->GetNumberOfTRDslices());
 762         //Double_t momentum = fkEsdTrack->GetTRDmomentum(iPlane);
 763         //printf("momentum %f, slide %f\n",momentum,slide);
 764         if(fkEsdTrack->GetTRDslice(iPlane) > 0.0)
 765           fAbsoluteGain->Fill(fkEsdTrack->GetTRDslice(iPlane)*8.0/100.0,
 766                               fkEsdTrack->GetTRDmomentum(iPlane));
 767       }
 768     }
 769
 770
 771     if(!fTrdTrack) continue;
 772
 773     if(good && fOnInstance) {
 774       //cout << "good" << endl;
 775       fTRDCalibraFillHisto->UpdateHistogramsV1(fTrdTrack);
 776       //printf("Fill fTRDCalibraFillHisto\n");
 777     }
 778
 779
 780
 781     //////////////////////////////////
 782     // Debug
 783     ////////////////////////////////
 784
 785     if(fDebug > 0) {
 786
 787       //printf("Enter debug\n");
 788
 789       Int_t nbtracklets = 0;
 790
 791       // Check some stuff
 792       Bool_t standalonetracklet = kFALSE;
 793       const AliTRDseedV1 *tracklet = 0x0;
 794       //////////////////////////////////////
 795       // Loop tracklets
 796       /////////////////////////////////////
 797       Int_t nbclusters=0;
 798       Double_t phtb[AliTRDseedV1::kNtb];
 799       memset(phtb, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
 800       Double_t sum = 0.0;
 801       Float_t normalisation = 1.13;
 802       Int_t detector = 0;
 803       Int_t sector = 0;
 804       for(Int_t itr = 0; itr < 6; ++itr){
 805
 806         if(!(tracklet = fTrdTrack->GetTracklet(itr))) continue;
 807         if(!tracklet->IsOK()) continue;
 808         ++nbtracklets;
 809         standalonetracklet = kFALSE;
 810         if(tracklet->IsStandAlone()) standalonetracklet = kTRUE;
 811
 812         nbclusters = 0;
 813         memset(phtb, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
 814         sum = 0.0;
 815         detector = 0;
 816         sector = 0;
 817         //Int_t crossrow = 0;
 818
 819           // Check no shared clusters
 820           //for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
 821           //  if((fcl = tracklet->GetClusters(icc)))  crossrow = 1;
 822           // }
 823
 824           // Loop on clusters
 825         Int_t time = 0;
 826         Float_t ch = 0;
 827         Float_t qcl = 0;
 828         for(int ic=0; ic<AliTRDseedV1::kNtb; ++ic){
 829
 830           if(!(fCl = tracklet->GetClusters(ic))) continue;
 831           ++nbclusters;
 832           time = fCl->GetPadTime();
 833           //ch =  tracklet->GetdQdl(ic);
 834           ch =  tracklet->GetQperTB(ic);
 835           qcl = TMath::Abs(fCl->GetQ());
 836           detector = fCl->GetDetector();
 837           // Add the charge if shared cluster
 838           if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
 839             if((fCl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) {
 840               qcl += TMath::Abs(fCl->GetQ());
 841               //printf("Add the cluster charge\n");
 842             }
 843           }
 844           if((time>-1) && (time<fNbTimeBins)) phtb[time]=qcl;
 845           if((fCalDetGain) && (fCalDetGain->GetValue(detector) > 0.0)) sum += ch*fCalDetGain->GetValue(detector)/normalisation;
 846           else sum += ch/normalisation;
 847
 848           if(fDebug > 1) {
 849             fNbTimeBin->Fill(time);
 850             if(tracklet->IsStandAlone()) fNbTimeBinStandalone->Fill(time);
 851             else fNbTimeBinOffline->Fill(time);
 852           }
 853         }
 854         sector = AliTRDgeometry::GetSector(detector);
 855
 856         if(fDebug > 1) {
 857           fNbTracklets->Fill(detector);
 858           if(tracklet->IsStandAlone()) fNbTrackletsStandalone->Fill(detector);
 859           else fNbTrackletsOffline->Fill(detector);
 860
 861           fNbClusters->Fill(nbclusters);
 862           if(tracklet->IsStandAlone())  fNbClustersStandalone->Fill(nbclusters);
 863           else  fNbClustersOffline->Fill(nbclusters);
 864         }
 865
 866         if((nbclusters > fLow) && (nbclusters < fHigh)){
 867           if(fRelativeScale > 0.0) sum = sum/fRelativeScale;
 868           fCH2dTest->Fill(sum,detector+0.5);
 869           fCH2dSM->Fill(sum,sector+0.5);
 870           fCH2dSum->Fill(sum,0.5);
 871           Bool_t checknoise = kTRUE;
 872           if(fMaxCluster > 0) {
 873             if(phtb[0] > fMaxCluster) checknoise = kFALSE;
 874             if(fNbTimeBins > fNbMaxCluster) {
 875               for(Int_t k = (fNbTimeBins-fNbMaxCluster); k < fNbTimeBins; k++){
 876                 if(phtb[k] > fMaxCluster) checknoise = kFALSE;
 877               }
 878             }
 879           }
 880           if(checknoise) {
 881             for(int ic=0; ic<fNbTimeBins; ic++){
 882               if(fFillZero) {
 883                 fPH2dTest->Fill((Double_t)(ic/10.0),detector+0.5,(Double_t)phtb[ic]);
 884                 fPH2dSum->Fill((Double_t)(ic/10.0),0.5,(Double_t)phtb[ic]);
 885                 fPH2dSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
 886               }
 887               else {
 888                 if(phtb[ic] > 0.0) {
 889                   fPH2dTest->Fill((Double_t)(ic/10.0),detector+0.5,(Double_t)phtb[ic]);
 890                   fPH2dSum->Fill((Double_t)(ic/10.0),0.0,(Double_t)phtb[ic]);
 891                   fPH2dSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
 892                 }
 893               }
 894             }
 895           }
 896         }
 897         if(detector == 0) FindP1TrackPHtrackletV1Test(tracklet,nbclusters);
 898
 899       } // loop on tracklets
 900
 901
 902     } // debug
 903
 904     if(nTRDtrackV1 > 0) {
 905       ++nbTrdTracks;
 906       if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {
 907         ++nbTrdTracksStandalone;
 908       }
 909       if((status&(AliESDtrack::kTRDin))) {
 910         ++nbTrdTracksOffline;
 911       }
 912     }
 913     //delete fFriendTrack;
 914   } // loop ESD track
 915
 916   if(fDebug > 1) {
 917     fNbTRDTrack->Fill(nbTrdTracks);
 918     fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
 919     fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
 920     fNbTPCTRDtrack->Fill(nbTrdTracks,nbtrackTPC);
 921   }
 922
 923   // Post output data
 924   PostData(1, fListHist);
 925   //cout << "AliTRDCalibTask::Exec() OUT" << endl;
 926 }
 927
 928 //________________________________________________________________________
 929 void AliTRDCalibTask::Terminate(Option_t *)
 930 {
 931   //
 932   // Terminate
 933   //
 934
 935   if(fTRDCalibraFillHisto) fTRDCalibraFillHisto->DestroyDebugStreamer();
 936
 937
 938 }
 939 //_______________________________________________________
 940 Bool_t AliTRDCalibTask::Load(const Char_t *filename)
 941 {
 942   //
 943   // Generic container loader
 944   //
 945
 946   if(!TFile::Open(filename)){
 947     //AliWarning(Form("Couldn't open file %s.", filename));
 948     return kFALSE;
 949   }
 950   TList *o = 0x0;
 951   if(!(o = (TList*)gFile->Get(GetName()))){
 952     //AliWarning("Missing histogram container.");
 953     return kFALSE;
 954   }
 955   fListHist = (TList*)o->Clone(GetName());
 956   gFile->Close();
 957   return kTRUE;
 958 }
 959 //_______________________________________________________
 960 Bool_t AliTRDCalibTask::Load(TList *lister)
 961 {
 962   //
 963   // Generic container loader
 964   //
 965
 966   fListHist = (TList*)lister->Clone(GetName());
 967   return kTRUE;
 968 }
 969 //_______________________________________________________________________________________
 970 void  AliTRDCalibTask::AddTask(const AliTRDCalibTask * calibTask) {
 971
 972   //
 973   // Add stats
 974   //
 975
 976   TList *listcalibTask = calibTask->GetList();
 977   if(!listcalibTask) return;
 978
 979   THnSparseI *histoEntries = (THnSparseI *) listcalibTask->FindObject("NumberOfEntries");
 980
 981   TH1I *nEvents  = (TH1I *) listcalibTask->FindObject(Form("NEvents_%s",(const char*)calibTask->GetName()));
 982   TH1I *nEventsInput  = (TH1I *) listcalibTask->FindObject(Form("NEventsInput_%s",(const char*)calibTask->GetName()));
 983   TH2F *absoluteGain  = (TH2F *) listcalibTask->FindObject(Form("AbsoluteGain_%s",(const char*)calibTask->GetName()));
 984
 985   TH1F *trdTrack = (TH1F *) listcalibTask->FindObject(Form("TRDTrack_%s",(const char*)calibTask->GetName()));
 986   TH1F *trdTrackOffline = (TH1F *) listcalibTask->FindObject(Form("TRDTrackOffline_%s",(const char*)calibTask->GetName()));
 987   TH1F *trdTrackStandalone = (TH1F *) listcalibTask->FindObject(Form("TRDTrackStandalone_%s",(const char*)calibTask->GetName()));
 988
 989   TH2F *tpctrdTrack = (TH2F *) listcalibTask->FindObject(Form("NbTPCTRDtrack_%s",(const char*)calibTask->GetName()));
 990
 991   TH1F *nbTimeBin = (TH1F *) listcalibTask->FindObject(Form("NbTimeBin_%s",(const char*)calibTask->GetName()));
 992   TH1F *nbTimeBinOffline = (TH1F *) listcalibTask->FindObject(Form("NbTimeBinOffline_%s",(const char*)calibTask->GetName()));
 993   TH1F *nbTimeBinStandalone = (TH1F *) listcalibTask->FindObject(Form("NbTimeBinStandalone_%s",(const char*)calibTask->GetName()));
 994
 995   TH1F *nbClusters = (TH1F *) listcalibTask->FindObject(Form("NbClusters_%s",(const char*)calibTask->GetName()));
 996   TH1F *nbClustersOffline = (TH1F *) listcalibTask->FindObject(Form("NbClustersOffline_%s",(const char*)calibTask->GetName()));
 997   TH1F *nbClustersStandalone = (TH1F *) listcalibTask->FindObject(Form("NbClustersStandalone_%s",(const char*)calibTask->GetName()));
 998
 999   TH1F *nbTracklets = (TH1F *) listcalibTask->FindObject(Form("NbTracklets_%s",(const char*)calibTask->GetName()));
1000   TH1F *nbTrackletsOffline = (TH1F *) listcalibTask->FindObject(Form("NbTrackletsOffline_%s",(const char*)calibTask->GetName()));
1001   TH1F *nbTrackletsStandalone = (TH1F *) listcalibTask->FindObject(Form("NbTrackletsStandalone_%s",(const char*)calibTask->GetName()));
1002
1003   TH2I *ch2d = (TH2I *) listcalibTask->FindObject("CH2d");
1004   TProfile2D *ph2d = (TProfile2D *) listcalibTask->FindObject("PH2d");
1005   TProfile2D *prf2d = (TProfile2D *) listcalibTask->FindObject("PRF2d");
1006
1007   TH2I *ch2dSum = (TH2I *) listcalibTask->FindObject(Form("CH2dSum_%s",(const char*)calibTask->GetName()));
1008   TProfile2D *ph2dSum = (TProfile2D *) listcalibTask->FindObject(Form("PH2dSum_%s",(const char*)calibTask->GetName()));
1009
1010   TH2I *ch2dSM = (TH2I *) listcalibTask->FindObject(Form("CH2dSM_%s",(const char*)calibTask->GetName()));
1011   TProfile2D *ph2dSM = (TProfile2D *) listcalibTask->FindObject(Form("PH2dSM_%s",(const char*)calibTask->GetName()));
1012
1013   AliTRDCalibraVdriftLinearFit *linearfit = (AliTRDCalibraVdriftLinearFit *) listcalibTask->FindObject("AliTRDCalibraVdriftLinearFit");
1014   AliTRDCalibraExbAltFit *exbaltfit = (AliTRDCalibraExbAltFit *) listcalibTask->FindObject("AliTRDCalibraExbAltFit");
1015   AliTRDCalibraVector *calibraVector = (AliTRDCalibraVector *) listcalibTask->FindObject("AliTRDCalibraVector");
1016
1017   //
1018
1019   THnSparseI *inhistoEntries = (THnSparseI *) fListHist->FindObject("NumberOfEntries");
1020
1021   TH1I *inEventsInput  = (TH1I *) fListHist->FindObject(Form("NEventsInput_%s",(const char*)fName));
1022   TH1I *inEvents  = (TH1I *) fListHist->FindObject(Form("NEvents_%s",(const char*)fName));
1023   TH2F *iabsoluteGain  = (TH2F *) fListHist->FindObject(Form("AbsoluteGain_%s",(const char*)fName));
1024
1025   TH1F *itrdTrack = (TH1F *) fListHist->FindObject(Form("TRDTrack_%s",(const char*)fName));
1026   TH1F *itrdTrackOffline = (TH1F *) fListHist->FindObject(Form("TRDTrackOffline_%s",(const char*)fName));
1027   TH1F *itrdTrackStandalone = (TH1F *) fListHist->FindObject(Form("TRDTrackStandalone_%s",(const char*)fName));
1028
1029   TH2F *itpctrdTrack = (TH2F *) fListHist->FindObject(Form("NbTPCTRDtrack_%s",(const char*)fName));
1030
1031   TH1F *inbTimeBin = (TH1F *) fListHist->FindObject(Form("NbTimeBin_%s",(const char*)fName));
1032   TH1F *inbTimeBinOffline = (TH1F *) fListHist->FindObject(Form("NbTimeBinOffline_%s",(const char*)fName));
1033   TH1F *inbTimeBinStandalone = (TH1F *) fListHist->FindObject(Form("NbTimeBinStandalone_%s",(const char*)fName));
1034
1035   TH1F *inbClusters = (TH1F *) fListHist->FindObject(Form("NbClusters_%s",(const char*)fName));
1036   TH1F *inbClustersOffline = (TH1F *) fListHist->FindObject(Form("NbClustersOffline_%s",(const char*)fName));
1037   TH1F *inbClustersStandalone = (TH1F *) fListHist->FindObject(Form("NbClustersStandalone_%s",(const char*)fName));
1038
1039   TH1F *inbTracklets = (TH1F *) fListHist->FindObject(Form("NbTracklets_%s",(const char*)fName));
1040   TH1F *inbTrackletsOffline = (TH1F *) fListHist->FindObject(Form("NbTrackletsOffline_%s",(const char*)fName));
1041   TH1F *inbTrackletsStandalone = (TH1F *) fListHist->FindObject(Form("NbTrackletsStandalone_%s",(const char*)fName));
1042
1043   TH2I *ich2d = (TH2I *) fListHist->FindObject("CH2d");
1044   TProfile2D *iph2d = (TProfile2D *) fListHist->FindObject("PH2d");
1045   TProfile2D *iprf2d = (TProfile2D *) fListHist->FindObject("PRF2d");
1046
1047   TH2I *ich2dSum = (TH2I *) fListHist->FindObject(Form("CH2dSum_%s",(const char*)fName));
1048   TProfile2D *iph2dSum = (TProfile2D *) fListHist->FindObject(Form("PH2dSum_%s",(const char*)fName));
1049
1050   TH2I *ich2dSM = (TH2I *) fListHist->FindObject(Form("CH2dSM_%s",(const char*)fName));
1051   TProfile2D *iph2dSM = (TProfile2D *) fListHist->FindObject(Form("PH2dSM_%s",(const char*)fName));
1052
1053   AliTRDCalibraVdriftLinearFit *ilinearfit = (AliTRDCalibraVdriftLinearFit *) fListHist->FindObject("AliTRDCalibraVdriftLinearFit");
1054   AliTRDCalibraExbAltFit *iexbaltfit = (AliTRDCalibraExbAltFit *) fListHist->FindObject("AliTRDCalibraExbAltFit");
1055   AliTRDCalibraVector *icalibraVector = (AliTRDCalibraVector *) fListHist->FindObject("AliTRDCalibraVector");
1056
1057
1058   // Add
1059
1060   if(histoEntries) {
1061     if(inhistoEntries) {
1062       inhistoEntries->Add(histoEntries);
1063       //printf("Add Events\n");
1064     }
1065     else {
1066       //printf("Create new Events\n");
1067       inhistoEntries = (THnSparseI *) histoEntries->Clone();
1068       fListHist->Add(inhistoEntries);
1069     }
1070   }
1071
1072   if(nEventsInput) {
1073     if(inEventsInput) {
1074       inEventsInput->Add(nEventsInput);
1075       //printf("Add Events\n");
1076     }
1077     else {
1078       //printf("Create new Events\n");
1079       inEventsInput = new TH1I(*nEventsInput);
1080       fListHist->Add(inEventsInput);
1081     }
1082   }
1083
1084   if(nEvents) {
1085     if(inEvents) {
1086       inEvents->Add(nEvents);
1087       //printf("Add Events\n");
1088     }
1089     else {
1090       //printf("Create new Events\n");
1091       inEvents = new TH1I(*nEvents);
1092       fListHist->Add(inEvents);
1093     }
1094   }
1095
1096   if(absoluteGain) {
1097     if(iabsoluteGain) iabsoluteGain->Add(absoluteGain);
1098     else {
1099       iabsoluteGain = new TH2F(*absoluteGain);
1100       fListHist->Add(iabsoluteGain);
1101     }
1102   }
1103
1104   if(trdTrack) {
1105     if(itrdTrack) itrdTrack->Add(trdTrack);
1106     else {
1107      itrdTrack = new TH1F(*trdTrack);
1108      fListHist->Add(itrdTrack);
1109     }
1110   }
1111
1112   if(trdTrackOffline) {
1113     if(itrdTrackOffline) itrdTrackOffline->Add(trdTrackOffline);
1114     else {
1115       itrdTrackOffline = new TH1F(*trdTrackOffline);
1116       fListHist->Add(itrdTrackOffline);
1117     }
1118   }
1119
1120   if(trdTrackStandalone) {
1121     if(itrdTrackStandalone) itrdTrackStandalone->Add(trdTrackStandalone);
1122     else {
1123       itrdTrackStandalone = new TH1F(*trdTrackStandalone);
1124       fListHist->Add(itrdTrackStandalone);
1125     }
1126   }
1127
1128   if(tpctrdTrack) {
1129     if(itpctrdTrack) itpctrdTrack->Add(tpctrdTrack);
1130     else {
1131       itpctrdTrack = new TH2F(*tpctrdTrack);
1132       fListHist->Add(itpctrdTrack);
1133     }
1134   }
1135
1136   if(nbTimeBin) {
1137     if(inbTimeBin) inbTimeBin->Add(nbTimeBin);
1138     else {
1139       inbTimeBin = new TH1F(*inbTimeBin);
1140       fListHist->Add(inbTimeBin);
1141     }
1142   }
1143
1144   if(nbTimeBinOffline) {
1145     if(inbTimeBinOffline) inbTimeBinOffline->Add(nbTimeBinOffline);
1146     else {
1147       inbTimeBinOffline = new TH1F(*nbTimeBinOffline);
1148       fListHist->Add(inbTimeBinOffline);
1149     }
1150   }
1151
1152   if(nbTimeBinStandalone) {
1153     if(inbTimeBinStandalone) inbTimeBinStandalone->Add(nbTimeBinStandalone);
1154     else {
1155       inbTimeBinStandalone = new TH1F(*nbTimeBinStandalone);
1156       fListHist->Add(inbTimeBinStandalone);
1157     }
1158   }
1159
1160   if(nbClusters) {
1161     if(inbClusters) inbClusters->Add(nbClusters);
1162     else {
1163       inbClusters = new TH1F(*nbClusters);
1164       fListHist->Add(inbClusters);
1165     }
1166   }
1167
1168   if(nbClustersOffline) {
1169     if(inbClustersOffline) inbClustersOffline->Add(nbClustersOffline);
1170     else {
1171       inbClustersOffline = new TH1F(*nbClustersOffline);
1172       fListHist->Add(inbClustersOffline);
1173     }
1174   }
1175
1176   if(nbClustersStandalone) {
1177     if(inbClustersStandalone) inbClustersStandalone->Add(nbClustersStandalone);
1178     else {
1179       inbClustersStandalone = new TH1F(*nbClustersStandalone);
1180       fListHist->Add(inbClustersStandalone);
1181     }
1182   }
1183
1184   if(nbTracklets) {
1185     if(inbTracklets) inbTracklets->Add(nbTracklets);
1186     else {
1187       inbTracklets = new TH1F(*nbTracklets);
1188       fListHist->Add(inbTracklets);
1189     }
1190   }
1191
1192   if(nbTrackletsOffline) {
1193     if(inbTrackletsOffline) inbTrackletsOffline->Add(nbTrackletsOffline);
1194     else {
1195       inbTrackletsOffline = new TH1F(*nbTrackletsOffline);
1196       fListHist->Add(inbTrackletsOffline);
1197     }
1198   }
1199
1200   if(nbTrackletsStandalone) {
1201     if(inbTrackletsStandalone) inbTrackletsStandalone->Add(nbTrackletsStandalone);
1202     else {
1203       inbTrackletsStandalone = new TH1F(*nbTrackletsStandalone);
1204       fListHist->Add(inbTrackletsStandalone);
1205     }
1206   }
1207
1208   if(ch2d) {
1209     if(ich2d) ich2d->Add(ch2d);
1210     else {
1211       ich2d = new TH2I(*ch2d);
1212       fListHist->Add(ich2d);
1213     }
1214   }
1215
1216   if(ph2d) {
1217     if(iph2d) iph2d->Add(ph2d);
1218     else {
1219       iph2d = new TProfile2D(*ph2d);
1220       fListHist->Add(iph2d);
1221     }
1222   }
1223
1224   if(prf2d) {
1225     if(iprf2d) iprf2d->Add(prf2d);
1226     else {
1227       iprf2d = new TProfile2D(*prf2d);
1228       fListHist->Add(iprf2d);
1229     }
1230   }
1231
1232   if(ch2dSum) {
1233     if(ich2dSum) ich2dSum->Add(ch2dSum);
1234     else {
1235       ich2dSum = new TH2I(*ch2dSum);
1236       fListHist->Add(ich2dSum);
1237     }
1238   }
1239
1240   if(ph2dSum) {
1241     if(iph2dSum) iph2dSum->Add(ph2dSum);
1242     else {
1243       iph2dSum = new TProfile2D(*ph2dSum);
1244       fListHist->Add(iph2dSum);
1245     }
1246   }
1247
1248   if(ch2dSM) {
1249     if(ich2dSM) ich2dSM->Add(ch2dSM);
1250     else {
1251       ich2dSM = new TH2I(*ch2dSM);
1252       fListHist->Add(ich2dSM);
1253     }
1254   }
1255
1256   if(ph2dSM) {
1257     if(iph2dSM) iph2dSM->Add(ph2dSM);
1258     else {
1259       iph2dSM = new TProfile2D(*ph2dSM);
1260       fListHist->Add(iph2dSM);
1261     }
1262   }
1263
1264   if(linearfit) {
1265     if(ilinearfit) ilinearfit->Add(linearfit);
1266     else {
1267       ilinearfit = new AliTRDCalibraVdriftLinearFit(*linearfit);
1268       fListHist->Add(ilinearfit);
1269     }
1270   }
1271
1272   if(exbaltfit) {
1273     if(iexbaltfit) iexbaltfit->Add(exbaltfit);
1274     else {
1275       iexbaltfit = new AliTRDCalibraExbAltFit(*exbaltfit);
1276       fListHist->Add(iexbaltfit);
1277     }
1278   }
1279
1280   if(calibraVector) {
1281     if(icalibraVector) icalibraVector->Add(calibraVector);
1282     else {
1283       icalibraVector = new AliTRDCalibraVector(*calibraVector);
1284       fListHist->Add(icalibraVector);
1285     }
1286   }
1287
1288 }
1289 //________________________________________________________________________________
1290 Long64_t AliTRDCalibTask::Merge(TCollection *li) {
1291
1292   //
1293   // merge component
1294   //
1295
1296   TIterator* iter = li->MakeIterator();
1297   AliTRDCalibTask* cal = 0;
1298
1299   while ((cal = (AliTRDCalibTask*)iter->Next())) {
1300     if (!cal->InheritsFrom(AliTRDCalibTask::Class())) {
1301       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
1302       return -1;
1303     }
1304
1305     // add histograms here...
1306     this->AddTask(cal);
1307
1308   }
1309
1310   return 0;
1311
1312 }
1313 //_____________________________________________________
1314 Bool_t AliTRDCalibTask::SetVersionSubversion(){
1315   //
1316   // Load Chamber Gain factors into the Tender supply
1317   //
1318
1319   printf("SetVersionSubversion\n");
1320
1321   //find previous entry from the UserInfo
1322   TTree *tree=((TChain*)GetInputData(0))->GetTree();
1323   if (!tree) {
1324     AliError("Tree not found in ESDhandler");
1325     return kFALSE;
1326   }
1327
1328   TList *userInfo=(TList*)tree->GetUserInfo();
1329   if (!userInfo) {
1330     AliError("No UserInfo found in tree");
1331     return kFALSE;
1332   }
1333
1334   TList *cdbList=(TList*)userInfo->FindObject("cdbList");
1335   if (!cdbList) {
1336     AliError("No cdbList found in UserInfo");
1337     if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print();
1338     return kFALSE;
1339   }
1340
1341   TIter nextCDB(cdbList);
1342   TObjString *os=0x0;
1343   while ( (os=(TObjString*)nextCDB()) ){
1344     if(os->GetString().Contains("TRD/Calib/ChamberGainFactor")){
1345       // Get Old gain calibration
1346       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1347       fFirstRunGain = id->GetFirstRun();
1348       fVersionGainUsed = id->GetVersion();
1349       fSubVersionGainUsed = id->GetSubVersion();
1350     } else if(os->GetString().Contains("TRD/Calib/ChamberVdrift")){
1351       // Get Old drift velocity calibration
1352       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1353       fFirstRunVdrift = id->GetFirstRun();
1354       fVersionVdriftUsed = id->GetVersion();
1355       fSubVersionVdriftUsed = id->GetSubVersion();
1356     } else if(os->GetString().Contains("TRD/Calib/LocalGainFactor")){
1357       // Get Old drift velocity calibration
1358       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1359       fFirstRunGainLocal = id->GetFirstRun();
1360       fVersionGainLocalUsed = id->GetVersion();
1361       fSubVersionGainLocalUsed = id->GetSubVersion();
1362     } else if(os->GetString().Contains("TRD/Calib/ChamberExB")){
1363       // Get Old drift velocity calibration
1364       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1365       fFirstRunExB = id->GetFirstRun();
1366       fVersionExBUsed = id->GetVersion();
1367       fSubVersionExBUsed = id->GetSubVersion();
1368     }
1369   }
1370
1371   //printf("VersionGain %d, SubversionGain %d, VersionLocalGain %d, Subversionlocalgain %d, Versionvdrift %d, Subversionvdrift %d\n",fVersionGainUsed,fSubVersionGainUsed,fVersionGainLocalUsed,fSubVersionGainLocalUsed,fVersionVdriftUsed,fSubVersionVdriftUsed);
1372
1373   // Check
1374   if((fFirstRunGain < 0)            ||
1375      (fFirstRunGainLocal < 0)       ||
1376      (fFirstRunVdrift < 0)          ||
1377      (fVersionGainUsed < 0)         ||
1378      (fVersionGainLocalUsed < 0)    ||
1379      (fSubVersionGainUsed < 0)      ||
1380      (fSubVersionGainLocalUsed < 0) ||
1381      (fVersionVdriftUsed < 0)       ||
1382      (fSubVersionVdriftUsed < 0)) {
1383     AliError("No recent calibration found");
1384     return kFALSE;
1385   }
1386   else return kTRUE;
1387
1388 }
1389 //_________________________________________________________________________________________________________________________
1390 Bool_t AliTRDCalibTask::ParticleGood(int i) const {
1391
1392   //
1393   // Definition of good tracks
1394   //
1395
1396
1397   AliESDtrack *track = fESD->GetTrack(i);
1398   if (!track->IsOn(AliESDtrack::kTPCrefit)) return 0;        // TPC refit
1399   if (track->GetTPCNcls() < 90) return 0;                    // number of TPC clusters
1400   if (fabs(track->Eta())>0.8) return 0;                         // fiducial pseudorapidity
1401   Float_t r,z;
1402   track->GetImpactParametersTPC(r,z);
1403   if (fabs(z)>2.0) return 0;                          // impact parameter in z
1404   if (fabs(r)>2.0) return 0;                          // impact parameter in xy
1405   if (r==0) return 0;
1406   return 1;
1407
1408
1409 }
1410 //______________________________________________________________________________________________________________________
1411 Bool_t AliTRDCalibTask::FindP1TrackPHtrackletV1Test(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1412 {
1413   //
1414   // Drift velocity calibration:
1415   // Fit the clusters with a straight line
1416   // From the slope find the drift velocity
1417   //
1418
1419   ////////////////////////////////////////////////
1420   //Number of points: if less than 3 return kFALSE
1421   /////////////////////////////////////////////////
1422   if(nbclusters <= 2) return kFALSE;
1423
1424   ////////////
1425   //Variables
1426   ////////////
1427   // results of the linear fit
1428   Double_t dydt                       = 0.0;                                // dydt tracklet after straight line fit
1429   Double_t errorpar                   = 0.0;                                // error after straight line fit on dy/dt
1430   Double_t pointError                 = 0.0;                                // error after straight line fit
1431   // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
1432   Int_t    crossrow                   = 0;                                  // if it crosses a pad row
1433   Int_t    rowp                       = -1;                                 // if it crosses a pad row
1434   Float_t  tnt                        = tracklet->GetTilt();                // tan tiltingangle
1435   TLinearFitter linearFitterTracklet(2,"pol1");
1436   linearFitterTracklet.StoreData(kTRUE);
1437
1438
1439   ///////////////////////////////////////////
1440   // Take the parameters of the track
1441   //////////////////////////////////////////
1442   // take now the snp, tnp and tgl from the track
1443   Double_t snp = tracklet->GetSnp();             // sin dy/dx at the end of the chamber
1444   Double_t tnp = 0.0;                            // dy/dx at the end of the chamber
1445   if( TMath::Abs(snp) <  1.){
1446     tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1447   }
1448   Double_t tgl  = tracklet->GetTgl();           // dz/dl
1449   Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp);   // dz/dx calculated from dz/dl
1450   // at the entrance
1451   //Double_t tnp = tracklet->GetYref(1);      // dy/dx at the entrance of the chamber
1452   //Double_t tgl = tracklet->GetZref(1);      // dz/dl at the entrance of the chamber
1453   //Double_t dzdx = tgl;                      //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1454   // at the end with correction due to linear fit
1455   //Double_t tnp = tracklet->GetYfit(1);      // dy/dx at the end of the chamber after fit correction
1456   //Double_t tgl = tracklet->GetZfit(1);      // dz/dl at the end of the chamber after fit correction
1457
1458
1459   ////////////////////////////
1460   // loop over the clusters
1461   ////////////////////////////
1462   Int_t  nbli = 0;
1463   AliTRDcluster *cl                   = 0x0;
1464   //////////////////////////////
1465   // Check no shared clusters
1466   //////////////////////////////
1467   for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
1468     cl = tracklet->GetClusters(icc);
1469     if(cl)  crossrow = 1;
1470   }
1471   //////////////////////////////////
1472   // Loop clusters
1473   //////////////////////////////////
1474   for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1475     if(!(cl = tracklet->GetClusters(ic))) continue;
1476     //if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
1477
1478     Double_t ycluster                 = cl->GetY();
1479     Int_t time                        = cl->GetPadTime();
1480     Double_t timeis                   = time/10.0;
1481     //See if cross two pad rows
1482     Int_t    row                      = cl->GetPadRow();
1483     if(rowp==-1) rowp                 = row;
1484     if(row != rowp) crossrow          = 1;
1485
1486     linearFitterTracklet.AddPoint(&timeis,ycluster,1);
1487     nbli++;
1488
1489
1490   }
1491
1492   ////////////////////////////////////
1493   // Do the straight line fit now
1494   ///////////////////////////////////
1495   if(nbli <= 2){
1496     linearFitterTracklet.ClearPoints();
1497     return kFALSE;
1498   }
1499   TVectorD pars;
1500   linearFitterTracklet.Eval();
1501   linearFitterTracklet.GetParameters(pars);
1502   pointError  =  TMath::Sqrt(linearFitterTracklet.GetChisquare()/(nbli-2));
1503   errorpar    =  linearFitterTracklet.GetParError(1)*pointError;
1504   dydt        = pars[1];
1505   //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",linearFitterTracklet->GetChisquare(),nbli,pointError,linearFitterTracklet->GetParError(1),errorpar);
1506   linearFitterTracklet.ClearPoints();
1507
1508   /////////////////////////
1509   // Cuts quality
1510   ////////////////////////
1511
1512   if(nbclusters < fLow) return kFALSE;
1513   if(nbclusters > fHigh) return kFALSE;
1514   if(pointError >= 0.3) return kFALSE;
1515   if(crossrow == 1) return kTRUE;
1516
1517   ///////////////////////
1518   // Fill
1519   //////////////////////
1520
1521   if(fDebug > 0){
1522     //Add to the linear fitter of the detector
1523     if( TMath::Abs(snp) <  1.){
1524       Double_t x = tnp-dzdx*tnt;
1525       //if(!fLinearVdriftTest) printf("Not there\n");
1526       Double_t nbentries = fLinearVdriftTest->GetEntries();
1527       if(nbentries < (5.0*32767)) fLinearVdriftTest->Fill(x,dydt);
1528     }
1529   }
1530
1531   return kTRUE;
1532 }
1533