PWG4/JetTasks/AliAnalysisTaskJetCorrections.cxx

   1 // **************************************
   2 // Task used for estimating a charged to neutral correction
   3 // sona.pochybova@cern.ch
   4 // *******************************************
   5
   6
   7 /**************************************************************************
   8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   9  *                                                                        *
  10  * Author: The ALICE Off-line Project.                                    *
  11  * Contributors are mentioned in the code where appropriate.              *
  12  *                                                                        *
  13  * Permission to use, copy, modify and distribute this software and its   *
  14  * documentation strictly for non-commercial purposes is hereby granted   *
  15  * without fee, provided that the above copyright notice appears in all   *
  16  * copies and that both the copyright notice and this permission notice   *
  17  * appear in the supporting documentation. The authors make no claims     *
  18  * about the suitability of this software for any purpose. It is          *
  19  * provided "as is" without express or implied warranty.                  *
  20  **************************************************************************/
  21
  22
  23 #include <TFile.h>
  24 #include <TH1F.h>
  25 #include <TH2F.h>
  26 #include <TList.h>
  27 #include <TTree.h>
  28 #include <TLorentzVector.h>
  29 #include <TClonesArray.h>
  30 #include <TRefArray.h>
  31 #include <TVector3.h>
  32 #include <TProfile.h>
  33
  34 #include "AliAnalysisTaskJetCorrections.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliAODEvent.h"
  37 #include "AliAODVertex.h"
  38 #include "AliAODHandler.h"
  39 #include "AliAODTrack.h"
  40 #include "AliAODJet.h"
  41 //#include "AliMCEvent.h"
  42
  43 #include "AliAnalysisHelperJetTasks.h"
  44
  45 //
  46 //
  47 // corrections to jet energy by sona
  48 //
  49 //
  50
  51 ClassImp(AliAnalysisTaskJetCorrections)
  52
  53   AliAnalysisTaskJetCorrections::AliAnalysisTaskJetCorrections() : AliAnalysisTaskSE(),
  54                                                                    fAOD(0x0),
  55
  56                                                                    fBranchRec(""),
  57                                                                    fBranchGen(""),
  58
  59                                                                    fUseAODInput(kFALSE),
  60
  61                                                                    fR(0x0),
  62                                                                    fList(0x0),
  63
  64                                                                    fGlobVar(1),
  65                                                                    fXsection(1),
  66
  67                                                                    fhEGen(0x0),
  68                                                                    fhERec(0x0),
  69
  70                                                                    fhEGenRest(0x0),
  71                                                                    fhERecRest(0x0),
  72
  73                                                                    fhEsumGenRest(0x0),
  74                                                                    fhEsumRecRest(0x0),
  75
  76                                                                    fhE2vsE1Gen(0x0),
  77                                                                    fhE2vsE1Rec(0x0),
  78                                                                    fhE2E1vsEsumGen(0x0),
  79                                                                    fhE2E1vsEsumRec(0x0),
  80                                                                    fhE2E1vsE1Gen(0x0),
  81                                                                    fhE2E1vsE1Rec(0x0),
  82                                                                    fhE2E1vsdPhiGen(0x0),
  83                                                                    fhE2E1vsdPhiRec(0x0),
  84
  85                                                                    fhTrackBalance2(0x0),
  86                                                                    fhTrackBalance3(0x0),
  87
  88                                                                    fhEt1Et22(0x0),
  89                                                                    fhEt1Et23(0x0)
  90
  91 {
  92   //
  93   // ctor
  94   //
  95   for (Int_t i = 0; i < 3; i++)
  96     {
  97       fhECorrJet10[i] = 0;
  98       fhECorrJet05[i] = 0;
  99       fhECorrJet01[i] = 0;
 100       fhECorrJet001[i] = 0;
 101       fhdEvsErec10[i] = 0;
 102       fhdEvsErec05[i] = 0;
 103       fhdEvsErec01[i] = 0;
 104       fhdEvsErec001[i] = 0;
 105       fhdPhidEta10[i] = 0;
 106       fhdPhidEta05[i] = 0;
 107       fhdPhidEta01[i] = 0;
 108       fhdPhidEta001[i] = 0;
 109       fhdPhidEtaPt10[i] = 0;
 110       fhdPhidEtaPt05[i] = 0;
 111       fhdPhidEtaPt01[i] = 0;
 112       fhdPhidEtaPt001[i] = 0;
 113     }
 114 }
 115
 116 AliAnalysisTaskJetCorrections::AliAnalysisTaskJetCorrections(const char * name):
 117   AliAnalysisTaskSE(name),
 118
 119   fAOD(0x0),
 120
 121   fBranchRec(""),
 122   fBranchGen(""),
 123
 124   fUseAODInput(kFALSE),
 125
 126   fR(0x0),
 127   fList(0x0),
 128
 129   fGlobVar(1),
 130   fXsection(1),
 131
 132   fhEGen(0x0),
 133   fhERec(0x0),
 134
 135   fhEGenRest(0x0),
 136   fhERecRest(0x0),
 137
 138   fhEsumGenRest(0x0),
 139   fhEsumRecRest(0x0),
 140
 141   fhE2vsE1Gen(0x0),
 142   fhE2vsE1Rec(0x0),
 143   fhE2E1vsEsumGen(0x0),
 144   fhE2E1vsEsumRec(0x0),
 145   fhE2E1vsE1Gen(0x0),
 146   fhE2E1vsE1Rec(0x0),
 147   fhE2E1vsdPhiGen(0x0),
 148   fhE2E1vsdPhiRec(0x0),
 149
 150   fhTrackBalance2(0x0),
 151   fhTrackBalance3(0x0),
 152
 153   fhEt1Et22(0x0),
 154   fhEt1Et23(0x0)
 155 {
 156   //
 157   // ctor
 158   //
 159   for (Int_t i = 0; i < 3; i++)
 160     {
 161       fhECorrJet10[i] = 0;
 162       fhECorrJet05[i] = 0;
 163       fhECorrJet01[i] = 0;
 164       fhECorrJet001[i] = 0;
 165       fhdEvsErec10[i] = 0;
 166       fhdEvsErec05[i] = 0;
 167       fhdEvsErec01[i] = 0;
 168       fhdEvsErec001[i] = 0;
 169       fhdPhidEta10[i] = 0;
 170       fhdPhidEta05[i] = 0;
 171       fhdPhidEta01[i] = 0;
 172       fhdPhidEta001[i] = 0;
 173       fhdPhidEtaPt10[i] = 0;
 174       fhdPhidEtaPt05[i] = 0;
 175       fhdPhidEtaPt01[i] = 0;
 176       fhdPhidEtaPt001[i] = 0;
 177     }
 178   DefineOutput(1, TList::Class());
 179 }
 180
 181
 182
 183 Bool_t AliAnalysisTaskJetCorrections::Notify()
 184 {
 185   //
 186   // Implemented Notify() to read the cross sections
 187   // and number of trials from pyxsec.root
 188   //
 189
 190
 191   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 192   Float_t xsection = 0;
 193   Float_t ftrials  = 1;
 194
 195   Float_t fAvgTrials = 1;
 196   if(tree){
 197     TFile *curfile = tree->GetCurrentFile();
 198     if (!curfile) {
 199       Error("Notify","No current file");
 200       return kFALSE;
 201     }
 202     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 203     // construct a poor man average trials
 204     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 205     if(ftrials>=nEntries)fAvgTrials = ftrials/nEntries; // CKB take this into account for normalisation
 206   }
 207
 208   if(xsection>0)fXsection  = xsection;
 209
 210   return kTRUE;
 211
 212 }
 213
 214
 215 //___________________________________________________________________________________________________________________________________
 216 void AliAnalysisTaskJetCorrections::UserCreateOutputObjects()
 217 {
 218   //
 219   // Create the output container
 220   //
 221   //  Printf("Analysing event  %s :: # %5d\n", gSystem->pwd(), (Int_t) fEntry);
 222
 223   if(fUseAODInput){
 224     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 225     if(!fAOD){
 226       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 227       return;
 228     }
 229   }
 230   else{
 231     //  assume that the AOD is in the general output...
 232     fAOD  = AODEvent();
 233     if(!fAOD){
 234       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 235       return;
 236     }
 237   }
 238
 239   printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
 240
 241   fList = new TList();
 242
 243   fhEGen = new TH1F("EGen", "", 100, 0, 200);
 244   fhEGen->Sumw2();
 245   fList->Add(fhEGen);
 246
 247   fhERec = new TH1F("ERec", "", 100, 0, 200);
 248   fhERec->Sumw2();
 249   fList->Add(fhERec);
 250
 251   fhEGenRest = new TH1F("EGenRest", "", 100, 0, 200);
 252   fhEGenRest->Sumw2();
 253   fList->Add(fhEGenRest);
 254
 255   fhERecRest = new TH1F("ERecRest", "", 100, 0, 200);
 256   fhERecRest->Sumw2();
 257   fList->Add(fhERecRest);
 258
 259   fhEsumGenRest = new TH1F("EsumGenRest", "", 100, 0, 200);
 260   fhEsumGenRest->Sumw2();
 261   fList->Add(fhEsumGenRest);
 262
 263   fhEsumRecRest = new TH1F("EsumRecRest", "", 100, 0, 200);
 264   fhEsumRecRest->Sumw2();
 265   fList->Add(fhEsumRecRest);
 266
 267   fhE2vsE1Gen = new TH2F("E2vsE1Gen", "", 100, 0, 200, 100, 0, 200);
 268   fhE2vsE1Gen->Sumw2();
 269   fList->Add(fhE2vsE1Gen);
 270
 271   fhE2vsE1Rec = new TH2F("E2vsE1Rec", "", 100, 0, 200, 100, 0, 200);
 272   fhE2vsE1Rec->Sumw2();
 273   fList->Add(fhE2vsE1Rec);
 274
 275   fhE2E1vsEsumGen = new TH2F("E2E1vsEsumGen", "", 100, 0, 200, 25, 0, 1);
 276   fhE2E1vsEsumGen->Sumw2();
 277   fList->Add(fhE2E1vsEsumGen);
 278
 279   fhE2E1vsEsumRec = new TH2F("E2E1vsEsumRec", "", 100, 0, 200, 25, 0, 1);
 280   fhE2E1vsEsumRec->Sumw2();
 281   fList->Add(fhE2E1vsEsumRec);
 282
 283   fhE2E1vsE1Gen = new TH2F("E2E1vsE1Gen", "", 100, 0, 200, 25, 0, 1);
 284   fhE2E1vsE1Gen->Sumw2();
 285   fList->Add(fhE2E1vsE1Gen);
 286
 287   fhE2E1vsE1Rec = new TH2F("E2E1vsE1Rec", "", 100, 0, 200, 25, 0, 1);
 288   fhE2E1vsE1Rec->Sumw2();
 289   fList->Add(fhE2E1vsE1Rec);
 290
 291   fhE2E1vsdPhiGen =  new TH2F("E2E1vsdPhiGen", "", 64, -3.20, 3.20, 25, 0, 1);
 292   fList->Add(fhE2E1vsdPhiGen);
 293
 294   fhE2E1vsdPhiRec =  new TH2F("E2E1vsdPhiRec", "", 64, -3.20, 3.20, 25, 0, 1);
 295   fList->Add(fhE2E1vsdPhiRec);
 296
 297   fhTrackBalance2 = new TH2F("TrackBalance2", "", 60, 0, 30, 60, 0, 30);
 298   fhTrackBalance2->Sumw2();
 299   fList->Add(fhTrackBalance2);
 300
 301   fhTrackBalance3 = new TH2F("TrackBalance3", "", 60, 0, 30, 60, 0, 30);
 302   fhTrackBalance3->Sumw2();
 303   fList->Add(fhTrackBalance3);
 304
 305   fhEt1Et22 = new TH2F("Et1Et22", "", 100, 0, 50, 100, 0, 50);
 306   fhEt1Et22->Sumw2();
 307   fList->Add(fhEt1Et22);
 308
 309   fhEt1Et23 = new TH2F("Et1Et23", "", 100, 0, 50, 100, 0, 50);
 310   fhEt1Et23->Sumw2();
 311   fList->Add(fhEt1Et23);
 312
 313   for(Int_t i = 0; i < 3; i++)
 314     {
 315       fhECorrJet10[i] = new TProfile(Form("ECorrJet10%d", i+1), "", 100, 0, 200, 0, 10);
 316       fhECorrJet10[i]->SetXTitle("E_{rec} [GeV]");
 317       fhECorrJet10[i]->SetYTitle("C=E_{gen}/E_{rec}");
 318       fhECorrJet10[i]->Sumw2();
 319
 320       fhECorrJet05[i] = new TProfile(Form("ECorrJet05%d", i+1), "", 100, 0, 200, 0, 10);
 321       fhECorrJet05[i]->SetXTitle("E_{rec} [GeV]");
 322       fhECorrJet05[i]->SetYTitle("C=E_{gen}/E_{rec}");
 323       fhECorrJet05[i]->Sumw2();
 324
 325       fhECorrJet01[i] = new TProfile(Form("ECorrJet01%d", i+1), "", 100, 0, 200, 0, 10);
 326       fhECorrJet01[i]->SetXTitle("E_{rec} [GeV]");
 327       fhECorrJet01[i]->SetYTitle("C=E_{gen}/E_{rec}");
 328       fhECorrJet01[i]->Sumw2();
 329
 330       fhECorrJet001[i] = new TProfile(Form("ECorrJet001%d", i+1), "", 100, 0, 200, 0, 10);
 331       fhECorrJet001[i]->SetXTitle("E_{rec} [GeV]");
 332       fhECorrJet001[i]->SetYTitle("C=E_{gen}/E_{rec}");
 333       fhECorrJet001[i]->Sumw2();
 334
 335       fhdEvsErec10[i] = new TProfile(Form("dEvsErec10_%d", i+1),"", 100, 0, 200, -1, 10);
 336       fhdEvsErec10[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
 337       fhdEvsErec10[i]->SetXTitle("E_{rec} [GeV]");
 338       fhdEvsErec10[i]->Sumw2();
 339
 340       fhdEvsErec05[i] = new TProfile(Form("dEvsErec05_%d", i+1),"", 100, 0, 200, -1, 10);
 341       fhdEvsErec05[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
 342       fhdEvsErec05[i]->SetXTitle("E_{rec} [GeV]");
 343       fhdEvsErec05[i]->Sumw2();
 344
 345       fhdEvsErec01[i] = new TProfile(Form("dEvsErec01_%d", i+1),"", 100, 0, 200, -1, 10);
 346       fhdEvsErec01[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
 347       fhdEvsErec01[i]->SetXTitle("E_{rec} [GeV]");
 348       fhdEvsErec01[i]->Sumw2();
 349
 350       fhdEvsErec001[i] = new TProfile(Form("dEvsErec001_%d", i+1),"", 100, 0, 200, -1, 10);
 351       fhdEvsErec001[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
 352       fhdEvsErec001[i]->SetXTitle("E_{rec} [GeV]");
 353       fhdEvsErec001[i]->Sumw2();
 354
 355       fhdPhidEta10[i] = new TH2F(Form("dPhidEta10_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 356       fhdPhidEta10[i]->SetXTitle("#phi [rad]");
 357       fhdPhidEta10[i]->SetYTitle("#eta");
 358       fhdPhidEta10[i]->Sumw2();
 359
 360       fhdPhidEta05[i] = new TH2F(Form("dPhidEta05_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 361       fhdPhidEta05[i]->SetXTitle("#phi [rad]");
 362       fhdPhidEta05[i]->SetYTitle("#eta");
 363       fhdPhidEta05[i]->Sumw2();
 364
 365       fhdPhidEta01[i] = new TH2F(Form("dPhidEta01_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 366       fhdPhidEta01[i]->SetXTitle("#phi [rad]");
 367       fhdPhidEta01[i]->SetYTitle("#eta");
 368       fhdPhidEta01[i]->Sumw2();
 369
 370       fhdPhidEta001[i] = new TH2F(Form("dPhidEta001_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 371       fhdPhidEta001[i]->SetXTitle("#phi [rad]");
 372       fhdPhidEta001[i]->SetYTitle("#eta");
 373       fhdPhidEta001[i]->Sumw2();
 374
 375       fhdPhidEtaPt10[i] = new TH2F(Form("dPhidEtaPt10_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 376       fhdPhidEtaPt10[i]->SetXTitle("#phi [rad]");
 377       fhdPhidEtaPt10[i]->SetYTitle("#eta");
 378       fhdPhidEtaPt10[i]->Sumw2();
 379
 380       fhdPhidEtaPt05[i] = new TH2F(Form("dPhidEtaPt05_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 381       fhdPhidEtaPt05[i]->SetXTitle("#phi [rad]");
 382       fhdPhidEtaPt05[i]->SetYTitle("#eta");
 383       fhdPhidEtaPt05[i]->Sumw2();
 384
 385       fhdPhidEtaPt01[i] = new TH2F(Form("dPhidEtaPt01_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 386       fhdPhidEtaPt01[i]->SetXTitle("#phi [rad]");
 387       fhdPhidEtaPt01[i]->SetYTitle("#eta");
 388       fhdPhidEtaPt01[i]->Sumw2();
 389
 390       fhdPhidEtaPt001[i] = new TH2F(Form("dPhidEtaPt001_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 391       fhdPhidEtaPt001[i]->SetXTitle("#phi [rad]");
 392       fhdPhidEtaPt001[i]->SetYTitle("#eta");
 393       fhdPhidEtaPt001[i]->Sumw2();
 394
 395       fList->Add(fhECorrJet10[i]);
 396       fList->Add(fhECorrJet05[i]);
 397       fList->Add(fhECorrJet01[i]);
 398       fList->Add(fhECorrJet001[i]);
 399       fList->Add(fhdEvsErec10[i]);
 400       fList->Add(fhdEvsErec05[i]);
 401       fList->Add(fhdEvsErec01[i]);
 402       fList->Add(fhdEvsErec001[i]);
 403       fList->Add(fhdPhidEta10[i]);
 404       fList->Add(fhdPhidEta05[i]);
 405       fList->Add(fhdPhidEta01[i]);
 406       fList->Add(fhdPhidEta001[i]);
 407       fList->Add(fhdPhidEtaPt10[i]);
 408       fList->Add(fhdPhidEtaPt05[i]);
 409       fList->Add(fhdPhidEtaPt01[i]);
 410       fList->Add(fhdPhidEtaPt001[i]);
 411     }
 412
 413   Printf("UserCreateOutputObjects finished\n");
 414 }
 415
 416 //__________________________________________________________________________________________________________________________________________
 417 void AliAnalysisTaskJetCorrections::Init()
 418 {
 419   printf("AliAnalysisJetCut::Init() \n");
 420 }
 421
 422 //____________________________________________________________________________________________________________________________________________
 423 void AliAnalysisTaskJetCorrections::UserExec(Option_t * )
 424 {
 425 //  if (fDebug > 1) printf("Analysing event # %5d\n", (Int_t) fEntry);
 426
 427
 428   //create an AOD handler
 429   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 430
 431   if(!aodH)
 432     {
 433       Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 434       return;
 435     }
 436
 437 //   AliMCEvent* mcEvent =MCEvent();
 438 //   if(!mcEvent){
 439 //     Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 440 //     return;
 441 //   }
 442
 443   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 444
 445  //primary vertex
 446   AliAODVertex * pvtx = dynamic_cast<AliAODVertex*>(fAOD->GetPrimaryVertex());
 447   if(!pvtx) return;
 448
 449  AliAODJet genJets[kMaxJets];
 450  Int_t nGenJets = 0;
 451
 452  AliAODJet recJets[kMaxJets];
 453  Int_t nRecJets = 0;
 454
 455   //array of reconstructed jets from the AOD input
 456  TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 457  if(!aodRecJets){
 458    Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 459    return;
 460  }
 461
 462  // reconstructed jets
 463  nRecJets = aodRecJets->GetEntries();
 464  nRecJets = TMath::Min(nRecJets, kMaxJets);
 465
 466  for(int ir = 0;ir < nRecJets;++ir)
 467    {
 468      AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 469      if(!tmp)continue;
 470      recJets[ir] = *tmp;
 471     }
 472
 473  // If we set a second branch for the input jets fetch this
 474  TClonesArray * aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 475
 476  if(!aodGenJets)
 477    {
 478      printf("NO MC jets branch with name %s  Found \n",fBranchGen.Data());
 479      return;
 480    }
 481
 482  //   //Generated jets
 483  nGenJets = aodGenJets->GetEntries();
 484  nGenJets = TMath::Min(nGenJets, kMaxJets);
 485
 486  for(Int_t ig =0 ; ig < nGenJets; ++ig)
 487    {
 488      AliAODJet * tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 489      if(!tmp)continue;
 490      genJets[ig] = * tmp;
 491    }
 492
 493  Double_t eRec[kMaxJets];
 494  Double_t eGen[kMaxJets];
 495
 496  Double_t eRecRest[kMaxJets];
 497  Double_t eGenRest[kMaxJets];
 498
 499 // AliAODJet jetRec[kMaxJets];
 500  AliAODJet jetGen[kMaxJets];
 501
 502  Int_t idxRec[kMaxJets];
 503  Int_t idxGen[kMaxJets];
 504
 505 // Double_t EsumRec = 0;
 506  // Double_t EsumGen =0;
 507
 508  TLorentzVector vRec[kMaxJets];
 509  TLorentzVector vGen[kMaxJets];
 510
 511  TLorentzVector vsumRec;
 512  TLorentzVector vsumGen;
 513
 514  TVector3 pRec[kMaxJets];
 515  TVector3 pGen[kMaxJets];
 516
 517  // HISTOS FOR MC
 518  Int_t nGenSel = 0;
 519  Int_t counter = 0;
 520  Int_t tag = 0;
 521
 522  AliAODJet selJets[kMaxJets];
 523
 524  // loop for applying the separation cut
 525  for(Int_t i = 0; i < nGenJets; i++)
 526    {
 527      if(nGenJets == 1)
 528        {
 529          selJets[nGenSel] = genJets[i];
 530          nGenSel++;
 531        }
 532      else
 533        {
 534          counter = 0;
 535          tag = 0;
 536          for(Int_t j = 0; j < nGenJets; j++)
 537            {
 538              if(i!=j)
 539                {
 540                  Double_t dRij = genJets[i].DeltaR(&genJets[j]);
 541                  counter++;
 542                  if(dRij > 2*fR) tag++;
 543                }
 544            }
 545          if(counter!=0)
 546            {
 547              if(tag/counter == 1)
 548                {
 549                  selJets[nGenSel] = genJets[i];
 550                  nGenSel++;
 551                }
 552            }
 553        }
 554    }
 555
 556  for (Int_t gj = 0; gj < nGenSel; gj++)
 557    {
 558      eGen[gj] = selJets[gj].E();
 559      fhEGen->Fill(eGen[gj], fXsection);
 560    }
 561
 562  TMath::Sort(nGenSel, eGen, idxGen);
 563  for (Int_t ig = 0; ig < nGenSel; ig++)
 564    jetGen[ig] = selJets[idxGen[ig]];
 565
 566  //rest frame MC jets
 567  for (Int_t i = 0; i < nGenSel; ++i)
 568    {
 569      vGen[i].SetPxPyPzE(jetGen[i].Px(), jetGen[i].Py(), jetGen[i].Pz(), jetGen[i].E());
 570      pGen[i].SetXYZ(vGen[i].Px(), vGen[i].Py(), vGen[i].Pz());
 571      vsumGen += vGen[i];
 572    }
 573
 574  if(nGenSel > 1 && pGen[0].DeltaPhi(pGen[1]) > 2.8)
 575    {
 576      fhE2vsE1Gen->Fill(jetGen[0].E(), jetGen[1].E(), fXsection);
 577      fhE2E1vsEsumGen->Fill(jetGen[0].E()+jetGen[1].E(), TMath::Abs(jetGen[0].E()-jetGen[1].E())/jetGen[0].E(), fXsection);
 578      fhE2E1vsE1Gen->Fill(jetGen[0].E(),  TMath::Abs(jetGen[0].E()-jetGen[1].E())/jetGen[0].E(), fXsection);
 579      Double_t deltaPhi = (jetGen[0].Phi()-jetGen[1].Phi());
 580      if(deltaPhi > TMath::Pi()) deltaPhi = deltaPhi - 2.*TMath::Pi();
 581      if(deltaPhi < (-1.*TMath::Pi())) deltaPhi = deltaPhi + 2.*TMath::Pi();
 582      fhE2E1vsdPhiGen->Fill(deltaPhi,  TMath::Abs(jetGen[0].E()-jetGen[1].E())/jetGen[0].E(), fXsection);
 583    }
 584
 585  Double_t fPxGen = vsumGen.Px();
 586  Double_t fPyGen = vsumGen.Py();
 587  Double_t fPzGen = vsumGen.Pz();
 588  Double_t fEGen = vsumGen.E();
 589
 590  Double_t eSumGenRest = 0;
 591  for (Int_t j = 0; j < nGenSel; j++)
 592    {
 593      vGen[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
 594      eGenRest[j] = vGen[j].E();
 595      if(nGenSel > 1)
 596        fhEGenRest->Fill(eGenRest[j], fXsection);
 597      eSumGenRest += eGenRest[j];
 598    }
 599
 600  if(nGenSel > 1)
 601    fhEsumGenRest->Fill(eSumGenRest, fXsection);
 602
 603  //END VARIABLES FOR MC JETS ---------------
 604  //   }
 605
 606  //AOD JET VARIABLES
 607  Int_t nRecSel = 0;
 608  Int_t counter1 = 0;
 609  Int_t tag1 = 0;
 610
 611  AliAODJet recSelJets[kMaxJets];
 612
 613  for(Int_t i = 0; i < nRecJets; i++)
 614    {
 615      if(nRecJets == 1)
 616        {
 617          recSelJets[nRecSel] = recJets[i];
 618          nRecSel++;
 619        }
 620      else
 621        {
 622          counter1 = 0;
 623          tag1 = 0;
 624          for(Int_t j = 0; j < nRecJets; j++)
 625            {
 626              if(i!=j)
 627                {
 628                  Double_t dRij = recJets[i].DeltaR(&recJets[j]);
 629                  counter1++;
 630                  if(dRij > 2*fR) tag1++;
 631                }
 632            }
 633          if(counter1!=0)
 634            {
 635              if(tag1/counter1 == 1)
 636                {
 637                  recSelJets[nRecSel] = recJets[i];
 638                  nRecSel++;
 639                }
 640            }
 641        }
 642    }
 643
 644  if(nRecSel == 0) return;
 645  Printf("******NUMBER OF JETS AFTER DELTA R CUT : %d **********\n", nRecSel);
 646  //sort rec/gen jets by energy in C.M.S
 647  AliAODJet jetRecTmp[kMaxJets];
 648  Int_t nAccJets = 0;
 649  Double_t jetTrackPt[kTracks];
 650  TLorentzVector jetTrackTmp[kTracks];
 651  Int_t nTracks = 0;
 652  for (Int_t rj = 0; rj < nRecSel; rj++)
 653    {
 654      TRefArray * jetTracksAOD = dynamic_cast<TRefArray*>(recSelJets[rj].GetRefTracks());
 655      if(!jetTracksAOD) continue;
 656      if(jetTracksAOD->GetEntries() < 3) continue;
 657      Int_t nJetTracks = 0;
 658      for(Int_t j = 0; j < jetTracksAOD->GetEntries(); j++)
 659        {
 660          AliAODTrack * track = dynamic_cast<AliAODTrack*>(jetTracksAOD->At(j));
 661          if(!track) continue;
 662          Double_t cv[21];
 663          track->GetCovarianceXYZPxPyPz(cv);
 664          if(cv[14] > 1000.) continue;
 665          jetTrackPt[nTracks] = track->Pt();
 666          jetTrackTmp[nTracks].SetPxPyPzE(track->Px(),track->Py(),track->Pz(),track->E());
 667          nTracks++;
 668          nJetTracks++;
 669        }
 670      if(nJetTracks < 4) continue;
 671      jetRecTmp[nAccJets] = recSelJets[rj];
 672      eRec[nAccJets] = recSelJets[rj].E();
 673      fhERec->Fill(eRec[nAccJets], fXsection);
 674      nAccJets++;
 675    }
 676
 677  if(nAccJets == 0) return;
 678  if(nTracks == 0) return;
 679
 680  Printf(" ************ Number of accepted jets : %d ************ \n", nAccJets);
 681
 682  AliAODJet jetRecAcc[kMaxJets];
 683  TMath::Sort(nAccJets, eRec, idxRec);
 684  for (Int_t rj = 0; rj < nAccJets; rj++)
 685    jetRecAcc[rj] = jetRecTmp[idxRec[rj]];
 686
 687  //rest frame for reconstructed jets
 688  for (Int_t i = 0; i < nAccJets; i++)
 689    {
 690      vRec[i].SetPxPyPzE(jetRecAcc[i].Px(), jetRecAcc[i].Py(), jetRecAcc[i].Pz(), jetRecAcc[i].E());
 691      pRec[i].SetXYZ(vRec[i].Px(), vRec[i].Py(), vRec[i].Pz());
 692      vsumRec += vRec[i];
 693    }
 694
 695  //check balance of two leading hadrons, deltaPhi > 2.
 696  Int_t idxTrack[kTracks];
 697  TMath::Sort(nTracks, jetTrackPt, idxTrack);
 698
 699  TLorentzVector jetTrack[kTracks];
 700  for(Int_t iTr = 0; iTr < nTracks; iTr++)
 701    jetTrack[iTr] = jetTrackTmp[idxTrack[iTr]];
 702
 703  Int_t n = 1;
 704  while(jetTrack[0].DeltaPhi(jetTrack[n]) < 2.8)
 705    n++;
 706
 707  Double_t et1 = 0;
 708  Double_t et2 = 0;
 709  for(Int_t iTr = 0; iTr < nTracks; iTr++)
 710    {
 711      if(TMath::Abs(jetTrack[0].DeltaPhi(jetTrack[iTr]) < 1.) && iTr != 0)
 712        et1 += jetTrack[iTr].Et();
 713
 714      if(TMath::Abs(jetTrack[n].DeltaPhi(jetTrack[iTr]) < 1.) && iTr != n)
 715        et2 += jetTrack[iTr].Et();
 716    }
 717
 718  if(nAccJets == 2)
 719    {
 720      fhTrackBalance2->Fill(jetTrack[0].Et(), jetTrack[n].Et());
 721      fhEt1Et22->Fill(et1, et2);
 722    }
 723  if(nAccJets == 3)
 724    {
 725      fhTrackBalance3->Fill(jetTrack[0].Et(), jetTrack[n].Et());
 726      fhEt1Et23->Fill(et1, et2);
 727    }
 728
 729  if(nAccJets > 1 && pRec[0].DeltaPhi(pRec[1]) > 2.8)
 730    {
 731      fhE2vsE1Rec->Fill(jetRecAcc[0].E(), jetRecAcc[1].E(), fXsection);
 732      fhE2E1vsEsumRec->Fill(jetRecAcc[0].E()+jetRecAcc[1].E(), TMath::Abs(jetRecAcc[0].E()-jetRecAcc[1].E())/jetRecAcc[0].E(), fXsection);
 733      fhE2E1vsE1Rec->Fill(jetRecAcc[0].E(), TMath::Abs(jetRecAcc[0].E()-jetRecAcc[1].E())/jetRecAcc[0].E(), fXsection);
 734      Double_t deltaPhi = (jetRecAcc[0].Phi()-jetRecAcc[1].Phi());
 735      if(deltaPhi > TMath::Pi()) deltaPhi = deltaPhi - 2.*TMath::Pi();
 736      if(deltaPhi < (-1.*TMath::Pi())) deltaPhi = deltaPhi + 2.*TMath::Pi();
 737      fhE2E1vsdPhiRec->Fill(-1*deltaPhi, TMath::Abs(jetRecAcc[0].E()-jetRecAcc[1].E())/jetRecAcc[0].E(), fXsection);
 738    }
 739
 740  Double_t fPx = vsumRec.Px();
 741  Double_t fPy = vsumRec.Py();
 742  Double_t fPz = vsumRec.Pz();
 743  Double_t fE = vsumRec.E();
 744
 745  Double_t eSumRecRest = 0;
 746  for (Int_t j = 0; j < nAccJets; j++)
 747    {
 748      vRec[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
 749      eRecRest[j] = vRec[j].E();
 750      if(nAccJets > 1)
 751        fhERecRest->Fill(eRecRest[j], fXsection);
 752      eSumRecRest += eRecRest[j];
 753    }
 754  if(nAccJets > 1)
 755    fhEsumRecRest->Fill(eSumRecRest, fXsection);
 756
 757  // Relate the jets
 758  Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 759  Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 760
 761  for(int i = 0;i<kMaxJets;++i){
 762    iGenIndex[i] = iRecIndex[i] = -1;
 763  }
 764  Double_t dR = 1.4;
 765  if(nAccJets == nGenSel)
 766    {
 767     AliAnalysisHelperJetTasks::GetClosestJets(jetGen,nGenSel,jetRecAcc,nAccJets,
 768                                                  iGenIndex,iRecIndex,0);
 769
 770      for(int ir = 0;ir < nAccJets;ir++){
 771        Int_t ig = iGenIndex[ir];
 772        if(ig>=0&&ig<nGenSel){
 773          Double_t dPhi = TMath::Abs(jetRecAcc[ir].Phi()-jetGen[ig].Phi());
 774          if(dPhi > TMath::Pi()) dPhi = dPhi - 2.*TMath::Pi();
 775          if(dPhi < (-1.*TMath::Pi())) dPhi = dPhi + 2.*TMath::Pi();
 776          Double_t sigma = TMath::Abs(jetGen[ig].E()-jetRecAcc[ir].E())/jetGen[ig].E();
 777          dR = jetRecAcc[ir].DeltaR(&jetGen[ig]);
 778          if(dR < 2*fR && dR >= fR)
 779            {
 780              switch(nAccJets)
 781                {
 782                case 1:
 783                  {
 784                    fhdEvsErec10[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 785                    fhECorrJet10[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 786                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 787                      {
 788                        fhdPhidEtaPt10[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 789                        fhdPhidEta10[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 790                      }
 791                  }
 792                  break;
 793                case 2:
 794                  {
 795                    fhdEvsErec10[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 796                    fhECorrJet10[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 797                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 798                      {
 799                        fhdPhidEtaPt10[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 800                        fhdPhidEta10[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 801                      }
 802                  }
 803                  break;
 804                case 3:
 805                  {
 806                    fhdEvsErec10[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 807                    fhECorrJet10[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 808                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 809                      {
 810                        fhdPhidEtaPt10[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 811                        fhdPhidEta10[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 812                      }
 813                  }
 814                  break;
 815                }
 816            }
 817          if(dR < fR && dR >= 0.1)
 818            {
 819              switch(nAccJets)
 820                {
 821                case 1:
 822                  {
 823                    fhdEvsErec05[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 824                    fhECorrJet05[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 825                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 826                      {
 827                        fhdPhidEtaPt05[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 828                        fhdPhidEta05[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 829                      }
 830                  }
 831                  break;
 832                case 2:
 833                  {
 834                    fhdEvsErec05[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 835                    fhECorrJet05[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 836                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 837                      {
 838                        fhdPhidEtaPt05[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 839                        fhdPhidEta05[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 840                      }
 841                  }
 842                  break;
 843                case 3:
 844                  {
 845                    fhdEvsErec05[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 846                    fhECorrJet05[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 847                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 848                      {
 849                        fhdPhidEtaPt05[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 850                        fhdPhidEta05[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 851                      }
 852                  }
 853                  break;
 854                }
 855            }
 856          if(dR < 0.1 && dR >= 0.01)
 857            {
 858              switch(nAccJets)
 859                {
 860                case 1:
 861                  {
 862                    fhdEvsErec01[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 863                    fhECorrJet01[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 864                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 865                      {
 866                        fhdPhidEtaPt01[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 867                        fhdPhidEta01[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 868                      }
 869                  }
 870                  break;
 871                case 2:
 872                  {
 873                    fhdEvsErec01[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 874                    fhECorrJet01[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 875                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 876                      {
 877                        fhdPhidEtaPt01[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 878                        fhdPhidEta01[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 879                      }
 880                  }
 881                  break;
 882                case 3:
 883                  {
 884                    fhdEvsErec01[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 885                    fhECorrJet01[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 886                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 887                      {
 888                        fhdPhidEtaPt01[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 889                        fhdPhidEta01[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 890                      }
 891                  }
 892                  break;
 893                }
 894            }
 895          if(dR > 0.01) continue;
 896          switch(nAccJets)
 897            {
 898            case 1:
 899              {
 900                fhECorrJet001[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 901                fhdEvsErec001[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 902                for(Int_t iTr = 0; iTr < nTracks; iTr++)
 903                  {
 904                    fhdPhidEtaPt001[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 905                    fhdPhidEta001[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 906                  }
 907              }
 908              break;
 909            case 2:
 910              {
 911                fhECorrJet001[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 912                fhdEvsErec001[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 913                for(Int_t iTr = 0; iTr < nTracks; iTr++)
 914                  {
 915                    fhdPhidEtaPt001[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 916                    fhdPhidEta001[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 917                  }
 918              }
 919              break;
 920            case 3:
 921              {
 922                fhECorrJet001[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 923                fhdEvsErec001[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 924                for(Int_t iTr = 0; iTr < nTracks; iTr++)
 925                  {
 926                    fhdPhidEtaPt001[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 927                    fhdPhidEta001[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 928                  }
 929              }
 930              break;
 931            }
 932        }
 933      }
 934      // loop over reconstructed jets
 935    }
 936
 937  Printf("%s:%d",(char*)__FILE__,__LINE__);
 938
 939  PostData(1, fList);
 940
 941  Printf("%s:%d Data Posted",(char*)__FILE__,__LINE__);
 942
 943 }
 944
 945
 946 //__________________________________________________________________________________________________________________________________________________
 947 void AliAnalysisTaskJetCorrections::Terminate(Option_t *)
 948 {
 949   printf("AnalysisJetCorrections::Terminate()");
 950
 951 }
 952
 953 //_______________________________________User defined functions_____________________________________________________________________________________
 954
 955