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   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 190   UInt_t   ntrials  = 0;
 191   if(tree){
 192     TFile *curfile = tree->GetCurrentFile();
 193     if (!curfile) {
 194       Error("Notify","No current file");
 195       return kFALSE;
 196     }
 197
 198     TString fileName(curfile->GetName());
 199     if(fileName.Contains("AliESDs.root")){
 200         fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
 201     }
 202     else if(fileName.Contains("AliAOD.root")){
 203         fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
 204     }
 205     else if(fileName.Contains("galice.root")){
 206         // for running with galice and kinematics alone...
 207         fileName.ReplaceAll("galice.root", "pyxsec.root");
 208     }
 209     TFile *fxsec = TFile::Open(fileName.Data());
 210     if(!fxsec){
 211       Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
 212       // no a severe condition
 213       return kTRUE;
 214     }
 215     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 216     if(!xtree){
 217       Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
 218     }
 219     xtree->SetBranchAddress("xsection",&fXsection);
 220     xtree->SetBranchAddress("ntrials",&ntrials);
 221     xtree->GetEntry(0);
 222   }
 223
 224   return kTRUE;
 225 }
 226
 227
 228 //___________________________________________________________________________________________________________________________________
 229 void AliAnalysisTaskJetCorrections::UserCreateOutputObjects()
 230 {
 231   //
 232   // Create the output container
 233   //
 234   //  Printf("Analysing event  %s :: # %5d\n", gSystem->pwd(), (Int_t) fEntry);
 235
 236   if(fUseAODInput){
 237     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 238     if(!fAOD){
 239       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 240       return;
 241     }
 242   }
 243   else{
 244     //  assume that the AOD is in the general output...
 245     fAOD  = AODEvent();
 246     if(!fAOD){
 247       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 248       return;
 249     }
 250   }
 251
 252   printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
 253
 254   fList = new TList();
 255
 256   fhEGen = new TH1F("EGen", "", 100, 0, 200);
 257   fhEGen->Sumw2();
 258   fList->Add(fhEGen);
 259
 260   fhERec = new TH1F("ERec", "", 100, 0, 200);
 261   fhERec->Sumw2();
 262   fList->Add(fhERec);
 263
 264   fhEGenRest = new TH1F("EGenRest", "", 100, 0, 200);
 265   fhEGenRest->Sumw2();
 266   fList->Add(fhEGenRest);
 267
 268   fhERecRest = new TH1F("ERecRest", "", 100, 0, 200);
 269   fhERecRest->Sumw2();
 270   fList->Add(fhERecRest);
 271
 272   fhEsumGenRest = new TH1F("EsumGenRest", "", 100, 0, 200);
 273   fhEsumGenRest->Sumw2();
 274   fList->Add(fhEsumGenRest);
 275
 276   fhEsumRecRest = new TH1F("EsumRecRest", "", 100, 0, 200);
 277   fhEsumRecRest->Sumw2();
 278   fList->Add(fhEsumRecRest);
 279
 280   fhE2vsE1Gen = new TH2F("E2vsE1Gen", "", 100, 0, 200, 100, 0, 200);
 281   fhE2vsE1Gen->Sumw2();
 282   fList->Add(fhE2vsE1Gen);
 283
 284   fhE2vsE1Rec = new TH2F("E2vsE1Rec", "", 100, 0, 200, 100, 0, 200);
 285   fhE2vsE1Rec->Sumw2();
 286   fList->Add(fhE2vsE1Rec);
 287
 288   fhE2E1vsEsumGen = new TH2F("E2E1vsEsumGen", "", 100, 0, 200, 25, 0, 1);
 289   fhE2E1vsEsumGen->Sumw2();
 290   fList->Add(fhE2E1vsEsumGen);
 291
 292   fhE2E1vsEsumRec = new TH2F("E2E1vsEsumRec", "", 100, 0, 200, 25, 0, 1);
 293   fhE2E1vsEsumRec->Sumw2();
 294   fList->Add(fhE2E1vsEsumRec);
 295
 296   fhE2E1vsE1Gen = new TH2F("E2E1vsE1Gen", "", 100, 0, 200, 25, 0, 1);
 297   fhE2E1vsE1Gen->Sumw2();
 298   fList->Add(fhE2E1vsE1Gen);
 299
 300   fhE2E1vsE1Rec = new TH2F("E2E1vsE1Rec", "", 100, 0, 200, 25, 0, 1);
 301   fhE2E1vsE1Rec->Sumw2();
 302   fList->Add(fhE2E1vsE1Rec);
 303
 304   fhE2E1vsdPhiGen =  new TH2F("E2E1vsdPhiGen", "", 64, -3.20, 3.20, 25, 0, 1);
 305   fList->Add(fhE2E1vsdPhiGen);
 306
 307   fhE2E1vsdPhiRec =  new TH2F("E2E1vsdPhiRec", "", 64, -3.20, 3.20, 25, 0, 1);
 308   fList->Add(fhE2E1vsdPhiRec);
 309
 310   fhTrackBalance2 = new TH2F("TrackBalance2", "", 60, 0, 30, 60, 0, 30);
 311   fhTrackBalance2->Sumw2();
 312   fList->Add(fhTrackBalance2);
 313
 314   fhTrackBalance3 = new TH2F("TrackBalance3", "", 60, 0, 30, 60, 0, 30);
 315   fhTrackBalance3->Sumw2();
 316   fList->Add(fhTrackBalance3);
 317
 318   fhEt1Et22 = new TH2F("Et1Et22", "", 100, 0, 50, 100, 0, 50);
 319   fhEt1Et22->Sumw2();
 320   fList->Add(fhEt1Et22);
 321
 322   fhEt1Et23 = new TH2F("Et1Et23", "", 100, 0, 50, 100, 0, 50);
 323   fhEt1Et23->Sumw2();
 324   fList->Add(fhEt1Et23);
 325
 326   for(Int_t i = 0; i < 3; i++)
 327     {
 328       fhECorrJet10[i] = new TProfile(Form("ECorrJet10%d", i+1), "", 100, 0, 200, 0, 10);
 329       fhECorrJet10[i]->SetXTitle("E_{rec} [GeV]");
 330       fhECorrJet10[i]->SetYTitle("C=E_{gen}/E_{rec}");
 331       fhECorrJet10[i]->Sumw2();
 332
 333       fhECorrJet05[i] = new TProfile(Form("ECorrJet05%d", i+1), "", 100, 0, 200, 0, 10);
 334       fhECorrJet05[i]->SetXTitle("E_{rec} [GeV]");
 335       fhECorrJet05[i]->SetYTitle("C=E_{gen}/E_{rec}");
 336       fhECorrJet05[i]->Sumw2();
 337
 338       fhECorrJet01[i] = new TProfile(Form("ECorrJet01%d", i+1), "", 100, 0, 200, 0, 10);
 339       fhECorrJet01[i]->SetXTitle("E_{rec} [GeV]");
 340       fhECorrJet01[i]->SetYTitle("C=E_{gen}/E_{rec}");
 341       fhECorrJet01[i]->Sumw2();
 342
 343       fhECorrJet001[i] = new TProfile(Form("ECorrJet001%d", i+1), "", 100, 0, 200, 0, 10);
 344       fhECorrJet001[i]->SetXTitle("E_{rec} [GeV]");
 345       fhECorrJet001[i]->SetYTitle("C=E_{gen}/E_{rec}");
 346       fhECorrJet001[i]->Sumw2();
 347
 348       fhdEvsErec10[i] = new TProfile(Form("dEvsErec10_%d", i+1),"", 100, 0, 200, -1, 10);
 349       fhdEvsErec10[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
 350       fhdEvsErec10[i]->SetXTitle("E_{rec} [GeV]");
 351       fhdEvsErec10[i]->Sumw2();
 352
 353       fhdEvsErec05[i] = new TProfile(Form("dEvsErec05_%d", i+1),"", 100, 0, 200, -1, 10);
 354       fhdEvsErec05[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
 355       fhdEvsErec05[i]->SetXTitle("E_{rec} [GeV]");
 356       fhdEvsErec05[i]->Sumw2();
 357
 358       fhdEvsErec01[i] = new TProfile(Form("dEvsErec01_%d", i+1),"", 100, 0, 200, -1, 10);
 359       fhdEvsErec01[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
 360       fhdEvsErec01[i]->SetXTitle("E_{rec} [GeV]");
 361       fhdEvsErec01[i]->Sumw2();
 362
 363       fhdEvsErec001[i] = new TProfile(Form("dEvsErec001_%d", i+1),"", 100, 0, 200, -1, 10);
 364       fhdEvsErec001[i]->SetYTitle("|E_{rec}-E_{rec}|/E_{rec}");
 365       fhdEvsErec001[i]->SetXTitle("E_{rec} [GeV]");
 366       fhdEvsErec001[i]->Sumw2();
 367
 368       fhdPhidEta10[i] = new TH2F(Form("dPhidEta10_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 369       fhdPhidEta10[i]->SetXTitle("#phi [rad]");
 370       fhdPhidEta10[i]->SetYTitle("#eta");
 371       fhdPhidEta10[i]->Sumw2();
 372
 373       fhdPhidEta05[i] = new TH2F(Form("dPhidEta05_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 374       fhdPhidEta05[i]->SetXTitle("#phi [rad]");
 375       fhdPhidEta05[i]->SetYTitle("#eta");
 376       fhdPhidEta05[i]->Sumw2();
 377
 378       fhdPhidEta01[i] = new TH2F(Form("dPhidEta01_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 379       fhdPhidEta01[i]->SetXTitle("#phi [rad]");
 380       fhdPhidEta01[i]->SetYTitle("#eta");
 381       fhdPhidEta01[i]->Sumw2();
 382
 383       fhdPhidEta001[i] = new TH2F(Form("dPhidEta001_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 384       fhdPhidEta001[i]->SetXTitle("#phi [rad]");
 385       fhdPhidEta001[i]->SetYTitle("#eta");
 386       fhdPhidEta001[i]->Sumw2();
 387
 388       fhdPhidEtaPt10[i] = new TH2F(Form("dPhidEtaPt10_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 389       fhdPhidEtaPt10[i]->SetXTitle("#phi [rad]");
 390       fhdPhidEtaPt10[i]->SetYTitle("#eta");
 391       fhdPhidEtaPt10[i]->Sumw2();
 392
 393       fhdPhidEtaPt05[i] = new TH2F(Form("dPhidEtaPt05_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 394       fhdPhidEtaPt05[i]->SetXTitle("#phi [rad]");
 395       fhdPhidEtaPt05[i]->SetYTitle("#eta");
 396       fhdPhidEtaPt05[i]->Sumw2();
 397
 398       fhdPhidEtaPt01[i] = new TH2F(Form("dPhidEtaPt01_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 399       fhdPhidEtaPt01[i]->SetXTitle("#phi [rad]");
 400       fhdPhidEtaPt01[i]->SetYTitle("#eta");
 401       fhdPhidEtaPt01[i]->Sumw2();
 402
 403       fhdPhidEtaPt001[i] = new TH2F(Form("dPhidEtaPt001_%d", i+1), "", 63, (-1)*TMath::Pi(), TMath::Pi(), 18, -0.9, 0.9);
 404       fhdPhidEtaPt001[i]->SetXTitle("#phi [rad]");
 405       fhdPhidEtaPt001[i]->SetYTitle("#eta");
 406       fhdPhidEtaPt001[i]->Sumw2();
 407
 408       fList->Add(fhECorrJet10[i]);
 409       fList->Add(fhECorrJet05[i]);
 410       fList->Add(fhECorrJet01[i]);
 411       fList->Add(fhECorrJet001[i]);
 412       fList->Add(fhdEvsErec10[i]);
 413       fList->Add(fhdEvsErec05[i]);
 414       fList->Add(fhdEvsErec01[i]);
 415       fList->Add(fhdEvsErec001[i]);
 416       fList->Add(fhdPhidEta10[i]);
 417       fList->Add(fhdPhidEta05[i]);
 418       fList->Add(fhdPhidEta01[i]);
 419       fList->Add(fhdPhidEta001[i]);
 420       fList->Add(fhdPhidEtaPt10[i]);
 421       fList->Add(fhdPhidEtaPt05[i]);
 422       fList->Add(fhdPhidEtaPt01[i]);
 423       fList->Add(fhdPhidEtaPt001[i]);
 424     }
 425
 426   Printf("UserCreateOutputObjects finished\n");
 427 }
 428
 429 //__________________________________________________________________________________________________________________________________________
 430 void AliAnalysisTaskJetCorrections::Init()
 431 {
 432   printf("AliAnalysisJetCut::Init() \n");
 433 }
 434
 435 //____________________________________________________________________________________________________________________________________________
 436 void AliAnalysisTaskJetCorrections::UserExec(Option_t * )
 437 {
 438 //  if (fDebug > 1) printf("Analysing event # %5d\n", (Int_t) fEntry);
 439
 440
 441   //create an AOD handler
 442   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 443
 444   if(!aodH)
 445     {
 446       Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 447       return;
 448     }
 449
 450   AliMCEvent* mcEvent =MCEvent();
 451   if(!mcEvent){
 452     Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 453     return;
 454   }
 455
 456   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 457
 458  //primary vertex
 459   AliAODVertex * pvtx = dynamic_cast<AliAODVertex*>(fAOD->GetPrimaryVertex());
 460   if(!pvtx) return;
 461
 462  AliAODJet genJets[kMaxJets];
 463  Int_t nGenJets = 0;
 464
 465  AliAODJet recJets[kMaxJets];
 466  Int_t nRecJets = 0;
 467
 468   //array of reconstructed jets from the AOD input
 469  TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 470  if(!aodRecJets){
 471    Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 472    return;
 473  }
 474
 475  // reconstructed jets
 476  nRecJets = aodRecJets->GetEntries();
 477  nRecJets = TMath::Min(nRecJets, kMaxJets);
 478
 479  for(int ir = 0;ir < nRecJets;++ir)
 480    {
 481      AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 482      if(!tmp)continue;
 483      recJets[ir] = *tmp;
 484     }
 485
 486  // If we set a second branch for the input jets fetch this
 487  TClonesArray * aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 488
 489  if(!aodGenJets)
 490    {
 491      printf("NO MC jets branch with name %s  Found \n",fBranchGen.Data());
 492      return;
 493    }
 494
 495  //   //Generated jets
 496  nGenJets = aodGenJets->GetEntries();
 497  nGenJets = TMath::Min(nGenJets, kMaxJets);
 498
 499  for(Int_t ig =0 ; ig < nGenJets; ++ig)
 500    {
 501      AliAODJet * tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 502      if(!tmp)continue;
 503      genJets[ig] = * tmp;
 504    }
 505
 506  Double_t eRec[kMaxJets];
 507  Double_t eGen[kMaxJets];
 508
 509  Double_t eRecRest[kMaxJets];
 510  Double_t eGenRest[kMaxJets];
 511
 512 // AliAODJet jetRec[kMaxJets];
 513  AliAODJet jetGen[kMaxJets];
 514
 515  Int_t idxRec[kMaxJets];
 516  Int_t idxGen[kMaxJets];
 517
 518 // Double_t EsumRec = 0;
 519  // Double_t EsumGen =0;
 520
 521  TLorentzVector vRec[kMaxJets];
 522  TLorentzVector vGen[kMaxJets];
 523
 524  TLorentzVector vsumRec;
 525  TLorentzVector vsumGen;
 526
 527  TVector3 pRec[kMaxJets];
 528  TVector3 pGen[kMaxJets];
 529
 530  // HISTOS FOR MC
 531  Int_t nGenSel = 0;
 532  Int_t counter = 0;
 533  Int_t tag = 0;
 534
 535  AliAODJet selJets[kMaxJets];
 536
 537  // loop for applying the separation cut
 538  for(Int_t i = 0; i < nGenJets; i++)
 539    {
 540      if(nGenJets == 1)
 541        {
 542          selJets[nGenSel] = genJets[i];
 543          nGenSel++;
 544        }
 545      else
 546        {
 547          counter = 0;
 548          tag = 0;
 549          for(Int_t j = 0; j < nGenJets; j++)
 550            {
 551              if(i!=j)
 552                {
 553                  Double_t dRij = genJets[i].DeltaR(&genJets[j]);
 554                  counter++;
 555                  if(dRij > 2*fR) tag++;
 556                }
 557            }
 558          if(counter!=0)
 559            {
 560              if(tag/counter == 1)
 561                {
 562                  selJets[nGenSel] = genJets[i];
 563                  nGenSel++;
 564                }
 565            }
 566        }
 567    }
 568
 569  for (Int_t gj = 0; gj < nGenSel; gj++)
 570    {
 571      eGen[gj] = selJets[gj].E();
 572      fhEGen->Fill(eGen[gj], fXsection);
 573    }
 574
 575  TMath::Sort(nGenSel, eGen, idxGen);
 576  for (Int_t ig = 0; ig < nGenSel; ig++)
 577    jetGen[ig] = selJets[idxGen[ig]];
 578
 579  //rest frame MC jets
 580  for (Int_t i = 0; i < nGenSel; ++i)
 581    {
 582      vGen[i].SetPxPyPzE(jetGen[i].Px(), jetGen[i].Py(), jetGen[i].Pz(), jetGen[i].E());
 583      pGen[i].SetXYZ(vGen[i].Px(), vGen[i].Py(), vGen[i].Pz());
 584      vsumGen += vGen[i];
 585    }
 586
 587  if(nGenSel > 1 && pGen[0].DeltaPhi(pGen[1]) > 2.8)
 588    {
 589      fhE2vsE1Gen->Fill(jetGen[0].E(), jetGen[1].E(), fXsection);
 590      fhE2E1vsEsumGen->Fill(jetGen[0].E()+jetGen[1].E(), TMath::Abs(jetGen[0].E()-jetGen[1].E())/jetGen[0].E(), fXsection);
 591      fhE2E1vsE1Gen->Fill(jetGen[0].E(),  TMath::Abs(jetGen[0].E()-jetGen[1].E())/jetGen[0].E(), fXsection);
 592      Double_t deltaPhi = (jetGen[0].Phi()-jetGen[1].Phi());
 593      if(deltaPhi > TMath::Pi()) deltaPhi = deltaPhi - 2.*TMath::Pi();
 594      if(deltaPhi < (-1.*TMath::Pi())) deltaPhi = deltaPhi + 2.*TMath::Pi();
 595      fhE2E1vsdPhiGen->Fill(deltaPhi,  TMath::Abs(jetGen[0].E()-jetGen[1].E())/jetGen[0].E(), fXsection);
 596    }
 597
 598  Double_t fPxGen = vsumGen.Px();
 599  Double_t fPyGen = vsumGen.Py();
 600  Double_t fPzGen = vsumGen.Pz();
 601  Double_t fEGen = vsumGen.E();
 602
 603  Double_t eSumGenRest = 0;
 604  for (Int_t j = 0; j < nGenSel; j++)
 605    {
 606      vGen[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
 607      eGenRest[j] = vGen[j].E();
 608      if(nGenSel > 1)
 609        fhEGenRest->Fill(eGenRest[j], fXsection);
 610      eSumGenRest += eGenRest[j];
 611    }
 612
 613  if(nGenSel > 1)
 614    fhEsumGenRest->Fill(eSumGenRest, fXsection);
 615
 616  //END VARIABLES FOR MC JETS ---------------
 617  //   }
 618
 619  //AOD JET VARIABLES
 620  Int_t nRecSel = 0;
 621  Int_t counter1 = 0;
 622  Int_t tag1 = 0;
 623
 624  AliAODJet recSelJets[kMaxJets];
 625
 626  for(Int_t i = 0; i < nRecJets; i++)
 627    {
 628      if(nRecJets == 1)
 629        {
 630          recSelJets[nRecSel] = recJets[i];
 631          nRecSel++;
 632        }
 633      else
 634        {
 635          counter1 = 0;
 636          tag1 = 0;
 637          for(Int_t j = 0; j < nRecJets; j++)
 638            {
 639              if(i!=j)
 640                {
 641                  Double_t dRij = recJets[i].DeltaR(&recJets[j]);
 642                  counter1++;
 643                  if(dRij > 2*fR) tag1++;
 644                }
 645            }
 646          if(counter1!=0)
 647            {
 648              if(tag1/counter1 == 1)
 649                {
 650                  recSelJets[nRecSel] = recJets[i];
 651                  nRecSel++;
 652                }
 653            }
 654        }
 655    }
 656
 657  if(nRecSel == 0) return;
 658  Printf("******NUMBER OF JETS AFTER DELTA R CUT : %d **********\n", nRecSel);
 659  //sort rec/gen jets by energy in C.M.S
 660  AliAODJet jetRecTmp[kMaxJets];
 661  Int_t nAccJets = 0;
 662  Double_t jetTrackPt[kTracks];
 663  TLorentzVector jetTrackTmp[kTracks];
 664  Int_t nTracks = 0;
 665  for (Int_t rj = 0; rj < nRecSel; rj++)
 666    {
 667      TRefArray * jetTracksAOD = dynamic_cast<TRefArray*>(recSelJets[rj].GetRefTracks());
 668      if(!jetTracksAOD) continue;
 669      if(jetTracksAOD->GetEntries() < 3) continue;
 670      Int_t nJetTracks = 0;
 671      for(Int_t j = 0; j < jetTracksAOD->GetEntries(); j++)
 672        {
 673          AliAODTrack * track = dynamic_cast<AliAODTrack*>(jetTracksAOD->At(j));
 674          if(!track) continue;
 675          Double_t cv[21];
 676          track->GetCovarianceXYZPxPyPz(cv);
 677          if(cv[14] > 1000.) continue;
 678          jetTrackPt[nTracks] = track->Pt();
 679          jetTrackTmp[nTracks].SetPxPyPzE(track->Px(),track->Py(),track->Pz(),track->E());
 680          nTracks++;
 681          nJetTracks++;
 682        }
 683      if(nJetTracks < 4) continue;
 684      jetRecTmp[nAccJets] = recSelJets[rj];
 685      eRec[nAccJets] = recSelJets[rj].E();
 686      fhERec->Fill(eRec[nAccJets], fXsection);
 687      nAccJets++;
 688    }
 689
 690  if(nAccJets == 0) return;
 691  if(nTracks == 0) return;
 692
 693  Printf(" ************ Number of accepted jets : %d ************ \n", nAccJets);
 694
 695  AliAODJet jetRecAcc[kMaxJets];
 696  TMath::Sort(nAccJets, eRec, idxRec);
 697  for (Int_t rj = 0; rj < nAccJets; rj++)
 698    jetRecAcc[rj] = jetRecTmp[idxRec[rj]];
 699
 700  //rest frame for reconstructed jets
 701  for (Int_t i = 0; i < nAccJets; i++)
 702    {
 703      vRec[i].SetPxPyPzE(jetRecAcc[i].Px(), jetRecAcc[i].Py(), jetRecAcc[i].Pz(), jetRecAcc[i].E());
 704      pRec[i].SetXYZ(vRec[i].Px(), vRec[i].Py(), vRec[i].Pz());
 705      vsumRec += vRec[i];
 706    }
 707
 708  //check balance of two leading hadrons, deltaPhi > 2.
 709  Int_t idxTrack[kTracks];
 710  TMath::Sort(nTracks, jetTrackPt, idxTrack);
 711
 712  TLorentzVector jetTrack[kTracks];
 713  for(Int_t iTr = 0; iTr < nTracks; iTr++)
 714    jetTrack[iTr] = jetTrackTmp[idxTrack[iTr]];
 715
 716  Int_t n = 1;
 717  while(jetTrack[0].DeltaPhi(jetTrack[n]) < 2.8)
 718    n++;
 719
 720  Double_t et1 = 0;
 721  Double_t et2 = 0;
 722  for(Int_t iTr = 0; iTr < nTracks; iTr++)
 723    {
 724      if(TMath::Abs(jetTrack[0].DeltaPhi(jetTrack[iTr]) < 1.) && iTr != 0)
 725        et1 += jetTrack[iTr].Et();
 726
 727      if(TMath::Abs(jetTrack[n].DeltaPhi(jetTrack[iTr]) < 1.) && iTr != n)
 728        et2 += jetTrack[iTr].Et();
 729    }
 730
 731  if(nAccJets == 2)
 732    {
 733      fhTrackBalance2->Fill(jetTrack[0].Et(), jetTrack[n].Et());
 734      fhEt1Et22->Fill(et1, et2);
 735    }
 736  if(nAccJets == 3)
 737    {
 738      fhTrackBalance3->Fill(jetTrack[0].Et(), jetTrack[n].Et());
 739      fhEt1Et23->Fill(et1, et2);
 740    }
 741
 742  if(nAccJets > 1 && pRec[0].DeltaPhi(pRec[1]) > 2.8)
 743    {
 744      fhE2vsE1Rec->Fill(jetRecAcc[0].E(), jetRecAcc[1].E(), fXsection);
 745      fhE2E1vsEsumRec->Fill(jetRecAcc[0].E()+jetRecAcc[1].E(), TMath::Abs(jetRecAcc[0].E()-jetRecAcc[1].E())/jetRecAcc[0].E(), fXsection);
 746      fhE2E1vsE1Rec->Fill(jetRecAcc[0].E(), TMath::Abs(jetRecAcc[0].E()-jetRecAcc[1].E())/jetRecAcc[0].E(), fXsection);
 747      Double_t deltaPhi = (jetRecAcc[0].Phi()-jetRecAcc[1].Phi());
 748      if(deltaPhi > TMath::Pi()) deltaPhi = deltaPhi - 2.*TMath::Pi();
 749      if(deltaPhi < (-1.*TMath::Pi())) deltaPhi = deltaPhi + 2.*TMath::Pi();
 750      fhE2E1vsdPhiRec->Fill(-1*deltaPhi, TMath::Abs(jetRecAcc[0].E()-jetRecAcc[1].E())/jetRecAcc[0].E(), fXsection);
 751    }
 752
 753  Double_t fPx = vsumRec.Px();
 754  Double_t fPy = vsumRec.Py();
 755  Double_t fPz = vsumRec.Pz();
 756  Double_t fE = vsumRec.E();
 757
 758  Double_t eSumRecRest = 0;
 759  for (Int_t j = 0; j < nAccJets; j++)
 760    {
 761      vRec[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
 762      eRecRest[j] = vRec[j].E();
 763      if(nAccJets > 1)
 764        fhERecRest->Fill(eRecRest[j], fXsection);
 765      eSumRecRest += eRecRest[j];
 766    }
 767  if(nAccJets > 1)
 768    fhEsumRecRest->Fill(eSumRecRest, fXsection);
 769
 770  // Relate the jets
 771  Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 772  Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 773
 774  for(int i = 0;i<kMaxJets;++i){
 775    iGenIndex[i] = iRecIndex[i] = -1;
 776  }
 777  Double_t dR = 1.4;
 778  if(nAccJets == nGenSel)
 779    {
 780     AliAnalysisHelperJetTasks::GetClosestJets(jetGen,nGenSel,jetRecAcc,nAccJets,
 781                                                  iGenIndex,iRecIndex,0);
 782
 783      for(int ir = 0;ir < nAccJets;ir++){
 784        Int_t ig = iGenIndex[ir];
 785        if(ig>=0&&ig<nGenSel){
 786          Double_t dPhi = TMath::Abs(jetRecAcc[ir].Phi()-jetGen[ig].Phi());
 787          if(dPhi > TMath::Pi()) dPhi = dPhi - 2.*TMath::Pi();
 788          if(dPhi < (-1.*TMath::Pi())) dPhi = dPhi + 2.*TMath::Pi();
 789          Double_t sigma = TMath::Abs(jetGen[ig].E()-jetRecAcc[ir].E())/jetGen[ig].E();
 790          dR = jetRecAcc[ir].DeltaR(&jetGen[ig]);
 791          if(dR < 2*fR && dR >= fR)
 792            {
 793              switch(nAccJets)
 794                {
 795                case 1:
 796                  {
 797                    fhdEvsErec10[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 798                    fhECorrJet10[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 799                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 800                      {
 801                        fhdPhidEtaPt10[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 802                        fhdPhidEta10[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 803                      }
 804                  }
 805                  break;
 806                case 2:
 807                  {
 808                    fhdEvsErec10[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 809                    fhECorrJet10[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 810                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 811                      {
 812                        fhdPhidEtaPt10[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 813                        fhdPhidEta10[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 814                      }
 815                  }
 816                  break;
 817                case 3:
 818                  {
 819                    fhdEvsErec10[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 820                    fhECorrJet10[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 821                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 822                      {
 823                        fhdPhidEtaPt10[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 824                        fhdPhidEta10[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 825                      }
 826                  }
 827                  break;
 828                }
 829            }
 830          if(dR < fR && dR >= 0.1)
 831            {
 832              switch(nAccJets)
 833                {
 834                case 1:
 835                  {
 836                    fhdEvsErec05[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 837                    fhECorrJet05[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 838                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 839                      {
 840                        fhdPhidEtaPt05[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 841                        fhdPhidEta05[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 842                      }
 843                  }
 844                  break;
 845                case 2:
 846                  {
 847                    fhdEvsErec05[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 848                    fhECorrJet05[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 849                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 850                      {
 851                        fhdPhidEtaPt05[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 852                        fhdPhidEta05[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 853                      }
 854                  }
 855                  break;
 856                case 3:
 857                  {
 858                    fhdEvsErec05[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 859                    fhECorrJet05[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 860                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 861                      {
 862                        fhdPhidEtaPt05[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 863                        fhdPhidEta05[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 864                      }
 865                  }
 866                  break;
 867                }
 868            }
 869          if(dR < 0.1 && dR >= 0.01)
 870            {
 871              switch(nAccJets)
 872                {
 873                case 1:
 874                  {
 875                    fhdEvsErec01[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 876                    fhECorrJet01[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 877                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 878                      {
 879                        fhdPhidEtaPt01[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 880                        fhdPhidEta01[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 881                      }
 882                  }
 883                  break;
 884                case 2:
 885                  {
 886                    fhdEvsErec01[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 887                    fhECorrJet01[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 888                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 889                      {
 890                        fhdPhidEtaPt01[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 891                        fhdPhidEta01[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 892                      }
 893                  }
 894                  break;
 895                case 3:
 896                  {
 897                    fhdEvsErec01[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 898                    fhECorrJet01[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 899                    for(Int_t iTr = 0; iTr < nTracks; iTr++)
 900                      {
 901                        fhdPhidEtaPt01[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 902                        fhdPhidEta01[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 903                      }
 904                  }
 905                  break;
 906                }
 907            }
 908          if(dR > 0.01) continue;
 909          switch(nAccJets)
 910            {
 911            case 1:
 912              {
 913                fhECorrJet001[0]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 914                fhdEvsErec001[0]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 915                for(Int_t iTr = 0; iTr < nTracks; iTr++)
 916                  {
 917                    fhdPhidEtaPt001[0]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 918                    fhdPhidEta001[0]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 919                  }
 920              }
 921              break;
 922            case 2:
 923              {
 924                fhECorrJet001[1]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 925                fhdEvsErec001[1]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 926                for(Int_t iTr = 0; iTr < nTracks; iTr++)
 927                  {
 928                    fhdPhidEtaPt001[1]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 929                    fhdPhidEta001[1]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 930                  }
 931              }
 932              break;
 933            case 3:
 934              {
 935                fhECorrJet001[2]->Fill(jetRecAcc[ir].E(), jetGen[ig].E()/jetRecAcc[ir].E(), fXsection);
 936                fhdEvsErec001[2]->Fill(jetRecAcc[ir].E(), sigma, fXsection);
 937                for(Int_t iTr = 0; iTr < nTracks; iTr++)
 938                  {
 939                    fhdPhidEtaPt001[2]->Fill(jetTrack[iTr].Phi(), jetTrack[iTr].Eta(), jetTrack[iTr].Pt());
 940                    fhdPhidEta001[2]->Fill(jetTrack[iTr].Phi(), jetTrack[ir].Eta());
 941                  }
 942              }
 943              break;
 944            }
 945        }
 946      }
 947      // loop over reconstructed jets
 948    }
 949
 950  Printf("%s:%d",(char*)__FILE__,__LINE__);
 951
 952  PostData(1, fList);
 953
 954  Printf("%s:%d Data Posted",(char*)__FILE__,__LINE__);
 955
 956 }
 957
 958
 959 //__________________________________________________________________________________________________________________________________________________
 960 void AliAnalysisTaskJetCorrections::Terminate(Option_t *)
 961 {
 962   printf("AnalysisJetCorrelation::Terminate()");
 963
 964 }
 965
 966 //_______________________________________User defined functions_____________________________________________________________________________________
 967
 968 void AliAnalysisTaskJetCorrections::GetThrustAxis(TVector3 &n01, TVector3 * pTrack, const Int_t &nTracks)
 969 {
 970   //
 971   // fetch the thrust axis
 972   //
 973   TVector3 psum;
 974   Double_t psum1 = 0;
 975   Double_t psum2 = 0;
 976   Double_t thrust[kTracks];
 977   Double_t th = -3;
 978   Double_t tpom = -1;
 979   Int_t j = 0;
 980
 981
 982   //  for(Int_t j = 0; j < nTracks; j++)
 983   while(TMath::Abs(th-tpom) > 10e-7 && j < nTracks)
 984     {
 985       th = tpom;
 986       psum.SetXYZ(0., 0., 0.);
 987       psum1 = 0;
 988       psum2 = 0;
 989       for(Int_t i = 0; i < nTracks; i++)
 990         {
 991           psum1 += (TMath::Abs(n01.Dot(pTrack[i])));
 992           psum2 += pTrack[i].Mag();
 993
 994           if (n01.Dot(pTrack[i]) > 0) psum += pTrack[i];
 995           if (n01.Dot(pTrack[i]) < 0) psum -= pTrack[i];
 996         }
 997
 998       thrust[j] = psum1/psum2;
 999
1000       //      if(th == thrust[j])
1001       //        break;
1002
1003       tpom = thrust[j];
1004
1005       n01 = psum.Unit();
1006       j++;
1007     }
1008 }
1009 //______________________________________________________________________________________________________
1010
1011