TRD/AliTRDCalibTask.cxx

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