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