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