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   PostData(1,fListHist);
 412
 413   //cout << "AliTRDCalibTask::UserCreateOutputObjects() OUT" << endl;
 414
 415 }
 416
 417 //________________________________________________________________________
 418 void AliTRDCalibTask::UserExec(Option_t *)
 419 {
 420   //
 421   // Filling of the histos
 422   //
 423   //cout << "AliTRDCalibTask::Exec() IN" << endl;
 424
 425   // Init Versions and subversions used
 426   if((fFirstRunGain==-1) || (fVersionGainUsed==-1) || (fSubVersionGainUsed==-1) || (fFirstRunGainLocal==-1) || (fVersionGainLocalUsed==-1) || (fSubVersionGainLocalUsed==-1) || (fFirstRunVdrift==-1) || (fVersionVdriftUsed==-1) || (fSubVersionVdriftUsed==-1)) {
 427     if(!SetVersionSubversion()) {
 428       PostData(1, fListHist);
 429       return;
 430     }
 431   }
 432   if(fCounter==0) {
 433     fTRDCalibraFillHisto->SetFirstRunGain(fFirstRunGain); // Gain Used
 434     fTRDCalibraFillHisto->SetVersionGainUsed(fVersionGainUsed); // Gain Used
 435     fTRDCalibraFillHisto->SetSubVersionGainUsed(fSubVersionGainUsed); // Gain Used
 436     fTRDCalibraFillHisto->SetFirstRunGainLocal(fFirstRunGainLocal); // Gain Used
 437     fTRDCalibraFillHisto->SetVersionGainLocalUsed(fVersionGainLocalUsed); // Gain Used
 438     fTRDCalibraFillHisto->SetSubVersionGainLocalUsed(fSubVersionGainLocalUsed); // Gain Used
 439     fTRDCalibraFillHisto->SetFirstRunVdrift(fFirstRunVdrift); // Vdrift Used
 440     fTRDCalibraFillHisto->SetVersionVdriftUsed(fVersionVdriftUsed); // Vdrift Used
 441     fTRDCalibraFillHisto->SetSubVersionVdriftUsed(fSubVersionVdriftUsed); // Vdrift Used
 442     fTRDCalibraFillHisto->InitCalDet();
 443   }
 444
 445   //  AliLog::SetGlobalLogLevel(AliLog::kError);
 446   //  cout << "AliTRDCalibTask::Exec() 1" << endl;
 447   fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
 448   if(!fESD){
 449     AliError("ESD Event missing");
 450     PostData(1, fListHist);
 451     return;
 452   }
 453
 454   const char* type = fESD->GetBeamType();
 455
 456
 457   //printf("Counter %d\n",fCounter);
 458
 459   fCounter++;
 460   fNEventsInput->Fill(1);
 461
 462   //cout << "maxEvent = " << fMaxEvent << endl;
 463   //if(fCounter%100==0) cout << "fCounter = " << fCounter << endl;
 464   if((fMaxEvent != 0) && (fMaxEvent < fCounter)) {
 465     PostData(1, fListHist);
 466     return;
 467   }
 468   //if(fCounter%100==0) cout << "fCounter1 = " << fCounter << endl;
 469   //cout << "event = " << fCounter << endl;
 470
 471   //printf("Counter %d\n",fCounter);
 472
 473   ///////////////////
 474   // Check trigger
 475   ///////////////////
 476   Bool_t pass = kTRUE;
 477
 478   if (strstr(type,"p-p")) {
 479
 480     //printf("Will check the triggers\n");
 481
 482     Int_t numberOfTriggerSelected = fSelectedTrigger->GetEntriesFast();
 483     //printf("numberofTriggerSelected %d\n",numberOfTriggerSelected);
 484     if(fRejected) {
 485       pass = kTRUE;
 486       for(Int_t k = 0; k < numberOfTriggerSelected; k++){
 487         const TObjString *const obString=(TObjString*)fSelectedTrigger->At(k);
 488         const TString tString=obString->GetString();
 489         if(fESD->IsTriggerClassFired((const char*)tString)) {
 490           pass = kFALSE;
 491         }
 492       }
 493     }
 494     else {
 495       pass = kFALSE;
 496       for(Int_t k = 0; k < numberOfTriggerSelected; k++){
 497         const TObjString *const obString=(TObjString*)fSelectedTrigger->At(k);
 498         const TString tString=obString->GetString();
 499         if(fESD->IsTriggerClassFired((const char*)tString)) {
 500           pass = kTRUE;
 501         }
 502       }
 503     }
 504     if(!pass) {
 505       PostData(1, fListHist);
 506       return;
 507     }
 508
 509   }
 510
 511   //printf("Class Fired %s\n",(const char*)fESD->GetFiredTriggerClasses());
 512   //printf("Trigger passed\n");
 513
 514   ///////////////////////////////
 515   // Require a primary vertex
 516   //////////////////////////////
 517   if(fRequirePrimaryVertex) {
 518     const AliESDVertex* vtxESD = 0x0;
 519     if      (fVtxTPC) vtxESD = fESD->GetPrimaryVertexTPC() ;
 520     else if (fVtxSPD) vtxESD = fESD->GetPrimaryVertexSPD() ;
 521     else              vtxESD = fESD->GetPrimaryVertexTracks() ;
 522     if(!vtxESD){
 523       PostData(1, fListHist);
 524       return;
 525     }
 526     Int_t nCtrb = vtxESD->GetNContributors();
 527     if(nCtrb < fMinNbContributors) {
 528       PostData(1, fListHist);
 529       return;
 530     }
 531     Double_t zPosition = vtxESD->GetZ();
 532     if(TMath::Abs(zPosition) > fRangePrimaryVertexZ) {
 533       PostData(1, fListHist);
 534       return;
 535     }
 536
 537   }
 538
 539   //printf("Primary vertex passed\n");
 540
 541   //////////////////////////////////////
 542   // Requirement on number of good tracks
 543   //////////////////////////////////////
 544   Int_t nGoodParticles = 0;
 545   Double_t nbTracks = fESD->GetNumberOfTracks();
 546   for(Int_t itrack = 0; itrack < nbTracks; itrack++) {
 547     if(ParticleGood(itrack)) nGoodParticles++;
 548   }
 549   if(fDebug > 0)  {
 550     // Centrality
 551     AliCentrality *esdCentrality = fESD->GetCentrality();
 552     Float_t centrality = esdCentrality->GetCentralityPercentile("V0M");
 553     //Float_t centralityb = esdCentrality->GetCentralityPercentile("CL1");
 554     fNbGoodTracks->Fill(nGoodParticles,centrality);
 555     //printf("centrality %f, centralityb %f\n",centrality,centralityb);
 556   }
 557
 558   if (strstr(type,"Pb-Pb")) {
 559     //printf("Will check the number of good tracks\n");
 560     if((nGoodParticles < fMinNbTracks) || (nGoodParticles > fMaxNbTracks)) {
 561       PostData(1, fListHist);
 562       return;
 563     }
 564   }
 565
 566   fNEvents->Fill(1);
 567
 568   // In total
 569   Int_t nbTrdTracks = 0;
 570   // standalone
 571   Int_t nbTrdTracksStandalone = 0;
 572   // offline
 573   Int_t nbTrdTracksOffline = 0;
 574   // TPC
 575   Int_t nbtrackTPC = 0;
 576
 577
 578
 579   if (nbTracks <= 0.0) {
 580
 581     if(fDebug > 1) {
 582       fNbTRDTrack->Fill(nbTrdTracks);
 583       fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
 584       fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
 585     }
 586     PostData(1, fListHist);
 587     return;
 588   }
 589
 590
 591   fESDfriend = dynamic_cast<AliESDfriend*> (fESD->FindListObject("AliESDfriend"));
 592   if(!fESDfriend){
 593     AliError("fESDfriend not available");
 594     PostData(1, fListHist);
 595     return;
 596   }
 597
 598   //printf("has friends\n");
 599
 600   /*
 601   ////////////////////////////////////
 602    // Check the number of TPC tracks
 603    ///////////////////////////////////
 604    //printf("Nb of tracks %f\n",nbTracks);
 605    for(Int_t itrk = 0; itrk < nbTracks; itrk++){
 606      // Get ESD track
 607      fkEsdTrack = fESD->GetTrack(itrk);
 608      ULong_t status = fkEsdTrack->GetStatus();
 609      if(status&(AliESDtrack::kTPCout)) nbtrackTPC++;
 610      if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {
 611        nbTrdTracks++;
 612        nbTrdTracksStandalone++;
 613      }
 614      if((status&(AliESDtrack::kTRDin))) {
 615        nbTrdTracks++;
 616        nbTrdTracksOffline++;
 617      }
 618    }
 619
 620    if((nbtrackTPC>0) && (nbTrdTracks > (3.0*nbtrackTPC))) pass = kFALSE;
 621
 622    if(fDebug > 1) {
 623
 624      fNbTRDTrack->Fill(nbTrdTracks);
 625      fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
 626      fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
 627      fNbTPCTRDtrack->Fill(nbTrdTracks,nbtrackTPC);
 628
 629    }
 630
 631    if(!pass) {
 632      PostData(1, fListHist);
 633      return;
 634    }
 635   */
 636
 637   /////////////////////////////////////
 638   // Loop on AliESDtrack
 639   ////////////////////////////////////
 640   //printf("Nb of tracks %f\n",nbTracks);
 641   for(int itrk=0; itrk < nbTracks; ++itrk){
 642
 643     // Get ESD track
 644     fkEsdTrack = fESD->GetTrack(itrk);
 645     if(!fkEsdTrack) continue;
 646     ULong_t status = fkEsdTrack->GetStatus();
 647     if(status&(AliESDtrack::kTPCout)) ++nbtrackTPC;
 648
 649     // Quality cuts on the AliESDtrack
 650     if((fEsdTrackCuts) && (!fEsdTrackCuts->IsSelected((AliVParticle *)fkEsdTrack))) {
 651       //printf("Not a good track\n");
 652       continue;
 653     }
 654
 655     // First Absolute gain calibration
 656     Int_t trdNTracklets = (Int_t) fkEsdTrack->GetTRDntracklets();
 657     Int_t trdNTrackletsPID = (Int_t) fkEsdTrack->GetTRDntrackletsPID();
 658     if((trdNTracklets > 0) && (trdNTrackletsPID > 0)) {
 659       for(Int_t iPlane = 0; iPlane < 6; ++iPlane){
 660         //Double_t slide = fkEsdTrack->GetTRDslice(iPlane);
 661         //printf("Number of slide %d\n",fkEsdTrack->GetNumberOfTRDslices());
 662         //Double_t momentum = fkEsdTrack->GetTRDmomentum(iPlane);
 663         //printf("momentum %f, slide %f\n",momentum,slide);
 664         if(fkEsdTrack->GetTRDslice(iPlane) > 0.0)
 665           fAbsoluteGain->Fill(fkEsdTrack->GetTRDslice(iPlane)*8.0/100.0,
 666                               fkEsdTrack->GetTRDmomentum(iPlane));
 667       }
 668     }
 669
 670     // Other cuts
 671     Bool_t good = kTRUE;
 672     Bool_t standalonetrack = kFALSE;
 673     Bool_t offlinetrack = kFALSE;
 674     //ULong_t status = fkEsdTrack->GetStatus();
 675
 676     fFriendTrack = fESDfriend->GetTrack(itrk);
 677     if(!fFriendTrack)  {
 678       //printf("No friend track %d\n",itrk);
 679       continue;
 680     }
 681     //////////////////////////////////////
 682     // Loop on calibration objects
 683     //////////////////////////////////////
 684     Int_t icalib=0;
 685     Int_t nTRDtrackV1=0;
 686     while((fCalibObject = (TObject *)(fFriendTrack->GetCalibObject(icalib++)))){
 687       //printf("Name %s\n",fCalibObject->IsA()->GetName());
 688       if(strcmp(fCalibObject->IsA()->GetName(), "AliTRDtrackV1") != 0) continue;
 689       //printf("Find the calibration object\n");
 690       ++nTRDtrackV1;
 691
 692       if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {
 693         standalonetrack = kTRUE;
 694       }
 695       if((status&(AliESDtrack::kTRDin))) {
 696         offlinetrack = kTRUE;
 697       }
 698       if(fOfflineTracks){
 699         if(!offlinetrack){
 700           good = kFALSE;
 701         }
 702       }
 703       else if(fStandaloneTracks){
 704         if(!standalonetrack){
 705           good = kFALSE;
 706         }
 707       }
 708
 709       fTrdTrack = (AliTRDtrackV1 *)fCalibObject;
 710       if(good) {
 711         //cout << "good" << endl;
 712         fTRDCalibraFillHisto->UpdateHistogramsV1(fTrdTrack);
 713         //printf("Fill fTRDCalibraFillHisto\n");
 714       }
 715
 716       //////////////////////////////////
 717       // Debug
 718       ////////////////////////////////
 719
 720       if(fDebug > 0) {
 721
 722         //printf("Enter debug\n");
 723
 724         Int_t nbtracklets = 0;
 725
 726         // Check some stuff
 727         Bool_t standalonetracklet = kFALSE;
 728         const AliTRDseedV1 *tracklet = 0x0;
 729         //////////////////////////////////////
 730         // Loop tracklets
 731         /////////////////////////////////////
 732         Int_t nbclusters=0;
 733         Double_t phtb[AliTRDseedV1::kNtb];
 734         memset(phtb, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
 735         Double_t sum = 0.0;
 736         Float_t normalisation = 6.67;
 737         Int_t detector = 0;
 738         Int_t sector = 0;
 739         for(Int_t itr = 0; itr < 6; ++itr){
 740
 741           if(!(tracklet = fTrdTrack->GetTracklet(itr))) continue;
 742           if(!tracklet->IsOK()) continue;
 743           ++nbtracklets;
 744           standalonetracklet = kFALSE;
 745           if(tracklet->IsStandAlone()) standalonetracklet = kTRUE;
 746
 747           nbclusters = 0;
 748           memset(phtb, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
 749           sum = 0.0;
 750           detector = 0;
 751           sector = 0;
 752           //Int_t crossrow = 0;
 753
 754           // Check no shared clusters
 755           //for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
 756           //  if((fcl = tracklet->GetClusters(icc)))  crossrow = 1;
 757           // }
 758
 759           // Loop on clusters
 760           Int_t time = 0;
 761           Float_t ch = 0;
 762           Float_t qcl = 0;
 763           for(int ic=0; ic<AliTRDseedV1::kNtb; ++ic){
 764
 765             if(!(fCl = tracklet->GetClusters(ic))) continue;
 766             ++nbclusters;
 767             time = fCl->GetPadTime();
 768             ch =  tracklet->GetdQdl(ic);
 769             qcl = TMath::Abs(fCl->GetQ());
 770             detector = fCl->GetDetector();
 771             // Add the charge if shared cluster
 772             if((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) {
 773               if((fCl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) {
 774                 qcl += TMath::Abs(fCl->GetQ());
 775                 //printf("Add the cluster charge\n");
 776               }
 777             }
 778             if((time>-1) && (time<fNbTimeBins)) phtb[time]=qcl;
 779             if((fCalDetGain) && (fCalDetGain->GetValue(detector) > 0.0)) sum += ch*fCalDetGain->GetValue(detector)/normalisation;
 780             else sum += ch/normalisation;
 781
 782             if(fDebug > 1) {
 783               fNbTimeBin->Fill(time);
 784               if(tracklet->IsStandAlone()) fNbTimeBinStandalone->Fill(time);
 785               else fNbTimeBinOffline->Fill(time);
 786             }
 787           }
 788           sector = AliTRDgeometry::GetSector(detector);
 789
 790           if(fDebug > 1) {
 791             fNbTracklets->Fill(detector);
 792             if(tracklet->IsStandAlone()) fNbTrackletsStandalone->Fill(detector);
 793             else fNbTrackletsOffline->Fill(detector);
 794
 795             fNbClusters->Fill(nbclusters);
 796             if(tracklet->IsStandAlone())  fNbClustersStandalone->Fill(nbclusters);
 797             else  fNbClustersOffline->Fill(nbclusters);
 798           }
 799
 800           if(fDebug > 0) {
 801             if((nbclusters > fLow) && (nbclusters < fHigh)){
 802               if(fRelativeScale > 0.0) sum = sum/fRelativeScale;
 803               fCH2dSM->Fill(sum,sector+0.5);
 804               fCH2dSum->Fill(sum,0.5);
 805               Bool_t checknoise = kTRUE;
 806               if(fMaxCluster > 0) {
 807                 if(phtb[0] > fMaxCluster) checknoise = kFALSE;
 808                 if(fNbTimeBins > fNbMaxCluster) {
 809                   for(Int_t k = (fNbTimeBins-fNbMaxCluster); k < fNbTimeBins; k++){
 810                     if(phtb[k] > fMaxCluster) checknoise = kFALSE;
 811                   }
 812                 }
 813               }
 814               if(checknoise) {
 815                 for(int ic=0; ic<fNbTimeBins; ic++){
 816                   if(fFillZero) {
 817                     fPH2dSum->Fill((Double_t)(ic/10.0),0.5,(Double_t)phtb[ic]);
 818                     fPH2dSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
 819                   }
 820                   else {
 821                     if(phtb[ic] > 0.0) {
 822                       fPH2dSum->Fill((Double_t)(ic/10.0),0.0,(Double_t)phtb[ic]);
 823                       fPH2dSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
 824                     }
 825                   }
 826                 }
 827               }
 828             }
 829           }
 830         } // loop on tracklets
 831
 832       } // debug
 833
 834     }// while calibration objects
 835     if(nTRDtrackV1 > 0) {
 836       ++nbTrdTracks;
 837       if((status&(AliESDtrack::kTRDout)) && (!(status&(AliESDtrack::kTRDin)))) {
 838         ++nbTrdTracksStandalone;
 839       }
 840       if((status&(AliESDtrack::kTRDin))) {
 841         ++nbTrdTracksOffline;
 842       }
 843     }
 844     //delete fFriendTrack;
 845   } // loop ESD track
 846
 847   if(fDebug > 1) {
 848     fNbTRDTrack->Fill(nbTrdTracks);
 849     fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
 850     fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
 851     fNbTPCTRDtrack->Fill(nbTrdTracks,nbtrackTPC);
 852   }
 853
 854   // Post output data
 855   PostData(1, fListHist);
 856   //cout << "AliTRDCalibTask::Exec() OUT" << endl;
 857 }
 858
 859 //________________________________________________________________________
 860 void AliTRDCalibTask::Terminate(Option_t *)
 861 {
 862   //
 863   // Terminate
 864   //
 865
 866   if(fTRDCalibraFillHisto) fTRDCalibraFillHisto->DestroyDebugStreamer();
 867
 868
 869 }
 870 //_______________________________________________________
 871 Bool_t AliTRDCalibTask::Load(const Char_t *filename)
 872 {
 873   //
 874   // Generic container loader
 875   //
 876
 877   if(!TFile::Open(filename)){
 878     //AliWarning(Form("Couldn't open file %s.", filename));
 879     return kFALSE;
 880   }
 881   TList *o = 0x0;
 882   if(!(o = (TList*)gFile->Get(GetName()))){
 883     //AliWarning("Missing histogram container.");
 884     return kFALSE;
 885   }
 886   fListHist = (TList*)o->Clone(GetName());
 887   gFile->Close();
 888   return kTRUE;
 889 }
 890 //_______________________________________________________
 891 Bool_t AliTRDCalibTask::Load(TList *lister)
 892 {
 893   //
 894   // Generic container loader
 895   //
 896
 897   fListHist = (TList*)lister->Clone(GetName());
 898   return kTRUE;
 899 }
 900 //________________________________________________________________________
 901 void AliTRDCalibTask::Plot()
 902 {
 903   //
 904   // Plot the histos stored in the TList
 905   //
 906
 907   if(!fListHist) return;
 908
 909   /////////////////////////////////////
 910   // Take the debug stuff
 911   /////////////////////////////////////
 912
 913   TH1I *nEvents  = (TH1I *) fListHist->FindObject("NEvents");
 914
 915   TH2F *absoluteGain  = (TH2F *) fListHist->FindObject("AbsoluteGain");
 916
 917   TH1F *trdTrack = (TH1F *) fListHist->FindObject("TRDTrack");
 918   TH1F *trdTrackOffline = (TH1F *) fListHist->FindObject("TRDTrackOffline");
 919   TH1F *trdTrackStandalone = (TH1F *) fListHist->FindObject("TRDTrackStandalone");
 920
 921   TH2F *tpctrdTrack = (TH2F *) fListHist->FindObject("NbTPCTRDtrack");
 922
 923   TH1F *nbTimeBin = (TH1F *) fListHist->FindObject("NbTimeBin");
 924   TH1F *nbTimeBinOffline = (TH1F *) fListHist->FindObject("NbTimeBinOffline");
 925   TH1F *nbTimeBinStandalone = (TH1F *) fListHist->FindObject("NbTimeBinStandalone");
 926
 927   TH1F *nbClusters = (TH1F *) fListHist->FindObject("NbClusters");
 928   TH1F *nbClustersOffline = (TH1F *) fListHist->FindObject("NbClustersOffline");
 929   TH1F *nbClustersStandalone = (TH1F *) fListHist->FindObject("NbClustersStandalone");
 930
 931   TH1F *nbTracklets = (TH1F *) fListHist->FindObject("NbTracklets");
 932   TH1F *nbTrackletsOffline = (TH1F *) fListHist->FindObject("NbTrackletsOffline");
 933   TH1F *nbTrackletsStandalone = (TH1F *) fListHist->FindObject("NbTrackletsStandalone");
 934
 935   /////////////////////////////////////
 936   // Take the calibration objects
 937   /////////////////////////////////////
 938
 939   TH2I *ch2d = (TH2I *) fListHist->FindObject("CH2d");
 940   TProfile2D *ph2d = (TProfile2D *) fListHist->FindObject("PH2d");
 941
 942   TH2I *ch2dSum = (TH2I *) fListHist->FindObject("CH2dSum");
 943   TProfile2D *ph2dSum = (TProfile2D *) fListHist->FindObject("PH2dSum");
 944
 945   TH2I *ch2dSM = (TH2I *) fListHist->FindObject("CH2dSM");
 946   TProfile2D *ph2dSM = (TProfile2D *) fListHist->FindObject("PH2dSM");
 947
 948   AliTRDCalibraVdriftLinearFit *linearfit = (AliTRDCalibraVdriftLinearFit *) fListHist->FindObject("AliTRDCalibraVdriftLinearFit");
 949
 950   ////////////////////////////////////////////////
 951   // Add the AliTRDCalibraVdriftLinearFit
 952   ///////////////////////////////////////////////
 953
 954   Int_t first = 0;
 955   TH2S *histolinearfitsum = 0x0;
 956
 957   if(linearfit) {
 958     for(Int_t det = 0; det < 540; det++) {
 959       if(linearfit->GetLinearFitterHisto(det)){
 960         if(TMath::Abs(first)<0.0001){
 961           histolinearfitsum = linearfit->GetLinearFitterHisto(det);
 962           first += 1;
 963         }
 964         else {
 965           if (histolinearfitsum) {
 966             histolinearfitsum->Add(linearfit->GetLinearFitterHisto(det));
 967           }
 968         }
 969       }
 970     }
 971   }
 972
 973   ///////////////////////////
 974   // Style
 975   //////////////////////////
 976
 977   gStyle->SetPalette(1);
 978   gStyle->SetOptStat(1111);
 979   gStyle->SetOptFit(1111);
 980   gStyle->SetPadBorderMode(0);
 981   gStyle->SetCanvasColor(10);
 982   gStyle->SetPadLeftMargin(0.13);
 983   gStyle->SetPadRightMargin(0.13);
 984
 985   /////////////////////////
 986   // Plot
 987   ////////////////////////
 988
 989  if(nEvents) {
 990
 991     TCanvas *debugEvents = new TCanvas("cNEvents","cNEvents",10,10,510,510);
 992     debugEvents->cd(1);
 993     if(nEvents) nEvents->Draw();
 994
 995   }
 996
 997  if(absoluteGain) {
 998
 999     TCanvas *debugAbsoluteGain = new TCanvas("cAbsoluteGain","cAbsoluteGain",10,10,510,510);
1000     debugAbsoluteGain->cd(1);
1001     if(absoluteGain) absoluteGain->Draw();
1002
1003   }
1004
1005   if(trdTrack || tpctrdTrack) {
1006
1007     TCanvas *debugtrdtpcTrack = new TCanvas("TRDtracktpctrdtrack","TRDtracktpctrdtrack",10,10,510,510);
1008     debugtrdtpcTrack->Divide(2,1);
1009     debugtrdtpcTrack->cd(1);
1010     if(trdTrack) trdTrack->Draw();
1011     if(trdTrackOffline) trdTrackOffline->Draw("same");
1012     if(trdTrackStandalone) trdTrackStandalone->Draw("same");
1013     TLegend *leg = new TLegend(0.4,0.6,0.89,0.89);
1014     if(trdTrack) leg->AddEntry(trdTrack,"All","p");
1015     if(trdTrackOffline) leg->AddEntry(trdTrackOffline,"Offline","p");
1016     if(trdTrackStandalone) leg->AddEntry(trdTrackStandalone,"Standalone","p");
1017     leg->Draw("same");
1018     debugtrdtpcTrack->cd(2);
1019     if(tpctrdTrack) tpctrdTrack->Draw();
1020     TLine *line = new TLine(0.0,0.0,100.0,100.0);
1021     line->Draw("same");
1022
1023   }
1024
1025   if(nbTimeBin || nbTracklets || nbClusters) {
1026
1027     TCanvas *debugTracklets = new TCanvas("TRDtimebintrackletcluster","TRDtimebintrackletcluster",10,10,510,510);
1028     debugTracklets->Divide(3,1);
1029     debugTracklets->cd(1);
1030     if(nbTimeBin) nbTimeBin->Draw();
1031     if(nbTimeBinOffline) nbTimeBinOffline->Draw("same");
1032     if(nbTimeBinStandalone) nbTimeBinStandalone->Draw("same");
1033     TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1034     if(nbTimeBin) lega->AddEntry(nbTimeBin,"All","p");
1035     if(nbTimeBinOffline) lega->AddEntry(nbTimeBinOffline,"Offline","p");
1036     if(nbTimeBinStandalone) lega->AddEntry(nbTimeBinStandalone,"Standalone","p");
1037     lega->Draw("same");
1038     debugTracklets->cd(2);
1039     if(nbTracklets) nbTracklets->Draw();
1040     if(nbTrackletsOffline) nbTrackletsOffline->Draw("same");
1041     if(nbTrackletsStandalone) nbTrackletsStandalone->Draw("same");
1042     TLegend *legb = new TLegend(0.4,0.6,0.89,0.89);
1043     if(nbTracklets) legb->AddEntry(nbTracklets,"All","p");
1044     if(nbTrackletsOffline) legb->AddEntry(nbTrackletsOffline,"Offline","p");
1045     if(nbTrackletsStandalone) legb->AddEntry(nbTrackletsStandalone,"Standalone","p");
1046     legb->Draw("same");
1047     debugTracklets->cd(3);
1048     if(nbClusters) nbClusters->Draw();
1049     if(nbClustersOffline) nbClustersOffline->Draw("same");
1050     if(nbClustersStandalone) nbClustersStandalone->Draw("same");
1051     TLegend *legc = new TLegend(0.4,0.6,0.89,0.89);
1052     if(nbClusters) legc->AddEntry(nbClusters,"All","p");
1053     if(nbClustersOffline) legc->AddEntry(nbClustersOffline,"Offline","p");
1054     if(nbClustersStandalone) legc->AddEntry(nbClustersStandalone,"Standalone","p");
1055     legc->Draw("same");
1056
1057   }
1058
1059   if(ch2dSum || ph2dSum || histolinearfitsum) {
1060
1061     TCanvas *debugSum = new TCanvas("SumCalibrationObjects","SumCalibrationObjects",10,10,510,510);
1062     debugSum->Divide(3,1);
1063     debugSum->cd(1);
1064     if(ch2dSum) ch2dSum->Draw("lego");
1065     debugSum->cd(2);
1066     if(ph2dSum) ph2dSum->Draw("lego");
1067     debugSum->cd(3);
1068     if(histolinearfitsum) histolinearfitsum->Draw();
1069
1070   }
1071
1072   if(ch2dSM || ph2dSM) {
1073
1074     TCanvas *debugSM = new TCanvas("SMCalibrationObjects","SMCalibrationObjects",10,10,510,510);
1075     debugSM->Divide(2,1);
1076     debugSM->cd(1);
1077     if(ch2dSM) ch2dSM->Draw("lego");
1078     debugSM->cd(2);
1079     if(ph2dSM) ph2dSM->Draw("lego");
1080
1081   }
1082
1083   if(ch2d || ph2d) {
1084
1085     TCanvas *debug = new TCanvas("CalibrationObjects","CalibrationObjects",10,10,510,510);
1086     debug->Divide(2,1);
1087     debug->cd(1);
1088     if(ch2d) ch2d->Draw("lego");
1089     debug->cd(2);
1090     if(ph2d) ph2d->Draw("lego");
1091
1092   }
1093
1094 }
1095 //_______________________________________________________________________________________
1096 void  AliTRDCalibTask::AddTask(const AliTRDCalibTask * calibTask) {
1097
1098   //
1099   // Add stats
1100   //
1101
1102   TList *listcalibTask = calibTask->GetList();
1103   if(!listcalibTask) return;
1104
1105   TH1I *nEvents  = (TH1I *) listcalibTask->FindObject("NEvents");
1106   TH1I *nEventsInput  = (TH1I *) listcalibTask->FindObject("NEventsInput");
1107   TH2F *absoluteGain  = (TH2F *) listcalibTask->FindObject("AbsoluteGain");
1108
1109   TH1F *trdTrack = (TH1F *) listcalibTask->FindObject("TRDTrack");
1110   TH1F *trdTrackOffline = (TH1F *) listcalibTask->FindObject("TRDTrackOffline");
1111   TH1F *trdTrackStandalone = (TH1F *) listcalibTask->FindObject("TRDTrackStandalone");
1112
1113   TH2F *tpctrdTrack = (TH2F *) listcalibTask->FindObject("NbTPCTRDtrack");
1114
1115   TH1F *nbTimeBin = (TH1F *) listcalibTask->FindObject("NbTimeBin");
1116   TH1F *nbTimeBinOffline = (TH1F *) listcalibTask->FindObject("NbTimeBinOffline");
1117   TH1F *nbTimeBinStandalone = (TH1F *) listcalibTask->FindObject("NbTimeBinStandalone");
1118
1119   TH1F *nbClusters = (TH1F *) listcalibTask->FindObject("NbClusters");
1120   TH1F *nbClustersOffline = (TH1F *) listcalibTask->FindObject("NbClustersOffline");
1121   TH1F *nbClustersStandalone = (TH1F *) listcalibTask->FindObject("NbClustersStandalone");
1122
1123   TH1F *nbTracklets = (TH1F *) listcalibTask->FindObject("NbTracklets");
1124   TH1F *nbTrackletsOffline = (TH1F *) listcalibTask->FindObject("NbTrackletsOffline");
1125   TH1F *nbTrackletsStandalone = (TH1F *) listcalibTask->FindObject("NbTrackletsStandalone");
1126
1127   TH2I *ch2d = (TH2I *) listcalibTask->FindObject("CH2d");
1128   TProfile2D *ph2d = (TProfile2D *) listcalibTask->FindObject("PH2d");
1129   TProfile2D *prf2d = (TProfile2D *) listcalibTask->FindObject("PRF2d");
1130
1131   TH2I *ch2dSum = (TH2I *) listcalibTask->FindObject("CH2dSum");
1132   TProfile2D *ph2dSum = (TProfile2D *) listcalibTask->FindObject("PH2dSum");
1133
1134   TH2I *ch2dSM = (TH2I *) listcalibTask->FindObject("CH2dSM");
1135   TProfile2D *ph2dSM = (TProfile2D *) listcalibTask->FindObject("PH2dSM");
1136
1137   AliTRDCalibraVdriftLinearFit *linearfit = (AliTRDCalibraVdriftLinearFit *) listcalibTask->FindObject("AliTRDCalibraVdriftLinearFit");
1138   AliTRDCalibraVector *calibraVector = (AliTRDCalibraVector *) listcalibTask->FindObject("AliTRDCalibraVector");
1139
1140   //
1141
1142   TH1I *inEventsInput  = (TH1I *) fListHist->FindObject("NEventsInput");
1143   TH1I *inEvents  = (TH1I *) fListHist->FindObject("NEvents");
1144   TH2F *iabsoluteGain  = (TH2F *) fListHist->FindObject("AbsoluteGain");
1145
1146   TH1F *itrdTrack = (TH1F *) fListHist->FindObject("TRDTrack");
1147   TH1F *itrdTrackOffline = (TH1F *) fListHist->FindObject("TRDTrackOffline");
1148   TH1F *itrdTrackStandalone = (TH1F *) fListHist->FindObject("TRDTrackStandalone");
1149
1150   TH2F *itpctrdTrack = (TH2F *) fListHist->FindObject("NbTPCTRDtrack");
1151
1152   TH1F *inbTimeBin = (TH1F *) fListHist->FindObject("NbTimeBin");
1153   TH1F *inbTimeBinOffline = (TH1F *) fListHist->FindObject("NbTimeBinOffline");
1154   TH1F *inbTimeBinStandalone = (TH1F *) fListHist->FindObject("NbTimeBinStandalone");
1155
1156   TH1F *inbClusters = (TH1F *) fListHist->FindObject("NbClusters");
1157   TH1F *inbClustersOffline = (TH1F *) fListHist->FindObject("NbClustersOffline");
1158   TH1F *inbClustersStandalone = (TH1F *) fListHist->FindObject("NbClustersStandalone");
1159
1160   TH1F *inbTracklets = (TH1F *) fListHist->FindObject("NbTracklets");
1161   TH1F *inbTrackletsOffline = (TH1F *) fListHist->FindObject("NbTrackletsOffline");
1162   TH1F *inbTrackletsStandalone = (TH1F *) fListHist->FindObject("NbTrackletsStandalone");
1163
1164   TH2I *ich2d = (TH2I *) fListHist->FindObject("CH2d");
1165   TProfile2D *iph2d = (TProfile2D *) fListHist->FindObject("PH2d");
1166   TProfile2D *iprf2d = (TProfile2D *) fListHist->FindObject("PRF2d");
1167
1168   TH2I *ich2dSum = (TH2I *) fListHist->FindObject("CH2dSum");
1169   TProfile2D *iph2dSum = (TProfile2D *) fListHist->FindObject("PH2dSum");
1170
1171   TH2I *ich2dSM = (TH2I *) fListHist->FindObject("CH2dSM");
1172   TProfile2D *iph2dSM = (TProfile2D *) fListHist->FindObject("PH2dSM");
1173
1174   AliTRDCalibraVdriftLinearFit *ilinearfit = (AliTRDCalibraVdriftLinearFit *) fListHist->FindObject("AliTRDCalibraVdriftLinearFit");
1175   AliTRDCalibraVector *icalibraVector = (AliTRDCalibraVector *) fListHist->FindObject("AliTRDCalibraVector");
1176
1177
1178   // Add
1179
1180   if(nEventsInput) {
1181     if(inEventsInput) {
1182       inEventsInput->Add(nEventsInput);
1183       //printf("Add Events\n");
1184     }
1185     else {
1186       //printf("Create new Events\n");
1187       inEventsInput = new TH1I(*nEventsInput);
1188       fListHist->Add(inEventsInput);
1189     }
1190   }
1191
1192   if(nEvents) {
1193     if(inEvents) {
1194       inEvents->Add(nEvents);
1195       //printf("Add Events\n");
1196     }
1197     else {
1198       //printf("Create new Events\n");
1199       inEvents = new TH1I(*nEvents);
1200       fListHist->Add(inEvents);
1201     }
1202   }
1203
1204   if(absoluteGain) {
1205     if(iabsoluteGain) iabsoluteGain->Add(absoluteGain);
1206     else {
1207       iabsoluteGain = new TH2F(*absoluteGain);
1208       fListHist->Add(iabsoluteGain);
1209     }
1210   }
1211
1212   if(trdTrack) {
1213     if(itrdTrack) itrdTrack->Add(trdTrack);
1214     else {
1215      itrdTrack = new TH1F(*trdTrack);
1216      fListHist->Add(itrdTrack);
1217     }
1218   }
1219
1220   if(trdTrackOffline) {
1221     if(itrdTrackOffline) itrdTrackOffline->Add(trdTrackOffline);
1222     else {
1223       itrdTrackOffline = new TH1F(*trdTrackOffline);
1224       fListHist->Add(itrdTrackOffline);
1225     }
1226   }
1227
1228   if(trdTrackStandalone) {
1229     if(itrdTrackStandalone) itrdTrackStandalone->Add(trdTrackStandalone);
1230     else {
1231       itrdTrackStandalone = new TH1F(*trdTrackStandalone);
1232       fListHist->Add(itrdTrackStandalone);
1233     }
1234   }
1235
1236   if(tpctrdTrack) {
1237     if(itpctrdTrack) itpctrdTrack->Add(tpctrdTrack);
1238     else {
1239       itpctrdTrack = new TH2F(*tpctrdTrack);
1240       fListHist->Add(itpctrdTrack);
1241     }
1242   }
1243
1244   if(nbTimeBin) {
1245     if(inbTimeBin) inbTimeBin->Add(nbTimeBin);
1246     else {
1247       inbTimeBin = new TH1F(*inbTimeBin);
1248       fListHist->Add(inbTimeBin);
1249     }
1250   }
1251
1252   if(nbTimeBinOffline) {
1253     if(inbTimeBinOffline) inbTimeBinOffline->Add(nbTimeBinOffline);
1254     else {
1255       inbTimeBinOffline = new TH1F(*nbTimeBinOffline);
1256       fListHist->Add(inbTimeBinOffline);
1257     }
1258   }
1259
1260   if(nbTimeBinStandalone) {
1261     if(inbTimeBinStandalone) inbTimeBinStandalone->Add(nbTimeBinStandalone);
1262     else {
1263       inbTimeBinStandalone = new TH1F(*nbTimeBinStandalone);
1264       fListHist->Add(inbTimeBinStandalone);
1265     }
1266   }
1267
1268   if(nbClusters) {
1269     if(inbClusters) inbClusters->Add(nbClusters);
1270     else {
1271       inbClusters = new TH1F(*nbClusters);
1272       fListHist->Add(inbClusters);
1273     }
1274   }
1275
1276   if(nbClustersOffline) {
1277     if(inbClustersOffline) inbClustersOffline->Add(nbClustersOffline);
1278     else {
1279       inbClustersOffline = new TH1F(*nbClustersOffline);
1280       fListHist->Add(inbClustersOffline);
1281     }
1282   }
1283
1284   if(nbClustersStandalone) {
1285     if(inbClustersStandalone) inbClustersStandalone->Add(nbClustersStandalone);
1286     else {
1287       inbClustersStandalone = new TH1F(*nbClustersStandalone);
1288       fListHist->Add(inbClustersStandalone);
1289     }
1290   }
1291
1292   if(nbTracklets) {
1293     if(inbTracklets) inbTracklets->Add(nbTracklets);
1294     else {
1295       inbTracklets = new TH1F(*nbTracklets);
1296       fListHist->Add(inbTracklets);
1297     }
1298   }
1299
1300   if(nbTrackletsOffline) {
1301     if(inbTrackletsOffline) inbTrackletsOffline->Add(nbTrackletsOffline);
1302     else {
1303       inbTrackletsOffline = new TH1F(*nbTrackletsOffline);
1304       fListHist->Add(inbTrackletsOffline);
1305     }
1306   }
1307
1308   if(nbTrackletsStandalone) {
1309     if(inbTrackletsStandalone) inbTrackletsStandalone->Add(nbTrackletsStandalone);
1310     else {
1311       inbTrackletsStandalone = new TH1F(*nbTrackletsStandalone);
1312       fListHist->Add(inbTrackletsStandalone);
1313     }
1314   }
1315
1316   if(ch2d) {
1317     if(ich2d) ich2d->Add(ch2d);
1318     else {
1319       ich2d = new TH2I(*ch2d);
1320       fListHist->Add(ich2d);
1321     }
1322   }
1323
1324   if(ph2d) {
1325     if(iph2d) iph2d->Add(ph2d);
1326     else {
1327       iph2d = new TProfile2D(*ph2d);
1328       fListHist->Add(iph2d);
1329     }
1330   }
1331
1332   if(prf2d) {
1333     if(iprf2d) iprf2d->Add(prf2d);
1334     else {
1335       iprf2d = new TProfile2D(*prf2d);
1336       fListHist->Add(iprf2d);
1337     }
1338   }
1339
1340   if(ch2dSum) {
1341     if(ich2dSum) ich2dSum->Add(ch2dSum);
1342     else {
1343       ich2dSum = new TH2I(*ch2dSum);
1344       fListHist->Add(ich2dSum);
1345     }
1346   }
1347
1348   if(ph2dSum) {
1349     if(iph2dSum) iph2dSum->Add(ph2dSum);
1350     else {
1351       iph2dSum = new TProfile2D(*ph2dSum);
1352       fListHist->Add(iph2dSum);
1353     }
1354   }
1355
1356   if(ch2dSM) {
1357     if(ich2dSM) ich2dSM->Add(ch2dSM);
1358     else {
1359       ich2dSM = new TH2I(*ch2dSM);
1360       fListHist->Add(ich2dSM);
1361     }
1362   }
1363
1364   if(ph2dSM) {
1365     if(iph2dSM) iph2dSM->Add(ph2dSM);
1366     else {
1367       iph2dSM = new TProfile2D(*ph2dSM);
1368       fListHist->Add(iph2dSM);
1369     }
1370   }
1371
1372   if(linearfit) {
1373     if(ilinearfit) ilinearfit->Add(linearfit);
1374     else {
1375       ilinearfit = new AliTRDCalibraVdriftLinearFit(*linearfit);
1376       fListHist->Add(ilinearfit);
1377     }
1378   }
1379
1380   if(calibraVector) {
1381     if(icalibraVector) icalibraVector->Add(calibraVector);
1382     else {
1383       icalibraVector = new AliTRDCalibraVector(*calibraVector);
1384       fListHist->Add(icalibraVector);
1385     }
1386   }
1387
1388 }
1389 //________________________________________________________________________________
1390 Long64_t AliTRDCalibTask::Merge(TCollection *li) {
1391
1392   //
1393   // merge component
1394   //
1395
1396   TIterator* iter = li->MakeIterator();
1397   AliTRDCalibTask* cal = 0;
1398
1399   while ((cal = (AliTRDCalibTask*)iter->Next())) {
1400     if (!cal->InheritsFrom(AliTRDCalibTask::Class())) {
1401       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
1402       return -1;
1403     }
1404
1405     // add histograms here...
1406     this->AddTask(cal);
1407
1408   }
1409
1410   return 0;
1411
1412 }
1413 //_____________________________________________________
1414 Bool_t AliTRDCalibTask::SetVersionSubversion(){
1415   //
1416   // Load Chamber Gain factors into the Tender supply
1417   //
1418
1419   printf("SetVersionSubversion\n");
1420
1421   //find previous entry from the UserInfo
1422   TTree *tree=((TChain*)GetInputData(0))->GetTree();
1423   if (!tree) {
1424     AliError("Tree not found in ESDhandler");
1425     return kFALSE;
1426   }
1427
1428   TList *userInfo=(TList*)tree->GetUserInfo();
1429   if (!userInfo) {
1430     AliError("No UserInfo found in tree");
1431     return kFALSE;
1432   }
1433
1434   TList *cdbList=(TList*)userInfo->FindObject("cdbList");
1435   if (!cdbList) {
1436     AliError("No cdbList found in UserInfo");
1437     if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print();
1438     return kFALSE;
1439   }
1440
1441   TIter nextCDB(cdbList);
1442   TObjString *os=0x0;
1443   while ( (os=(TObjString*)nextCDB()) ){
1444     if(os->GetString().Contains("TRD/Calib/ChamberGainFactor")){
1445       // Get Old gain calibration
1446       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1447       fFirstRunGain = id->GetFirstRun();
1448       fVersionGainUsed = id->GetVersion();
1449       fSubVersionGainUsed = id->GetSubVersion();
1450     } else if(os->GetString().Contains("TRD/Calib/ChamberVdrift")){
1451       // Get Old drift velocity calibration
1452       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1453       fFirstRunVdrift = id->GetFirstRun();
1454       fVersionVdriftUsed = id->GetVersion();
1455       fSubVersionVdriftUsed = id->GetSubVersion();
1456     } else if(os->GetString().Contains("TRD/Calib/LocalGainFactor")){
1457       // Get Old drift velocity calibration
1458       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
1459       fFirstRunGainLocal = id->GetFirstRun();
1460       fVersionGainLocalUsed = id->GetVersion();
1461       fSubVersionGainLocalUsed = id->GetSubVersion();
1462     }
1463   }
1464
1465   //printf("VersionGain %d, SubversionGain %d, VersionLocalGain %d, Subversionlocalgain %d, Versionvdrift %d, Subversionvdrift %d\n",fVersionGainUsed,fSubVersionGainUsed,fVersionGainLocalUsed,fSubVersionGainLocalUsed,fVersionVdriftUsed,fSubVersionVdriftUsed);
1466
1467   // Check
1468   if((fFirstRunGain < 0)            ||
1469      (fFirstRunGainLocal < 0)       ||
1470      (fFirstRunVdrift < 0)          ||
1471      (fVersionGainUsed < 0)         ||
1472      (fVersionGainLocalUsed < 0)    ||
1473      (fSubVersionGainUsed < 0)      ||
1474      (fSubVersionGainLocalUsed < 0) ||
1475      (fVersionVdriftUsed < 0)       ||
1476      (fSubVersionVdriftUsed < 0)) {
1477     AliError("No recent calibration found");
1478     return kFALSE;
1479   }
1480   else return kTRUE;
1481
1482 }
1483 //_________________________________________________________________________________________________________________________
1484 Bool_t AliTRDCalibTask::ParticleGood(int i) const {
1485
1486   //
1487   // Definition of good tracks
1488   //
1489
1490
1491   AliESDtrack *track = fESD->GetTrack(i);
1492   if (!track->IsOn(AliESDtrack::kTPCrefit)) return 0;        // TPC refit
1493   if (track->GetTPCNcls() < 90) return 0;                    // number of TPC clusters
1494   if (fabs(track->Eta())>0.8) return 0;                         // fiducial pseudorapidity
1495   Float_t r,z;
1496   track->GetImpactParametersTPC(r,z);
1497   if (fabs(z)>2.0) return 0;                          // impact parameter in z
1498   if (fabs(r)>2.0) return 0;                          // impact parameter in xy
1499   if (r==0) return 0;
1500   return 1;
1501
1502
1503 }
1504
1505
1506