PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetJTJT.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 //
  16 //
  17 // Jet fragmentation transverse momentum (j_T) analysis task
  18 //
  19 // Author: T.Snellman
  20
  21 #include <TClonesArray.h>
  22 #include <TH1F.h>
  23 #include <TH2F.h>
  24 #include <TH3F.h>
  25 #include <TList.h>
  26 #include <TProfile.h>
  27 #include <TLorentzVector.h>
  28 #include <TVector.h>
  29 #include <TGraphErrors.h>
  30 #include <TGrid.h>
  31 #include <TSystem.h>
  32 #include <TFile.h>
  33
  34 #include "AliCentrality.h"
  35
  36
  37
  38 #include "AliVCluster.h"
  39 #include "AliAODCaloCluster.h"
  40 #include "AliESDCaloCluster.h"
  41 #include "AliVTrack.h"
  42 #include "AliEmcalJet.h"
  43 #include "AliRhoParameter.h"
  44 #include "AliLog.h"
  45 #include "AliJetContainer.h"
  46 #include "AliParticleContainer.h"
  47 #include "AliClusterContainer.h"
  48 #include "AliPicoTrack.h"
  49
  50 #include "AliAnalysisTaskJetJTJT.h"
  51
  52
  53 ClassImp(AliAnalysisTaskJetJTJT)
  54
  55         //________________________________________________________________________
  56         AliAnalysisTaskJetJTJT::AliAnalysisTaskJetJTJT() :
  57                 AliAnalysisTaskEmcalJet("AliAnalysisTaskJetJTJT", kTRUE),
  58                 fHistTracksPt(0),
  59                 fHistTracksJt(0),
  60                 fHistClustersPt(0),
  61                 fHistLeadingJetPt(0),
  62                 fHistJetsPt(0),
  63                 fHistBackgroundDone(0),
  64                 fHistJTPta(0),
  65                 fHistLogJTPta(0),
  66                 fHistJTPta_all(0),
  67                 fHistJTBg(0),
  68                 fHistLogJTBg(0),
  69                 fHistBgMulti(0),
  70                 fHistBgPt(0),
  71                 fHistJetEta(0),
  72                 fHistJetMulti(0),
  73                 fHistJetTracksPt(0),
  74                 fhTrackingEfficiency(0),
  75                 fNpttBins(1),
  76                 fNptaBins(1),
  77                 fJetsCont(0),
  78                 //fJetsConts(0),
  79                 //nJetsConts(0),
  80                 fTracksCont(0),
  81                 fCaloClustersCont(0),
  82                 fTracks(0),
  83                 fTrackArrayName("nonejk"),
  84                 runPeriod(""),
  85                 fEfficiency(0),
  86                 debug(0)
  87
  88
  89 {
  90         // Default constructor.
  91
  92         fHistTracksPt       = new TH1*[fNcentBins];
  93         fHistTracksJt       = new TH1*[fNcentBins];
  94         fHistClustersPt     = new TH1*[fNcentBins];
  95         fHistLeadingJetPt   = new TH1*[fNcentBins];
  96         fHistJetsPt         = new TH1**[fNcentBins];
  97         fHistBackgroundDone = new TH1**[fNcentBins];
  98         fHistJTPta          = new TH1***[fNcentBins];
  99         fHistLogJTPta       = new TH1***[fNcentBins];
 100         fHistJTPta_all      = new TH1***[fNcentBins];
 101         fHistJTBg           = new TH1***[fNcentBins];
 102         fHistLogJTBg        = new TH1***[fNcentBins];
 103         fHistBgMulti        = new TH1**[fNcentBins];
 104         fHistBgPt           = new TH1**[fNcentBins];
 105         fHistJetEta         = new TH1**[fNcentBins];
 106         fHistJetMulti       = new TH1**[fNcentBins];
 107         fHistJetTracksPt     = new TH1**[fNcentBins];
 108         fhTrackingEfficiency = new TProfile*[fNcentBins];
 109         //CentBinBorders      = new Double_t[10];
 110
 111
 112         for (Int_t i = 0; i < fNcentBins; i++) {
 113                 fHistJTPta[i] = 0;
 114                 fHistLogJTPta[i] = 0;
 115                 fHistJTPta_all[i] = 0;
 116                 fHistJTBg[i] = 0;
 117                 fHistLogJTBg[i] = 0;
 118                 fHistBackgroundDone[i] = 0;
 119                 fHistTracksPt[i] = 0;
 120                 fHistTracksJt[i] = 0;
 121                 fHistClustersPt[i] = 0;
 122                 fHistLeadingJetPt[i] = 0;
 123                 fHistJetsPt[i] = 0;
 124                 fHistBgMulti[i] = 0;
 125                 fHistBgPt[i] = 0;
 126                 fHistJetEta[i] = 0;
 127                 fHistJetMulti[i] = 0;
 128                 fHistJetTracksPt[i] = 0;
 129                 fhTrackingEfficiency[i] = 0;
 130         }
 131
 132         /*for(Int_t i = 0; i < nJetsConts; i++){
 133           fJetsConts[i] = 0;
 134           }*/
 135         SetMakeGeneralHistograms(kTRUE);
 136 }
 137
 138 //________________________________________________________________________
 139 AliAnalysisTaskJetJTJT::AliAnalysisTaskJetJTJT(const char *name) :
 140         AliAnalysisTaskEmcalJet(name, kTRUE),
 141         fHistTracksPt(0),
 142         fHistTracksJt(0),
 143         fHistClustersPt(0),
 144         fHistLeadingJetPt(0),
 145         fHistJetsPt(0),
 146         fHistBackgroundDone(0),
 147         fHistJTPta(0),
 148         fHistLogJTPta(0),
 149         fHistJTPta_all(0),
 150         fHistJTBg(0),
 151         fHistLogJTBg(0),
 152         fHistBgMulti(0),
 153         fHistBgPt(0),
 154         fHistJetEta(0),
 155         fHistJetMulti(0),
 156         fHistJetTracksPt(0),
 157         fhTrackingEfficiency(0),
 158         fNpttBins(1),
 159         fNptaBins(1),
 160         fJetsCont(0),
 161         //fJetsConts(0),
 162         //nJetsConts(0),
 163         fTracksCont(0),
 164         fCaloClustersCont(0),
 165         fTracks(0),
 166         fTrackArrayName("nonejk"),
 167         runPeriod(""),
 168         fEfficiency(0),
 169         debug(0)
 170 {
 171         // Standard constructor.
 172         fHistTracksPt       = new TH1*[fNcentBins];
 173         fHistTracksJt       = new TH1*[fNcentBins];
 174         fHistClustersPt     = new TH1*[fNcentBins];
 175         fHistLeadingJetPt   = new TH1*[fNcentBins];
 176         fHistJetsPt         = new TH1**[fNcentBins];
 177         fHistBackgroundDone = new TH1**[fNcentBins];
 178         fHistJTPta          = new TH1***[fNcentBins];
 179         fHistLogJTPta       = new TH1***[fNcentBins];
 180         fHistJTPta_all      = new TH1***[fNcentBins];
 181         fHistJTBg           = new TH1***[fNcentBins];
 182         fHistLogJTBg        = new TH1***[fNcentBins];
 183         fHistBgMulti        = new TH1**[fNcentBins];
 184         fHistBgPt           = new TH1**[fNcentBins];
 185         fHistJetEta         = new TH1**[fNcentBins];
 186         fHistJetMulti       = new TH1**[fNcentBins];
 187         fHistJetTracksPt     = new TH1**[fNcentBins];
 188         fhTrackingEfficiency = new TProfile*[fNcentBins];
 189         //CentBinBorders            = new Double_t[10];
 190
 191
 192         for (Int_t i = 0; i < fNcentBins; i++) {
 193                 fHistJTPta[i] = 0;
 194                 fHistLogJTPta[i] = 0;
 195                 fHistJTPta_all[i] = 0;
 196                 fHistJTBg[i] = 0;
 197                 fHistLogJTBg[i] = 0;
 198                 fHistBackgroundDone[i] = 0;
 199                 fHistTracksPt[i] = 0;
 200                 fHistTracksJt[i] = 0;
 201                 fHistClustersPt[i] = 0;
 202                 fHistLeadingJetPt[i] = 0;
 203                 fHistJetsPt[i] = 0;
 204                 fHistBgMulti[i] = 0;
 205                 fHistBgPt[i] = 0;
 206                 fHistJetEta[i] = 0;
 207                 fHistJetMulti[i] = 0;
 208                 fHistJetTracksPt[i] = 0;
 209                 fhTrackingEfficiency[i] = 0;
 210         }
 211
 212         /*for(Int_t i = 0; i < nJetsConts; i++){
 213           fJetsConts[i] = 0;
 214           }*/
 215         SetMakeGeneralHistograms(kTRUE);
 216 }
 217
 218 //________________________________________________________________________
 219 AliAnalysisTaskJetJTJT::~AliAnalysisTaskJetJTJT()
 220 {
 221         // Destructor.
 222 }
 223
 224
 225 void AliAnalysisTaskJetJTJT::setCentBinBorders( int n, Double_t *c){
 226         fNcentBins=n;
 227         if(debug > 0){
 228                 cout << "AliAnalysisTaskJetJTJT::setCentBinBorders: " << endl;
 229         }
 230         for(int i= 0 ; i < fNcentBins; i++){
 231                 CentBinBorders[i]= c[i];
 232                 if(debug > 0)
 233                         cout << CentBinBorders[i] << endl;
 234         }
 235 }
 236
 237 void AliAnalysisTaskJetJTJT::setTriggPtBorders( int n, Double_t *c){
 238         fNpttBins=n;
 239         if(debug > 0)
 240                 cout << "AliAnalysisTaskJetJTJT::setTriggPtBorders: " << endl;
 241         for(int i= 0 ; i < fNpttBins; i++){
 242                 TriggPtBorders[i]= c[i];
 243                 if(debug > 0)
 244                         cout << TriggPtBorders[i] << endl;
 245         }
 246 }
 247
 248 void AliAnalysisTaskJetJTJT::setAssocPtBorders( int n, Double_t *c){
 249         fNptaBins=n;
 250         if(debug > 0)
 251                 cout << "AliAnalysisTaskJetJTJT::setAssocPtBorders: " << endl;
 252         for(int i= 0 ; i < fNptaBins; i++){
 253                 AssocPtBorders[i]= c[i];
 254                 if(debug > 0)
 255                         cout << AssocPtBorders[i] << endl;
 256         }
 257 }
 258
 259
 260 //________________________________________________________________________
 261 void AliAnalysisTaskJetJTJT::UserCreateOutputObjects()
 262 {
 263         // Create user output.
 264
 265         AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 266         if(debug > 0)
 267                 cout << "Creating Histograms" << endl;
 268
 269         fJetsCont           = GetJetContainer(0);
 270         /*for(int i = 0; i < nJetsConts ; i++){
 271           fJetsConts[0]     = GetJetContainer(i);
 272           }*/
 273         /*if(fJetsCont) { //get particles and clusters connected to jets
 274           fTracksCont       = fJetsCont->GetParticleContainer();
 275           fCaloClustersCont = fJetsCont->GetClusterContainer();
 276           } else {        //no jets, just analysis tracks and clusters
 277           fTracksCont       = GetParticleContainer(0);
 278           fCaloClustersCont = GetClusterContainer(0);
 279           }*/
 280         fTracksCont       = GetParticleContainer(0);
 281         fCaloClustersCont = GetClusterContainer(0);
 282         fTracksCont->SetClassName("AliVTrack");
 283         fCaloClustersCont->SetClassName("AliAODCaloCluster");
 284
 285         TString histname;
 286         //Int_t fMinBinJt = 0;
 287         //Int_t fMaxBinJt = 5;
 288
 289         Int_t NBINSJt=150;
 290         double LogBinsJt[NBINSJt+1], LimLJt=0.01, LimHJt=10;
 291         double logBWJt = (TMath::Log(LimHJt)-TMath::Log(LimLJt))/(NBINSJt-1);
 292         LogBinsJt[0] = 0;
 293         for(int ij=1;ij<=NBINSJt;ij++) LogBinsJt[ij]=LimLJt*exp(ij*logBWJt);
 294
 295         int NBINSJtW=150;
 296         double LimLJtW=TMath::Log(0.01), LimHJtW=TMath::Log(10);
 297
 298         //==== Efficiency ====
 299         if(debug > 0)
 300                 cout << "AliAnalysisTaskJetJTJT::UserCreateOutputObjects: Creating efficiency" << endl;
 301         fEfficiency = new JTJTEfficiency;
 302         fEfficiency->SetMode( 1 ); // 0:NoEff, 1:Period 2:RunNum 3:Auto
 303         fEfficiency->SetDataPath("alien:///alice/cern.ch/user/d/djkim/legotrain/efficieny/data"); // Efficiency root file location local or alien
 304
 305
 306
 307         for (Int_t ic = 0; ic < fNcentBins; ic++) {
 308                 if (fParticleCollArray.GetEntriesFast()>0) {
 309                         histname = "fHistTracksPt_";
 310                         histname += ic;
 311                         fHistTracksPt[ic] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt / 4);
 312                         fHistTracksPt[ic]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
 313                         fHistTracksPt[ic]->GetYaxis()->SetTitle("counts");
 314                         fOutput->Add(fHistTracksPt[ic]);
 315                 }
 316
 317                 if (fParticleCollArray.GetEntriesFast()>0) {
 318                         histname = "fHistTracksJt_";
 319                         histname += ic;
 320                         fHistTracksJt[ic] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
 321                         fHistTracksJt[ic]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
 322                         fHistTracksJt[ic]->GetYaxis()->SetTitle("counts");
 323                         fOutput->Add(fHistTracksJt[ic]);
 324                 }
 325
 326                 if (fParticleCollArray.GetEntriesFast()>0) {
 327                         histname = "fhTrackingEfficiency_";
 328                         histname += ic;
 329                         fhTrackingEfficiency[ic] = new TProfile(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 330                         fhTrackingEfficiency[ic]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
 331                         fhTrackingEfficiency[ic]->GetYaxis()->SetTitle("counts");
 332                         fOutput->Add(fhTrackingEfficiency[ic]);
 333                 }
 334                 fHistJTPta[ic] = new TH1**[fNpttBins];
 335                 fHistLogJTPta[ic] = new TH1**[fNpttBins];
 336                 fHistJTPta_all[ic] = new TH1**[fNpttBins];
 337                 fHistJTBg[ic] = new TH1**[fNpttBins];
 338                 fHistLogJTBg[ic] = new TH1**[fNpttBins];
 339                 fHistJetsPt[ic] = new TH1*[fNpttBins];
 340                 fHistBackgroundDone[ic] = new TH1*[fNpttBins];
 341                 fHistBgMulti[ic] = new TH1*[fNpttBins];
 342                 fHistBgPt[ic] = new TH1*[fNpttBins];
 343                 fHistJetEta[ic] = new TH1*[fNpttBins];
 344                 fHistJetMulti[ic] = new TH1*[fNpttBins];
 345                 fHistJetTracksPt[ic] = new TH1*[fNpttBins];
 346                 for(Int_t j=0; j < fNpttBins; j++){
 347                         fHistJTPta[ic][j] = new TH1*[fNptaBins];
 348                         fHistLogJTPta[ic][j] = new TH1*[fNptaBins];
 349                         fHistJTPta_all[ic][j] = new TH1*[fNptaBins];
 350                         fHistJTBg[ic][j] = new TH1*[fNptaBins];
 351                         fHistLogJTBg[ic][j] = new TH1*[fNptaBins];
 352                         for(Int_t k=0; k < fNptaBins; k++){
 353                                 fHistJTPta[ic][j][k] = 0;
 354                                 fHistLogJTPta[ic][j][k] = 0;
 355                                 fHistJTPta_all[ic][j][k] = 0;
 356                                 fHistJTBg[ic][j][k] = 0;
 357                                 fHistLogJTBg[ic][j][k] = 0;
 358                         }
 359                         fHistJetsPt[ic][j] = 0;
 360                         fHistBackgroundDone[ic][j] = 0;
 361                         fHistBgMulti[ic][j] = 0;
 362                         fHistBgPt[ic][j] = 0;
 363                         fHistJetEta[ic][j] = 0;
 364                         fHistJetMulti[ic][j] =0;
 365                         fHistJetTracksPt[ic][j] = 0;
 366                 }
 367
 368
 369                 if (fParticleCollArray.GetEntriesFast()>0) {
 370                         for(Int_t iptt = 0 ; iptt <  fNpttBins; iptt++){
 371                                 for(Int_t ipta = 0 ; ipta < fNptaBins; ipta++){
 372                                         histname = "hJTPtaD00C";
 373                                         //histname += ic;
 374                                         histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
 375                                         if(debug > 1)
 376                                                 cout << histname << endl;
 377                                         fHistJTPta[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(),NBINSJt, LogBinsJt);
 378                                         fHistJTPta[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
 379                                         fHistJTPta[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
 380                                         fOutput->Add(fHistJTPta[ic][iptt][ipta]);
 381
 382                                         histname = "hLogJTPtaD00C";
 383                                         //histname += ic;
 384                                         histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
 385                                         if(debug > 1)
 386                                                 cout << histname << endl;
 387                                         fHistLogJTPta[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJtW, LimLJtW, LimHJtW);
 388                                         fHistLogJTPta[ic][iptt][ipta]->GetXaxis()->SetTitle("ln(J_{T,track}/ (GeV/c))");
 389                                         fHistLogJTPta[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
 390                                         fOutput->Add(fHistLogJTPta[ic][iptt][ipta]);
 391
 392                                         histname = "hJTPta_allD00C";
 393                                         //histname += ic;
 394                                         histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
 395                                         if(debug > 1)
 396                                                 cout << histname << endl;
 397                                         fHistJTPta_all[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
 398                                         fHistJTPta_all[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
 399                                         fHistJTPta_all[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
 400                                         fOutput->Add(fHistJTPta_all[ic][iptt][ipta]);
 401
 402                                         histname = "hJTBg";
 403                                         //histname += ic;
 404                                         histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
 405                                         if(debug > 1)
 406                                                 cout << histname << endl;
 407                                         fHistJTBg[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
 408                                         fHistJTBg[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
 409                                         fHistJTBg[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
 410                                         fOutput->Add(fHistJTBg[ic][iptt][ipta]);
 411
 412                                         histname = "hLogJTBg";
 413                                         //histname += ic;
 414                                         histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
 415                                         if(debug > 1)
 416                                                 cout << histname << endl;
 417                                         fHistLogJTBg[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJtW, LimLJtW, LimHJtW);
 418                                         fHistLogJTBg[ic][iptt][ipta]->GetXaxis()->SetTitle("ln(J_{T,track}/ (GeV/c))");
 419                                         fHistLogJTBg[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
 420                                         fOutput->Add(fHistLogJTBg[ic][iptt][ipta]);
 421
 422
 423                                 }
 424                         }
 425                 }
 426
 427                 if (fClusterCollArray.GetEntriesFast()>0) {
 428                         histname = "fHistClustersPt_";
 429                         histname += ic;
 430                         fHistClustersPt[ic] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 431                         fHistClustersPt[ic]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
 432                         fHistClustersPt[ic]->GetYaxis()->SetTitle("counts");
 433                         fOutput->Add(fHistClustersPt[ic]);
 434                 }
 435
 436                 if (fJetCollArray.GetEntriesFast()>0) {
 437                         histname = "fHistLeadingJetPt_";
 438                         histname += ic;
 439                         fHistLeadingJetPt[ic] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 440                         fHistLeadingJetPt[ic]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 441                         fHistLeadingJetPt[ic]->GetYaxis()->SetTitle("counts");
 442                         fOutput->Add(fHistLeadingJetPt[ic]);
 443
 444                         for(Int_t iptt = 0 ; iptt <  fNpttBins; iptt++){
 445
 446                                 histname = "fHistJetsPt_";
 447                                 histname += Form("C%02dT%02d", ic, iptt);
 448                                 fHistJetsPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 449                                 fHistJetsPt[ic][iptt]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 450                                 fHistJetsPt[ic][iptt]->GetYaxis()->SetTitle("counts");
 451                                 fOutput->Add(fHistJetsPt[ic][iptt]);
 452
 453                                 histname = "fHistBackgroundDone_";
 454                                 histname += Form("C%02dT%02d", ic, iptt);;
 455                                 fHistBackgroundDone[ic][iptt] = new TH1F(histname.Data(), histname.Data(), 2, -1, 2);
 456                                 fHistBackgroundDone[ic][iptt]->GetXaxis()->SetTitle("Number of jets");
 457                                 fHistBackgroundDone[ic][iptt]->GetYaxis()->SetTitle("0 = not done, 1 = done");
 458                                 fOutput->Add(fHistBackgroundDone[ic][iptt]);
 459
 460                                 histname = "fHistJetEta_";
 461                                 histname += Form("C%02dT%02d", ic, iptt);
 462                                 fHistJetEta[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, -5, 5);
 463                                 fHistJetEta[ic][iptt]->GetXaxis()->SetTitle("#eta");
 464                                 fHistJetEta[ic][iptt]->GetYaxis()->SetTitle("jets");
 465                                 fOutput->Add(fHistJetEta[ic][iptt]);
 466
 467                                 histname = "fHistJetMulti_";
 468                                 histname += Form("C%02dT%02d", ic, iptt);
 469                                 fHistJetMulti[ic][iptt] = new TH1F(histname.Data(), histname.Data(), 200, 0, 200);
 470                                 fHistJetMulti[ic][iptt]->GetXaxis()->SetTitle("Multiplicity");
 471                                 fHistJetMulti[ic][iptt]->GetYaxis()->SetTitle("jets");
 472                                 fOutput->Add(fHistJetMulti[ic][iptt]);
 473
 474                                 histname = "fHistBgMulti_";
 475                                 histname += Form("C%02dT%02d", ic, iptt);
 476                                 fHistBgMulti[ic][iptt] = new TH1F(histname.Data(), histname.Data(), 200, 0, 200);
 477                                 fHistBgMulti[ic][iptt]->GetXaxis()->SetTitle("Multiplicity");
 478                                 fHistBgMulti[ic][iptt]->GetYaxis()->SetTitle("Events");
 479                                 fOutput->Add(fHistBgMulti[ic][iptt]);
 480
 481                                 histname = "fHistBgPt_";
 482                                 histname += Form("C%02dT%02d", ic, iptt);
 483                                 fHistBgPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins,fMinBinPt, fMaxBinPt/4);
 484                                 fHistBgPt[ic][iptt]->GetXaxis()->SetTitle("Multiplicity");
 485                                 fHistBgPt[ic][iptt]->GetYaxis()->SetTitle("jets");
 486                                 fOutput->Add(fHistBgPt[ic][iptt]);
 487
 488                                 histname = "fHistJetTracksPt_";
 489                                 //histname += ic;
 490                                 histname += Form("C%02dT%02d", ic, iptt);
 491                                 if(debug > 1)
 492                                         cout << histname << endl;
 493                                 fHistJetTracksPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt/4);
 494                                 fHistJetTracksPt[ic][iptt]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
 495                                 fHistJetTracksPt[ic][iptt]->GetYaxis()->SetTitle("counts");
 496                                 fOutput->Add(fHistJetTracksPt[ic][iptt]);
 497                         }
 498
 499                         /*
 500                            if (!(GetJetContainer()->GetRhoName().IsNull())) {
 501                            histname = "fHistJetsCorrPtArea_";
 502                            histname += i;
 503                            fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
 504                            fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
 505                            fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
 506                            fOutput->Add(fHistJetsCorrPtArea[i]);
 507                            }
 508                            */
 509                 }
 510         }
 511
 512         PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 513 }
 514
 515 //________________________________________________________________________
 516 Bool_t AliAnalysisTaskJetJTJT::FillHistograms()
 517 {
 518         // Fill histograms.
 519
 520         AliCentrality *aliCent = InputEvent()->GetCentrality();
 521         fCentBin = 0;
 522         if (aliCent) {
 523                 //fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
 524                 fCent = aliCent->GetCentralityPercentile("V0M");
 525                 /*if(debug > 0){
 526                   cout << "Centrality " << fCent << endl;
 527                   }*/
 528                 for(int ic = 0; ic < fNcentBins; ic++){
 529                         /*if(debug > 0){
 530                           cout << "ic: " << ic << " / " << fNcentBins << endl;
 531                           cout << "Centrality bin " << fCentBin << endl;
 532                           cout << "Border: " << CentBinBorders[ic] << endl;
 533                           } */
 534                         if(fCent > CentBinBorders[ic]){
 535                                 fCentBin = ic;
 536                         }
 537                 }
 538                 //cout << "Centrality bin: " << fCentBin << endl;
 539         } else {
 540                 AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
 541                 fCentBin = 3;
 542         }
 543         int fHadronSelectionCut = 5; //5=Hybrid cut
 544         if (fTracksCont) {
 545                 AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
 546                 while(track) {
 547                         double ptt = track->Pt();
 548
 549                         //<<<<<<<<<<<< Efficiency >>>>>>>>>>>
 550                         //AliJBaseTrack *triggTr = (AliJBaseTrack*)fchargedHadronList->At(i);
 551
 552                         double effCorr = 1./fEfficiency->GetCorrection(ptt, fHadronSelectionCut, fCent);  // here you generate warning if ptt>30
 553                         //double effCorr = 1.;
 554                         fhTrackingEfficiency[fCentBin]->Fill( ptt, 1./effCorr );
 555                         //triggTr->SetTrackEff( 1./effCorr );
 556                         //<<<<<<<<<<<< Efficiency >>>>>>>>>>>
 557
 558                         if(ptt > 0 && 1.0/ptt > 0){
 559                                 fHistTracksPt[fCentBin]->Fill(ptt,1./ptt*effCorr);
 560                         }
 561
 562
 563                         track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
 564                 }
 565         }
 566
 567         if (fCaloClustersCont) {
 568                 AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
 569                 while(cluster) {
 570                         TLorentzVector nPart;
 571                         cluster->GetMomentum(nPart, fVertex);
 572                         fHistClustersPt[fCentBin]->Fill(nPart.Pt());
 573
 574                         cluster = fCaloClustersCont->GetNextAcceptCluster();
 575                 }
 576         }
 577
 578         Int_t fPttBin, fPtaBin;
 579         fPtaBin = 0;
 580
 581
 582         if (fJetsCont) {
 583                 //Int_t Njets = fJetsCont->GetNJets();
 584                 Int_t Njets = 0;
 585                 if(debug > 1){
 586                         cout << "Number of Jets: " << Njets << endl;
 587                 }
 588
 589                 //Make a array to hold jets to be tested in background jT
 590                 Float_t jetPhis[200] = {};
 591                 Float_t jetEtas[200] = {};
 592                 AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0);
 593                 Int_t ij = 0;
 594                 while(jet) {
 595                         //cout << "Jet found " << ij << " pt: " << jet->Pt() << endl;
 596                         if(jet->Pt() > 5){    //Only consider jets with pT > 5 GeV
 597                                 jetPhis[ij] = jet->Phi();
 598                                 jetEtas[ij] = jet->Eta();
 599                                 ij++;
 600                                 Njets++;
 601                                 if(debug > 1)
 602                                         cout << "i: " << ij << " jetPhi: " << jetPhis[ij] << " jetEta: " << jetEtas[ij] << endl;
 603                         }else{
 604                                 //jetPhis[ij] = 100;
 605                                 //jetEtas[ij] = 100;
 606                                 //Njets--;
 607                                 if(debug > 1)
 608                                         cout << "jetPt: " << jet->Pt() << " jetPhi: " << jet->Phi() << " jetEta: " << jet->Eta() << endl;
 609                         }
 610                         //i++;
 611                         jet = fJetsCont->GetNextAcceptJet();
 612                 }
 613
 614
 615                 //fTracks =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject( fTrackArrayName.Data() ));
 616                 jet = fJetsCont->GetNextAcceptJet(0);
 617                 while(jet) {
 618                         if(jet->Pt() > 5){
 619                                 if(jet->Eta() < -0.4 || jet->Eta() > 0.4){
 620                                         if(debug > 0)
 621                                                 cout << "Jet outside eta range, Eta: " << jet->Eta() << endl;
 622                                         jet = fJetsCont->GetNextAcceptJet();
 623                                         continue;
 624                                 }
 625                                 //cout << "Jet found " << ij << " pt: " << jet->Pt() << endl;
 626                                 //Get the trigger pT bin
 627                                 fPttBin = 0;
 628                                 for(int iptt = 0 ; iptt < fNpttBins; iptt++){
 629                                         if(jet->Pt() > TriggPtBorders[iptt]){
 630                                                 fPttBin = iptt;
 631                                         }
 632
 633                                 }
 634                                 fHistJetEta[fCentBin][fPttBin]->Fill(jet->Eta());
 635                                 if(jet->Pt() > 0 && 1.0/jet->Pt() > 0){
 636                                         fHistJetsPt[fCentBin][fPttBin]->Fill(jet->Pt(),1.0/jet->Pt());  //Fill jet dN/(pT dpT)
 637                                 }
 638
 639                                 /*if (fHistJetsCorrPtArea[fCentBin]) {
 640                                   Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
 641                                   fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
 642                                   }*/
 643                                 //Float_t jetp = sqrt(jet->Px()*jet->Px()+jet->Py()*jet->Py()+jet->Pz()*jet->Pz()); //Jet pT norm
 644
 645
 646                                 Int_t nTrack = jet->GetNumberOfTracks();
 647                                 if (debug > 0)
 648                                         cout << "Number of tracks " << nTrack << " Jet Pt: " << jet->Pt() << endl;
 649                                 fHistJetMulti[fCentBin][fPttBin]->Fill(nTrack);
 650                                 for(Int_t it = 0; it < nTrack; it++ ){
 651                                         AliVParticle *track = (AliVParticle*)jet->TrackAt( it, fTracks );
 652                                         if( !track ){
 653                                                 cout << "No Track found" << endl;
 654                                                 continue;
 655                                         }
 656                                         fPtaBin = 0; //Get the associated pT bin
 657                                         for(int ipta = 0 ; ipta < fNptaBins; ipta++){
 658                                                 if(track->Pt() > AssocPtBorders[ipta]){
 659                                                         fPtaBin = ipta;
 660                                                 }
 661                                         }
 662                                         fHistJetTracksPt[fCentBin][fPttBin]->Fill(track->Pt());
 663                                         if(debug > 2)
 664                                                 cout << "Filling fHistJetTracksPt C" << fCentBin << " T" << fPttBin << endl;
 665                                         /*Float_t dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz(); // p_track dot p_jet
 666                                         Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz()); // track pT norm
 667                                         Float_t normproduct = constp*jetp; // jet pT norm times track pT norm
 668                                         Float_t costheta2 = dotproduct/normproduct;
 669                                         //Float_t sintheta = sqrt(1-costheta2*costheta2);
 670                                         Float_t jt = constp*sqrt(1-costheta2*costheta2);*/
 671                                         Float_t jt = getJt(track,jet,0);
 672                                         double effCorr = 1./fEfficiency->GetCorrection(track->Pt(), fHadronSelectionCut, fCent);  // here you generate warning if ptt>30
 673                                         //double effCorr = 1.;
 674                                         if(jt > 0 && 1.0/jt > 0){
 675                                                 fHistTracksJt[fCentBin]->Fill(jt,1.0/jt*effCorr); //Fill dN/(djT jT)
 676                                                 fHistJTPta[fCentBin][fPttBin][fPtaBin]->Fill(jt,1.0/jt*effCorr); //Fill dN/(djT jT)
 677                                                 fHistLogJTPta[fCentBin][fPttBin][fPtaBin]->Fill(TMath::Log(jt),1.0/jt*effCorr); //Fill logarithmic dN/(djT jT)
 678                                         }
 679                                         if(debug > 1)
 680                                                 cout << "Filling JT C" << fCentBin << "T" <<  fPttBin << "A" << fPtaBin << " jt:" << jt << " with " << 1.0/jt*effCorr << endl;
 681                                 }
 682
 683                                 //Get Jet azimuth and rapidity of jet
 684                                 Float_t jetAngle = jet->Phi();
 685                                 Float_t jetRap = jet->Eta();
 686
 687                                 //Rotate jet angle for background cone
 688                                 Float_t rotatedAngle = jetAngle+TMath::Pi()/2;
 689                                 if(rotatedAngle > TMath::Pi()*2){
 690                                         rotatedAngle = rotatedAngle- TMath::Pi()*2;
 691                                 }
 692                                 Float_t jetArea = jet->Area();
 693                                 Float_t testRadius = TMath::Sqrt(jetArea/TMath::Pi());
 694
 695                                 Bool_t doBg = 1;
 696
 697                                 //Test if there are jets in the background test cone
 698                                 for(int i_j = 0; i_j < Njets; i_j++){
 699                                         Float_t diffR = TMath::Sqrt(TMath::Power(jetPhis[i_j]-rotatedAngle,2)+TMath::Power(jetEtas[i_j]-jetRap,2));
 700                                         if(debug > 1){
 701                                                 cout << "i_j: " << i_j << " JetPhi: " << jetPhis[i_j] << " jetEta: " << jetEtas[i_j] << endl;
 702                                                 cout << "DiffR: " << diffR << " doBG: " << doBg <<endl;
 703                                         }
 704                                         if(diffR < testRadius *2){ //Jets muts be at least 2*cone radius away from the background cone axis
 705                                                 doBg =0;
 706                                                 break;
 707                                         }
 708
 709                                 }
 710
 711                                 // Do jT for all particles in respect to jet axis
 712                                 if (fTracksCont) {
 713                                         int counter = 0;
 714                                         AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
 715                                         while(track) {
 716                                                 /*Float_t dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz();
 717                                                   Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
 718                                                   Float_t normproduct = constp*jetp;
 719                                                   Float_t costheta2 = dotproduct/normproduct;
 720                                                 //Float_t sintheta = sqrt(1-costheta2*costheta2);
 721                                                 Float_t jt = constp*sqrt(1-costheta2*costheta2);*/
 722                                                 Double_t jt = getJt(track,jet,0);
 723                                                 double effCorr = 1./fEfficiency->GetCorrection(track->Pt(), fHadronSelectionCut, fCent);  // here you generate warning if ptt>30
 724                                                 if(jt > 0 && 1.0/jt > 0){
 725                                                         fHistJTPta_all[fCentBin][fPttBin][fPtaBin]->Fill(jt,1.0/jt*effCorr);
 726                                                 }
 727                                                 for(int ipta = 0 ; ipta < fNptaBins; ipta++){
 728                                                         if(track->Pt() > AssocPtBorders[ipta]){
 729                                                                 fPtaBin = ipta;
 730                                                         }
 731                                                 }
 732
 733                                                 //If background is to be filled
 734                                                 if(doBg){
 735                                                         Float_t diffR = TMath::Sqrt(TMath::Power(track->Phi()-rotatedAngle,2)+TMath::Power(track->Eta()-jetRap,2));
 736                                                         //Particles in the rotated cone
 737                                                         if(diffR < testRadius){
 738                                                                 counter++;
 739                                                                 fHistBgPt[fCentBin][fPttBin]->Fill(track->Pt());
 740                                                                 /*dotproduct = -track->Px()*jet->Py()+track->Py()*jet->Px()+track->Pz()*jet->Pz();
 741                                                                   constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
 742                                                                   normproduct = constp*jetp;
 743                                                                   costheta2 = dotproduct/normproduct;
 744                                                                 //sintheta = sqrt(1-costheta2*costheta2);
 745                                                                 jt = constp*sqrt(1-costheta2*costheta2);*/
 746                                                                 jt = getJt(track,jet,1);
 747                                                                 if(jt > 0 && 1.0/jt > 0){
 748                                                                         fHistJTBg[fCentBin][fPttBin][fPtaBin]->Fill(jt,1.0/jt*effCorr);
 749                                                                         fHistLogJTBg[fCentBin][fPttBin][fPtaBin]->Fill(TMath::Log(jt),1.0/jt*effCorr);
 750                                                                 }
 751                                                                 if(debug > 1)
 752                                                                         cout << "Filling Background C" << fCentBin << "T" <<  fPttBin << "A" << fPtaBin << " jt:" << jt << " with " << 1.0/jt*effCorr << endl;
 753                                                                 //Fill background jT
 754
 755                                                         }
 756                                                 }
 757                                                 track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
 758                                         }
 759                                         if(doBg){
 760                                                 fHistBgMulti[fCentBin][fPttBin]->Fill(counter);
 761                                         }
 762                                 }
 763                                 if(doBg){
 764                                         fHistBackgroundDone[fCentBin][fPttBin]->Fill(1);
 765                                 }else{
 766                                         fHistBackgroundDone[fCentBin][fPttBin]->Fill(0);
 767                                 }
 768
 769                         }
 770                         jet = fJetsCont->GetNextAcceptJet();
 771
 772                 }
 773                 jet = fJetsCont->GetLeadingJet();
 774                 if(jet){
 775                         if(jet->Pt() > 0 && 1.0/jet->Pt() > 0){
 776                                 fHistLeadingJetPt[fCentBin]->Fill(jet->Pt(),1.0/jet->Pt());
 777                         }
 778                 }
 779         }
 780
 781         CheckClusTrackMatching();
 782
 783         return kTRUE;
 784 }
 785
 786
 787 //-----------------------------------------------------------------------
 788 Double_t AliAnalysisTaskJetJTJT::getJt(AliVTrack *track, AliEmcalJet *jet,int reverse){
 789         Float_t dotproduct = 0;
 790         Float_t jetp = sqrt(jet->Px()*jet->Px()+jet->Py()*jet->Py()+jet->Pz()*jet->Pz()); //Jet pT norm
 791         if(reverse){
 792                 dotproduct = -track->Px()*jet->Py()+track->Py()*jet->Px()+track->Pz()*jet->Pz();
 793         } else{
 794                 dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz();
 795         }
 796         Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
 797         Float_t normproduct = constp*jetp;
 798         Float_t costheta2 = dotproduct/normproduct;
 799         //Float_t sintheta = sqrt(1-costheta2*costheta2);
 800         Float_t jt = constp*sqrt(1-costheta2*costheta2);
 801         return jt;
 802 }
 803
 804 Double_t AliAnalysisTaskJetJTJT::getJt(AliVParticle *track, AliEmcalJet *jet,int reverse){
 805         Float_t dotproduct = 0;
 806         Float_t jetp = sqrt(jet->Px()*jet->Px()+jet->Py()*jet->Py()+jet->Pz()*jet->Pz()); //Jet pT norm
 807         if(reverse){
 808                 dotproduct = -track->Px()*jet->Py()+track->Py()*jet->Px()+track->Pz()*jet->Pz();
 809         } else{
 810                 dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz();
 811         }
 812         Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
 813         Float_t normproduct = constp*jetp;
 814         Float_t costheta2 = dotproduct/normproduct;
 815         //Float_t sintheta = sqrt(1-costheta2*costheta2);
 816         Float_t jt = constp*sqrt(1-costheta2*costheta2);
 817         return jt;
 818 }
 819
 820 //________________________________________________________________________
 821 void AliAnalysisTaskJetJTJT::CheckClusTrackMatching()
 822 {
 823
 824         if(!fTracksCont || !fCaloClustersCont)
 825                 return;
 826
 827         Double_t deta = 999;
 828         Double_t dphi = 999;
 829
 830         //Get closest cluster to track
 831         AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
 832         while(track) {
 833                 //Get matched cluster
 834                 Int_t emc1 = track->GetEMCALcluster();
 835                 if(fCaloClustersCont && emc1>=0) {
 836                         AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
 837                         if(clusMatch) {
 838                                 AliPicoTrack::GetEtaPhiDiff(track, clusMatch, dphi, deta);
 839                                 //fHistPtDEtaDPhiTrackClus->Fill(track->Pt(),deta,dphi);
 840                         }
 841                 }
 842                 track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
 843         }
 844
 845         //Get closest track to cluster
 846         AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
 847         while(cluster) {
 848                 TLorentzVector nPart;
 849                 cluster->GetMomentum(nPart, fVertex);
 850                 fHistClustersPt[fCentBin]->Fill(nPart.Pt());
 851
 852                 //Get matched track
 853                 AliVTrack *mt = NULL;
 854                 AliAODCaloCluster *acl = dynamic_cast<AliAODCaloCluster*>(cluster);
 855                 if(acl) {
 856                         if(acl->GetNTracksMatched()>1)
 857                                 mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
 858                 }
 859                 else {
 860                         AliESDCaloCluster *ecl = dynamic_cast<AliESDCaloCluster*>(cluster);
 861                         Int_t im = ecl->GetTrackMatchedIndex();
 862                         if(fTracksCont && im>=0) {
 863                                 mt = static_cast<AliVTrack*>(fTracksCont->GetParticle(im));
 864                         }
 865                 }
 866                 if(mt) {
 867                         AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
 868                         //fHistPtDEtaDPhiClusTrack->Fill(nPart.Pt(),deta,dphi);
 869
 870                         //debugging
 871                         /*
 872                            if(mt->IsEMCAL()) {
 873                            Int_t emc1 = mt->GetEMCALcluster();
 874                            Printf("current id: %d  emc1: %d",fCaloClustersCont->GetCurrentID(),emc1);
 875                            AliVCluster *clm = fCaloClustersCont->GetCluster(emc1);
 876                            AliPicoTrack::GetEtaPhiDiff(mt, clm, dphi, deta);
 877                            Printf("deta: %f dphi: %f",deta,dphi);
 878                            }
 879                            */
 880                 }
 881                 cluster = fCaloClustersCont->GetNextAcceptCluster();
 882         }
 883 }
 884
 885 //________________________________________________________________________
 886 /*AliJetContainer* AliAnalysisTaskJetJTJT::AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius) {
 887
 888   AliAnalysisTaskEmcalJet::ExecOnce();
 889   nJetsConts++;
 890   AliJetContainer *cont = 0x0;
 891   cont = AliAnalysisTaskEmcalJet::AddJetContainer(n,defaultCutType,jetRadius);
 892   return cont;
 893   }*/
 894
 895
 896 //________________________________________________________________________
 897 void AliAnalysisTaskJetJTJT::ExecOnce() {
 898
 899         if(debug > 0){
 900         cout << "AliAnalysisTaskJetJTJT::ExecOnce(): " << endl;
 901         cout << "Get fTracks from " << fTrackArrayName.Data() << endl;
 902         }
 903         fTracks =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject( fTrackArrayName.Data() ));
 904
 905         AliAnalysisTaskEmcalJet::ExecOnce();
 906         if(debug > 1)
 907                 cout << "Efficiency: Set Run Period Name " << runPeriod << endl;
 908         fEfficiency->SetPeriodName(runPeriod);
 909         if(debug > 1)
 910                 cout << "Efficiency: Set Run number " << InputEvent()->GetRunNumber() << endl;
 911         fEfficiency->SetRunNumber( InputEvent()->GetRunNumber() ); //TODO Get run Number
 912         if(debug > 1)
 913                 cout << "Efficiency: Load()" << endl;
 914         fEfficiency->Load();
 915         if(debug > 1)
 916                 cout << "fEfficiency loaded" << endl;
 917
 918
 919
 920         if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
 921         if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
 922         if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
 923
 924 }
 925
 926 //________________________________________________________________________
 927 Bool_t AliAnalysisTaskJetJTJT::Run()
 928 {
 929         // Run analysis code here, if needed. It will be executed before FillHistograms().
 930
 931         return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 932 }
 933
 934 //________________________________________________________________________
 935 void AliAnalysisTaskJetJTJT::Terminate(Option_t *)
 936 {
 937         // Called once at the end of the analysis.
 938 }
 939
 940
 941
 942 //________________________________________________________________________
 943 JTJTEfficiency::JTJTEfficiency():
 944         fMode(kAuto),
 945         fPeriod(-1),
 946         fDataPath(""),
 947         fName(""),
 948         fPeriodStr(""),
 949         fMCPeriodStr(""),
 950         fRunNumber(0),
 951         fTag(""),
 952         fInputRootName(""),
 953         fInputRoot(NULL),
 954         fCentBinAxis(0x0)
 955 {
 956   for (Int_t i=0; i < 3; i++) fEffDir[i] = 0;
 957
 958 }
 959
 960 JTJTEfficiency::JTJTEfficiency(const JTJTEfficiency& obj) :
 961   fMode(obj.fMode),
 962   fPeriod(obj.fPeriod),
 963   fDataPath(obj.fDataPath),
 964   fName(obj.fName),
 965   fPeriodStr(obj.fPeriodStr),
 966   fMCPeriodStr(obj.fMCPeriodStr),
 967   fRunNumber(obj.fRunNumber),
 968   fTag(obj.fTag),
 969   fInputRootName(obj.fInputRootName),
 970   fInputRoot(obj.fInputRoot),
 971   fCentBinAxis(obj.fCentBinAxis)
 972 {
 973         // copy constructor TODO: handling of pointer members
 974         JUNUSED(obj);
 975   for (Int_t i=0; i < 3; i++) fEffDir[i] = obj.fEffDir[i];
 976 }
 977
 978 JTJTEfficiency& JTJTEfficiency::operator=(const JTJTEfficiency& obj){
 979         // equal sign operator TODO: content
 980         JUNUSED(obj);
 981         return *this;
 982 }
 983
 984
 985 //________________________________________________________________________
 986 double JTJTEfficiency::GetCorrection( double pt, int icut , double cent ) const {
 987         if( fMode == kNotUse ) return 1;
 988         int icent = fCentBinAxis->FindBin( cent ) -1 ;
 989         if( icent < 0 || icent > fCentBinAxis->GetNbins()-1 ) {
 990                 cout<<"J_WARNING : Centrality "<<cent<<" is out of CentBinBorder"<<endl;
 991                 return 1;
 992         }
 993         // TODO error check for icent;
 994         int ivtx = 0;
 995         if( ! fCorrection[ivtx][icent][icut] ) {
 996                 cout<<"J_WARNING : No Eff Info "<<pt<<"\t"<<icut<<"\t"<<cent<<"\t"<<icent<<endl;
 997                 return 1;
 998         }
 999         TGraphErrors * gr = fCorrection[ivtx][icent][icut];
1000         //=== TEMPERORY SETTING. IT will be removed soon.
1001         if( pt > 30 ) pt = 30; // Getting eff of 30GeV for lager pt
1002         double cor = gr->Eval(pt);
1003         if ( cor < 0.2 ) cor = 0.2;
1004         return cor;
1005 }
1006
1007
1008 TString JTJTEfficiency::GetEffName() {
1009         /*
1010            1. kNotUse : no Load, efficiency is 1 always
1011            2. has fInputRootName : Load that or crash
1012            3. has fName : Load fName [+runnumber] or crash
1013            4. has runnumber : Find Good MC period from AliJRunTable, or crash
1014            3. has period : Find Good MC period from AliJRunTable, or crash
1015
1016            }
1017            */
1018 if(fPeriodStr == "LHC10b"){
1019         fInputRootName = "Eff--LHC10b-LHC10d1-0-.root";
1020 }
1021 if(fPeriodStr == "LHC10c"){
1022         fInputRootName = "Eff--LHC10c-LHC10d4-0-.root";
1023 }
1024 if(fPeriodStr == "LHC10d"){
1025         fInputRootName = "Eff--LHC10d-LHC10f6a-0-.root";
1026 }
1027 if(fPeriodStr == "LHC10e"){
1028         fInputRootName = "Eff--LHC10e-LHC10e20-0-.root";
1029 }
1030 if(fPeriodStr == "LHC10h"){
1031         fInputRootName = "Eff--LHC10h-LHC11a10a_bis-0-.root";
1032 }
1033 if(fPeriodStr == "LHC11a"){
1034         fInputRootName = "Eff--LHC11a-LHC11b10a-0-.root";
1035 }
1036 if(fPeriodStr == "LHC13b"){
1037         fInputRootName = "Eff--LHC13b-LHC13b2-efix_p1-0-.root";
1038 }
1039
1040 if(fPeriodStr == "LHC13c"){
1041         fInputRootName = "Eff--LHC13c-LHC13b2-efix_p1-0-.root";
1042 }
1043 if(fPeriodStr == "LHC13d"){
1044         fInputRootName = "Eff--LHC13d-LHC13b2-efix_p1-0-.root";
1045 }
1046 if(fPeriodStr == "LHC13e"){
1047         fInputRootName = "Eff--LHC13e-LHC13b2-efix_p1-0-.root";
1048 }
1049
1050 return fInputRootName;
1051 }
1052
1053 TString JTJTEfficiency::GetEffFullName() {
1054         GetEffName();
1055         fInputRootName = fDataPath + "/" + fInputRootName;
1056         return fInputRootName;
1057 }
1058
1059
1060 //________________________________________________________________________
1061 bool JTJTEfficiency::Load(){
1062         // Load Efficiency File based on fMode
1063         if( fMode == kNotUse ) {
1064                 cout<<"J_WARNING : Eff Mode is \"NOTUSE\". eff is 1 !!!"<<endl;
1065                 return true;
1066         }
1067         GetEffFullName();
1068         if (TString(fInputRootName).BeginsWith("alien:"))  TGrid::Connect("alien:");
1069         fInputRoot = TFile::Open( fInputRootName);
1070         //fInputRoot = new TFile( fInputRootName,"READ");
1071         if( !fInputRoot ) {
1072                 cout << "J_ERROR : %s does not exist" << fInputRootName << endl;
1073                 return false;
1074         }
1075
1076         //fEffDir[0] = (TDirectory*)fInputRoot->Get("EffRE");
1077         ///fEffDir[1] = (TDirectory*)fInputRoot->Get("EffMC");
1078         fEffDir[2] = (TDirectory*)fInputRoot->Get("Efficiency");
1079         //iif( fEffDir[0] && fEffDir[1] && fEffDir[2] )
1080         if( !fEffDir[2] )
1081         {
1082                 cout << "J_ERROR : Directory EFF is not exist"<<endl;
1083                 return false;
1084         }
1085
1086         fCentBinAxis = (TAxis*)fEffDir[2]->Get("CentralityBin");
1087         if( !fCentBinAxis ){
1088                 cout << "J_ERROR : No CentralityBin in directory" << endl;
1089                 gSystem->Exit(1);
1090         }
1091
1092
1093         int nVtx = 1;
1094         int nCentBin = fCentBinAxis->GetNbins();
1095         for( int ivtx=0;ivtx<nVtx;ivtx++ ){
1096                 for( int icent=0;icent<nCentBin;icent++ ){
1097                         for( int icut=0;icut<kJNTrackCuts;icut++ ){
1098                                 fCorrection[ivtx][icent][icut]
1099                                         = (TGraphErrors*) fEffDir[2]->Get(Form("gCor%02d%02d%02d", ivtx,icent,icut));
1100                                 //cout<<"J_LOG : Eff graph - "<<Form("gCor%02d%02d%02d", ivtx,icent,icut)<<" - "<<g<<endl;
1101                         }
1102                 }
1103         }
1104         cout<<"J_LOG : Eff file is "<<fInputRootName<<endl;
1105         cout<<"J_LOG : Eff Cent Bins are ";
1106         for( int i=0;i<=nCentBin;i++ ){
1107                 //cout<<fCentBinAxis->GetXbins()->At(i)<<" ";
1108         }
1109         //cout<<endl;
1110         return true;
1111 }
1112
1113
1114