PWG4/JetTasks/AliAnalyseUE.cxx

   1 /*************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: A.Abrahantes, E.Lopez, S.Vallero                               *
   5  * Version 1.0                                                            *
   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 #include <TROOT.h>
  16 #include <TBranch.h>
  17 #include <TCanvas.h>
  18 #include <TChain.h>
  19 #include <TFile.h>
  20 #include <TH1F.h>
  21 #include <TH1I.h>
  22 #include <TH2F.h>
  23 #include <TList.h>
  24 #include <TLorentzVector.h>
  25 #include <TMath.h>
  26 #include <TObjArray.h>
  27 #include <TProfile.h>
  28 #include <TRandom.h>
  29 #include <TSystem.h>
  30 #include <TTree.h>
  31 #include <TVector3.h>
  32
  33 #include "AliAnalyseUE.h"
  34 #include "AliAnalysisTaskUE.h"
  35 #include "AliAnalysisTask.h"
  36 #include "AliHistogramsUE.h"
  37
  38 #include "AliAnalysisManager.h"
  39 #include "AliAODEvent.h"
  40 #include "AliESDEvent.h"
  41 #include "AliAODHandler.h"
  42 #include "AliAODInputHandler.h"
  43 #include "AliAODJet.h"
  44 #include "AliAODMCParticle.h"
  45 #include "AliAODTrack.h"
  46 #include "AliESDtrack.h"
  47 #include "AliKFVertex.h"
  48 #include "AliMCEvent.h"
  49 #include "AliMCEventHandler.h"
  50 #include "AliStack.h"
  51
  52 #include "AliAnalysisHelperJetTasks.h"
  53 #include "AliGenPythiaEventHeader.h"
  54 #include "AliInputEventHandler.h"
  55 #include "AliLog.h"
  56 #include "AliStack.h"
  57
  58 ////////////////////////////////////////////////
  59 //---------------------------------------------
  60 // Class for transverse regions analysis
  61 //---------------------------------------------
  62 ////////////////////////////////////////////////
  63
  64
  65 using namespace std;
  66
  67 ClassImp(AliAnalyseUE)
  68
  69 //-------------------------------------------------------------------
  70 AliAnalyseUE::AliAnalyseUE() :
  71   TObject(),
  72   //fTaskUE(0),
  73   fkAOD(0x0),
  74   fkMC(0x0),
  75   fkESD(0x0),
  76   fDebug(0),
  77   fSimulateChJetPt(kFALSE),
  78   fStack(0x0),
  79   fAnaType(1),
  80   fAreaReg(1.5393), // Pi*0.7*0.7
  81   fConeRadius(0.7),
  82   fFilterBit(0xFF),
  83   fRegionType(1),
  84   fUseChargeHadrons(kFALSE),
  85   fUseChPartJet(kFALSE),
  86   fUsePositiveCharge(kTRUE),
  87   fUseSingleCharge(kFALSE),
  88   fOrdering(1),
  89   fJet1EtaCut(0.2),
  90   fJet2DeltaPhiCut(2.616),    // 150 degrees
  91   fJet2RatioPtCut(0.8),
  92   fJet3PtCut(15.),
  93   fTrackEtaCut(0.9),
  94   fTrackPtCut(0.),
  95   fHistos(0x0),
  96   fSumPtRegionPosit(0.),
  97   fSumPtRegionNegat(0.),
  98   fSumPtRegionForward(0.),
  99   fSumPtRegionBackward(0.),
 100   fMaxPartPtRegion(0.),
 101   fNTrackRegionPosit(0),
 102   fNTrackRegionNegat(0),
 103   fNTrackRegionForward(0),
 104   fNTrackRegionBackward(0),
 105   fSettingsTree(0x0),
 106   fLtLabel(-999),
 107   fLtMCLabel(-999)
 108 {
 109   // constructor
 110 }
 111
 112
 113 //-------------------------------------------------------------------
 114 AliAnalyseUE::AliAnalyseUE(const AliAnalyseUE & original) :
 115   TObject(original),
 116   //fTaskUE(original.fTaskUE),
 117   fkAOD(original.fkAOD),
 118   fkMC(original.fkMC),
 119   fkESD(original.fkESD),
 120   fDebug(original.fDebug),
 121   fSimulateChJetPt(original.fSimulateChJetPt),
 122   fStack(original.fStack),
 123   fAnaType(original.fAnaType),
 124   fAreaReg(original.fAreaReg),
 125   fConeRadius(original.fConeRadius),
 126   fFilterBit(original.fFilterBit),
 127   fRegionType(original.fRegionType),
 128   fUseChargeHadrons(original.fUseChargeHadrons),
 129   fUseChPartJet(original.fUseChPartJet),
 130   fUsePositiveCharge(original.fUsePositiveCharge),
 131   fUseSingleCharge(original.fUseSingleCharge),
 132   fOrdering(original.fOrdering),
 133   fJet1EtaCut(original.fJet1EtaCut),
 134   fJet2DeltaPhiCut(original.fJet2DeltaPhiCut),
 135   fJet2RatioPtCut(original.fJet2RatioPtCut),
 136   fJet3PtCut(original.fJet3PtCut),
 137   fTrackEtaCut(original.fTrackEtaCut),
 138   fTrackPtCut(original.fTrackPtCut),
 139   fHistos(original.fHistos),
 140   fSumPtRegionPosit(original.fSumPtRegionPosit),
 141   fSumPtRegionNegat(original.fSumPtRegionNegat),
 142   fSumPtRegionForward(original.fSumPtRegionForward),
 143   fSumPtRegionBackward(original.fSumPtRegionBackward),
 144   fMaxPartPtRegion(original.fMaxPartPtRegion),
 145   fNTrackRegionPosit(original.fNTrackRegionPosit),
 146   fNTrackRegionNegat(original.fNTrackRegionNegat),
 147   fNTrackRegionForward(original.fNTrackRegionForward),
 148   fNTrackRegionBackward(original.fNTrackRegionBackward),
 149   fSettingsTree(original.fSettingsTree),
 150   fLtLabel(original.fLtLabel),
 151   fLtMCLabel(original.fLtMCLabel)
 152 {
 153   //copy constructor
 154 }
 155
 156 //-------------------------------------------------------------------
 157 AliAnalyseUE & AliAnalyseUE::operator = (const AliAnalyseUE & /*source*/)
 158 {
 159   // assignment operator
 160   return *this;
 161 }
 162
 163
 164 //-------------------------------------------------------------------
 165 AliAnalyseUE::~AliAnalyseUE(){
 166
 167   //clear memory
 168
 169
 170
 171 }
 172
 173
 174 //-------------------------------------------------------------------
 175 void AliAnalyseUE::AnalyseMC(TVector3 *jetVect,AliMCEvent *mcEvent, AliGenPythiaEventHeader *pythiaGenHeader,Int_t conePosition, Bool_t useAliStack, Bool_t constrainDistance, Double_t minDistance){
 176
 177   // Execute the analysis in case of MC input
 178   fSumPtRegionPosit = 0.;
 179   fSumPtRegionNegat = 0.;
 180   fSumPtRegionForward = 0.;
 181   fSumPtRegionBackward = 0.;
 182   fMaxPartPtRegion = 0.;
 183   fNTrackRegionPosit = 0;
 184   fNTrackRegionNegat = 0;
 185   fNTrackRegionForward = 0;
 186   fNTrackRegionBackward = 0;
 187
 188   static Double_t const  kPI     = TMath::Pi();
 189   static Double_t const  k270rad = 270.*kPI/180.;
 190
 191   //Get Jets from MC header
 192   Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 193   AliAODJet pythiaGenJets[4];
 194   TVector3 jetVectnew[4];
 195   Int_t iCount = 0;
 196   for(int ip = 0;ip < nPythiaGenJets;++ip){
 197         if (iCount>3) break;
 198         Float_t p[4];
 199         pythiaGenHeader->TriggerJet(ip,p);
 200         TVector3 tempVect(p[0],p[1],p[2]);
 201         if ( TMath::Abs(tempVect.Eta())>fJet1EtaCut ) continue;
 202         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 203         jetVectnew[iCount].SetXYZ(pythiaGenJets[iCount].Px(), pythiaGenJets[iCount].Py(), pythiaGenJets[iCount].Pz());
 204         iCount++;
 205         }
 206
 207   if (!iCount) return;// no jet in eta acceptance
 208
 209   //Search the index of the nearest MC jet to the leading jet reconstructed from the input data
 210   Int_t index = 0;
 211   if (constrainDistance){
 212         Float_t deltaR = 0.;
 213         Float_t dRTemp = 0.;
 214         for (Int_t i=0; i<iCount; i++){
 215                 if (!i) {
 216                         dRTemp = jetVectnew[i].DeltaR(jetVect[0]);
 217                         index = i;
 218                         }
 219
 220                 deltaR = jetVectnew[i].DeltaR(jetVect[0]);
 221                 if (deltaR < dRTemp){
 222                         index = i;
 223                         dRTemp = deltaR;
 224                         }
 225                 }
 226
 227         if (jetVectnew[index].DeltaR(jetVect[0]) > minDistance) return;
 228         }
 229
 230   //Let's add some taste to jet and simulate pt of charged alone
 231   //eta and phi are kept as original
 232   //Play a Normal Distribution
 233   Float_t random = 1.;
 234   if (fSimulateChJetPt){
 235         while(1){
 236                 random = gRandom->Gaus(0.6,0.25);
 237                 if (random > 0. && random < 1. &&
 238                 (random * jetVectnew[index].Pt()>6.)) break;
 239                 }
 240         }
 241
 242   //Set new Pt & Fill histogram accordingly
 243   Double_t maxPtJet1 = random * jetVectnew[index].Pt();
 244
 245   fHistos->FillHistogram("hEleadingPt", maxPtJet1 );
 246
 247   if (useAliStack){//Try Stack Information to perform UE analysis
 248
 249         AliStack* mcStack = mcEvent->Stack();//Load Stack
 250         Int_t nTracksMC = mcStack->GetNtrack();
 251         for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
 252                 //Cuts
 253                 if(!(mcStack->IsPhysicalPrimary(iTracks))) continue;
 254
 255                 TParticle* mctrk = mcStack->Particle(iTracks);
 256
 257                 Double_t charge = mctrk->GetPDG()->Charge();
 258                 Double_t pT = mctrk->Pt();
 259                 Double_t eta = mctrk->Eta();
 260                 Int_t pdgCode = mctrk->GetPdgCode();
 261
 262                 if (!TrackMCSelected(charge, pT, eta, pdgCode))continue;
 263
 264                 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
 265                 Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
 266                 if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
 267                 fHistos->FillHistogram("hdNdEtaPhiDist",deltaPhi, maxPtJet1 );
 268
 269                 fHistos->FillHistogram("hFullRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
 270
 271                 //We are not interested on stack organization but don't loose track of info
 272
 273                 TVector3 tempVector =  jetVectnew[0];
 274                 jetVectnew[0] = jetVectnew[index];
 275                 jetVectnew[index] = tempVector;
 276
 277                 Int_t region = IsTrackInsideRegion( jetVectnew, &partVect, conePosition );
 278
 279                 if (region == 1) {
 280                         if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
 281                         fSumPtRegionPosit += mctrk->Pt();
 282                         fNTrackRegionPosit++;
 283                         fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
 284                         }
 285                 if (region == -1) {
 286                         if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
 287                         fSumPtRegionNegat += mctrk->Pt();
 288                         fNTrackRegionNegat++;
 289                         fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
 290                         }
 291                 if (region == 2){ //forward
 292                         fSumPtRegionForward += mctrk->Pt();
 293                         fNTrackRegionForward++;
 294                         fHistos->FillHistogram("hRegForwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
 295                         }
 296                 if (region == -2){ //backward
 297                         fSumPtRegionBackward += mctrk->Pt();
 298                         fNTrackRegionBackward++;
 299                         fHistos->FillHistogram("hRegBackwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
 300                         }
 301                 } //end loop on stack particles
 302     }else{//Try mc Particle
 303
 304       TClonesArray* farray = (TClonesArray*)fkAOD->FindListObject("mcparticles");
 305
 306       Int_t ntrks = farray->GetEntries();
 307       if (fDebug>1) AliInfo(Form("In UE MC analysis tracks %d \n",ntrks));
 308       for (Int_t i =0 ; i < ntrks; i++){
 309         AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
 310         //Cuts
 311         if (!(mctrk->IsPhysicalPrimary())) continue;
 312         //if (!(mctrk->IsPrimary())) continue;
 313
 314         Double_t charge = mctrk->Charge();
 315         Double_t pT = mctrk->Pt();
 316         Double_t eta = mctrk->Eta();
 317         Int_t pdgCode = mctrk->GetPdgCode();
 318
 319         if (!TrackMCSelected(charge, pT, eta, pdgCode))continue;
 320
 321         TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
 322
 323         Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
 324         if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
 325         fHistos->FillHistogram("hdNdEtaPhiDist", deltaPhi, maxPtJet1 );
 326
 327         fHistos->FillHistogram("hFullRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
 328
 329         //We are not interested on stack organization but don't loose track of info
 330         TVector3 tempVector =  jetVectnew[0];
 331         jetVectnew[0] = jetVectnew[index];
 332         jetVectnew[index] = tempVector;
 333
 334         Int_t region = IsTrackInsideRegion( jetVectnew, &partVect, conePosition );
 335
 336         if (region == 1) { //right
 337           if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
 338           fSumPtRegionPosit += mctrk->Pt();
 339           fNTrackRegionPosit++;
 340           fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
 341         }
 342         if (region == -1) { //left
 343           if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
 344           fSumPtRegionNegat += mctrk->Pt();
 345           fNTrackRegionNegat++;
 346           fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
 347         }
 348         if (region == 2){ //forward
 349           fSumPtRegionForward += mctrk->Pt();
 350           fNTrackRegionForward++;
 351           fHistos->FillHistogram("hRegForwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
 352         }
 353         if (region == -2){ //backward
 354           fSumPtRegionBackward += mctrk->Pt();
 355           fNTrackRegionBackward++;
 356           fHistos->FillHistogram("hRegBackwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
 357         }
 358
 359       }//end loop AliAODMCParticle tracks
 360    }
 361 }
 362
 363
 364
 365 //-------------------------------------------------------------------
 366 Bool_t AliAnalyseUE::AnaTypeSelection(TVector3 *jetVect ){
 367
 368   // Cut events by jets topology
 369   // anaType:
 370   //     1 = inclusive,
 371   //         - Jet1 |eta| < jet1EtaCut
 372   //     2 = back to back inclusive
 373   //         - fulfill case 1
 374   //         - |Jet1.Phi - Jet2.Phi| > jet2DeltaPhiCut
 375   //         - Jet2.Pt/Jet1Pt > jet2RatioPtCut
 376   //     3 = back to back exclusive
 377   //         - fulfill case 2
 378   //         - Jet3.Pt < jet3PtCut
 379
 380   Double_t eta=jetVect[0].Eta();
 381   if( TMath::Abs(eta) > fJet1EtaCut) {
 382         if( fDebug > 1 ) AliInfo("\n   Skipping Event...Jet1 |eta| > fJet1EtaCut");
 383         return kFALSE;
 384         }
 385   // back to back inclusive
 386   if( fAnaType > 1 && fAnaType < 4 && jetVect[1].Pt() < 0. ) {
 387         if( fDebug > 1 ) AliInfo("\n   Skipping Event... no second Jet found");
 388         return kFALSE;
 389         }
 390   if( fAnaType > 1 && fAnaType < 4 && jetVect[1].Pt() > 0. ) {
 391         if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
 392         jetVect[1].Pt()/jetVect[0].Pt() < fJet2RatioPtCut ) {
 393                 if( fDebug > 1 ) AliInfo("\n   Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
 394                 return kFALSE;
 395                 }
 396         }
 397   // back to back exclusive
 398   if( fAnaType > 2 && fAnaType < 4 && jetVect[2].Pt() > 0. ) {
 399         if( jetVect[2].Pt() > fJet3PtCut ) {
 400                 if( fDebug > 1 ) AliInfo("\n   Skipping Event... Jet3.Pt > fJet3PtCut ");
 401                 return kFALSE;
 402                 }
 403         }
 404   return kTRUE;
 405 }
 406
 407
 408 //-------------------------------------------------------------------
 409 void AliAnalyseUE::FillRegions(Bool_t isNorm2Area,  TVector3 *jetVect){
 410
 411   // Fill the different topological regions
 412   Double_t maxPtJet1 = jetVect[0].Pt();
 413   static Double_t const  kPI     = TMath::Pi();
 414   static Double_t const  k120rad = 120.*kPI/180.;
 415   Double_t const kMyTolerance = 0.0000001;
 416
 417   //Area for Normalization
 418   // Forward and backward
 419   Double_t normArea = 1.;
 420   // Transverse
 421   if (isNorm2Area) {
 422         SetRegionArea(jetVect);
 423         normArea =  2.*fTrackEtaCut*k120rad ;
 424         } else fAreaReg = 1.;
 425
 426   Double_t avePosRegion = (fNTrackRegionPosit) ? fSumPtRegionPosit/fNTrackRegionPosit : 0.;
 427   Double_t aveNegRegion = (fNTrackRegionNegat) ? fSumPtRegionNegat/fNTrackRegionNegat : 0.;
 428   if( avePosRegion > aveNegRegion ) {
 429      FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg );
 430   } else {
 431      FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg );
 432   }
 433
 434   //How quantities will be sorted before Fill Min and Max Histogram
 435   //  1=Plots will be CDF-like
 436   //  2=Plots will be Marchesini-like
 437   //  3=Minimum zone is selected as the one having lowest pt per track
 438   if( fOrdering == 1 ) {
 439     if( fSumPtRegionPosit > fSumPtRegionNegat ) {
 440       FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
 441     } else {
 442       FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
 443     }
 444     if (fNTrackRegionPosit > fNTrackRegionNegat ) {
 445       FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
 446     } else {
 447       FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
 448     }
 449   } else if( fOrdering == 2 ) {
 450     if (fSumPtRegionPosit > fSumPtRegionNegat) {
 451       FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
 452       FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
 453     } else {
 454       FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
 455       FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
 456     }
 457   } else if( fOrdering == 3 ){
 458      if (avePosRegion > aveNegRegion) {
 459         FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
 460         FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
 461      }else{
 462         FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
 463         FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
 464      }
 465   }
 466   fHistos->FillHistogram("hRegionMaxPartPtMaxVsEt",maxPtJet1, fMaxPartPtRegion);
 467
 468   // Compute pedestal like magnitudes
 469   fHistos->FillHistogram("hRegionDiffSumPtVsEt",maxPtJet1, (TMath::Abs(fSumPtRegionPosit-fSumPtRegionNegat)/(2.0*fAreaReg))+kMyTolerance);
 470   fHistos->FillHistogram("hRegionAveSumPtVsEt", maxPtJet1, (fSumPtRegionPosit+fSumPtRegionNegat)/(2.0*fAreaReg));
 471
 472   // Transverse as a whole
 473   fHistos->FillHistogram("hRegTransMult", maxPtJet1, fNTrackRegionPosit + fNTrackRegionNegat, (fNTrackRegionPosit + fNTrackRegionNegat)/(2.0*fAreaReg));
 474  fHistos->FillHistogram("hRegTransSumPtVsMult",maxPtJet1, fNTrackRegionPosit + fNTrackRegionNegat , (fSumPtRegionNegat + fSumPtRegionPosit)/(2.0 *fAreaReg));
 475
 476   // Fill Histograms for Forward and away side w.r.t. leading jet direction
 477   // Pt dependence
 478   //fHistos->FillHistogram("hRegForwardSumPtVsEt",maxPtJet1, fSumPtRegionForward/normArea );
 479   //fHistos->FillHistogram("hRegForwardMultVsEt",maxPtJet1, fNTrackRegionForward/normArea );
 480   //fHistos->FillHistogram("hRegBackwardSumPtVsEt",maxPtJet1, fSumPtRegionBackward/normArea );
 481   //fHistos->FillHistogram("hRegBackwardMultVsEt",maxPtJet1, fNTrackRegionBackward/normArea);
 482
 483   // Multiplicity dependence
 484   fHistos->FillHistogram("hRegForwardMult", maxPtJet1, fNTrackRegionForward, fNTrackRegionForward/normArea);
 485   fHistos->FillHistogram("hRegForwardSumPtvsMult", maxPtJet1, fNTrackRegionForward,fSumPtRegionForward/normArea);
 486   fHistos->FillHistogram("hRegBackwardMult", maxPtJet1, fNTrackRegionBackward, fNTrackRegionBackward/normArea );
 487   fHistos->FillHistogram("hRegBackwardSumPtvsMult", maxPtJet1, fNTrackRegionBackward,fSumPtRegionBackward/normArea);
 488 }
 489
 490
 491 //-------------------------------------------------------------------
 492 void AliAnalyseUE::FindMaxMinRegions(TVector3 *jetVect, Int_t conePosition, Int_t mctrue=0, Int_t eff=0){
 493
 494   // If mctrue = 1 consider branch AliAODMCParticles
 495   // If    eff = 1 track cuts for efficiency studies
 496
 497   // Identify the different topological zones
 498   fSumPtRegionPosit = 0.;
 499   fSumPtRegionNegat = 0.;
 500   fSumPtRegionForward = 0.;
 501   fSumPtRegionBackward = 0.;
 502   fMaxPartPtRegion = 0.;
 503   fNTrackRegionPosit = 0;
 504   fNTrackRegionNegat = 0;
 505   fNTrackRegionForward = 0;
 506   fNTrackRegionBackward = 0;
 507   static Double_t const  kPI     = TMath::Pi();
 508   static Double_t const  kTWOPI  = 2.*kPI;
 509   static Double_t const  k270rad = 270.*kPI/180.;
 510   Double_t const kMyTolerance = 0.0000001;
 511
 512     Int_t nTracks=0;
 513     TClonesArray *tca = 0;
 514     if (!mctrue) {
 515         nTracks = fkAOD->GetNTracks();
 516         if (fDebug > 1) AliInfo(Form(" ==== AOD tracks = %d \n ",nTracks));
 517     }else{
 518         tca = dynamic_cast<TClonesArray*>(fkAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 519         if(!tca){
 520           Printf("%s:%d No AODMC Branch found !!!",(char*)__FILE__,__LINE__);
 521           return;
 522         }
 523         nTracks = tca->GetEntriesFast();
 524         if (fDebug > 1) AliInfo(Form(" ==== AOD MC particles = %d \n ",nTracks));
 525         }
 526
 527     //If UE task d0 distribution is not filled
 528     Int_t flag_tmp=0;
 529     if (fHistos->GetHistograms()->FindObject("hDCAxy")) flag_tmp = 1;
 530
 531     for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 532     AliVParticle *part;
 533     AliAODTrack *partRECO=0;
 534     AliAODMCParticle *partMC=0;
 535     Double_t charge=-999.;
 536     Double_t pt=-999.;
 537     Double_t eta=-999.;
 538     Int_t pdgcode=-999;
 539     TString suffix("");
 540
 541     // get track reco or true
 542     if (!mctrue){
 543          partRECO = dynamic_cast<AliAODTrack*>(fkAOD->GetTrack(ipart));
 544          part = partRECO;
 545          }
 546     else {
 547         partMC = dynamic_cast<AliAODMCParticle*>(tca->At(ipart));
 548         part = partMC;
 549         if(!partMC)return;
 550         charge = partMC->Charge();
 551         pt = partMC->Pt();
 552         eta = partMC->Eta();
 553         pdgcode = partMC->GetPdgCode();
 554         suffix.Append("MC");
 555     }
 556
 557     if(!part)return;
 558
 559     if (fDebug > 1) AliInfo(Form(" ==== AOD track = %d pt = %f charge = %d \n ",ipart,part->Pt(),part->Charge()));
 560
 561     // track selection
 562     if (!mctrue && !eff ){
 563         if( !TrackSelected(partRECO)) continue; //track selection for data and MC reconstructed
 564         if (flag_tmp){
 565                 if (fkESD && fkESD->GetNumberOfTracks() ){
 566                         AliInfo("READING ESD *************************************************");
 567                         Int_t id = partRECO->GetID();
 568                         AliESDtrack *trackESD;
 569                         trackESD = (AliESDtrack*)fkESD->GetTrack(id);
 570                         Float_t d0;
 571                         Float_t z;
 572                         trackESD->GetImpactParameters(d0,z);
 573                         fHistos->FillHistogram("hDCAxy", d0, jetVect[0].Pt());
 574                 }else AliInfo("NO TRACKS ************************************************") ;
 575         }
 576         }
 577
 578     if (!mctrue && eff){
 579         if (!TrackSelectedEfficiency(partRECO )) continue; //track selection for MC reconstructed for efficiency studies
 580                 if(fkESD && fkESD->GetNumberOfTracks()){
 581                         Int_t id = partRECO->GetID();
 582                         AliESDtrack * trackESD = (AliESDtrack*) fkESD->GetTrack(id);
 583                         Float_t d0;
 584                         Float_t z;
 585                         trackESD->GetImpactParameters(d0,z);
 586                         fHistos->FillHistogram("hDCAxyPrimary", d0, jetVect[0].Pt());
 587                         }
 588         }
 589
 590     if (mctrue){
 591          if ( !(TrackMCSelected(charge, pt, eta, pdgcode) && partMC->IsPhysicalPrimary())) continue; //track selection for MC true
 592       }
 593
 594     TVector3 partVect(part->Px(), part->Py(), part->Pz());
 595     Bool_t isFlagPart = kTRUE;
 596     Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
 597     if( deltaPhi > kTWOPI )  deltaPhi-= kTWOPI;
 598     if (fAnaType != 4 ) fHistos->FillHistogram(Form("hdNdEtaPhiDist%s",suffix.Data()),deltaPhi, jetVect[0].Pt());
 599     else if (TMath::Abs(deltaPhi-k270rad) >= kMyTolerance && TMath::Abs(jetVect[0].Eta()-partVect.Eta()) >= kMyTolerance){
 600         fHistos->FillHistogram(Form("hdNdEtaPhiDist%s",suffix.Data()),deltaPhi, jetVect[0].Pt());
 601         isFlagPart = kFALSE;
 602         }
 603
 604     fHistos->FillHistogram(Form("hFullRegPartPtDistVsEt%s",suffix.Data()), part->Pt(), jetVect[0].Pt());
 605
 606     Int_t region = IsTrackInsideRegion( jetVect, &partVect, conePosition );
 607     if (region == 1) {
 608         if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
 609         fSumPtRegionPosit += part->Pt();
 610         fNTrackRegionPosit++;
 611         fHistos->FillHistogram(Form("hTransRegPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
 612         }
 613     if (region == -1) {
 614         if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
 615         fSumPtRegionNegat += part->Pt();
 616         fNTrackRegionNegat++;
 617         fHistos->FillHistogram(Form("hTransRegPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
 618         }
 619     if (region == 2){ //forward
 620         fSumPtRegionForward += part->Pt();
 621         fNTrackRegionForward++;
 622         fHistos->FillHistogram(Form("hRegForwardPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
 623         }
 624     if (region == -2){ //backward
 625         fSumPtRegionBackward += part->Pt();
 626         fNTrackRegionBackward++;
 627         fHistos->FillHistogram(Form("hRegBackwardPartPtDistVsEt%s",suffix.Data()),part->Pt(), jetVect[0].Pt());
 628         }
 629     }//end loop AOD tracks
 630
 631 }
 632
 633 //-------------------------------------------------------------------
 634 /*TVector3 AliAnalyseUE::GetLeadingTracksMC(AliMCEvent *mcEvent){
 635
 636   fkMC = mcEvent;
 637
 638   Double_t maxPtJet1 = 0.;
 639   Int_t    index1 = -1;
 640
 641   TVector3 jetVect[3];
 642
 643   jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
 644   jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
 645   jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
 646
 647   Int_t nJets = 0;
 648
 649   TObjArray* tracks = SortChargedParticlesMC();
 650     if( tracks ) {
 651         nJets = tracks->GetEntriesFast();
 652         if( nJets > 0 ) {
 653                 index1 = 0;
 654                 TParticle* jet = (TParticle*)tracks->At(0);
 655                 maxPtJet1 = jet->Pt();
 656                 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 657                 }
 658         tracks->Clear();
 659         delete tracks;
 660   }
 661   return *jetVect;
 662
 663 }
 664 */
 665
 666 //-------------------------------------------------------------------
 667 TVector3 AliAnalyseUE::GetLeadingTracksMC(AliMCEvent *mcEvent){
 668
 669   fkMC = mcEvent;
 670
 671   Double_t maxPtJet1 = 0.;
 672   Int_t    index1 = -1;
 673
 674   TVector3 jetVect[3];
 675
 676   jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
 677   jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
 678   jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
 679
 680   Int_t nJets = 0;
 681
 682   TObjArray* tracks = SortChargedParticlesMC();
 683     if( tracks ) {
 684         nJets = tracks->GetEntriesFast();
 685         if( nJets > 0 ) {
 686                 index1 = 0;
 687                 AliAODMCParticle* jet = (AliAODMCParticle*)tracks->At(0);
 688                 fLtMCLabel = jet->GetLabel();
 689                 maxPtJet1 = jet->Pt();
 690                 jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 691                 }
 692         tracks->Clear();
 693         delete tracks;
 694   }
 695   return *jetVect;
 696
 697 }
 698
 699 //-------------------------------------------------------------------
 700 TVector3 AliAnalyseUE::GetOrderedClusters(TString aodBranch, Bool_t chargedJets, Double_t chJetPtMin){
 701
 702   // jets from AOD, on-the-fly or leading particle
 703   Double_t maxPtJet1 = 0.;
 704   Int_t    index1 = -1;
 705   Double_t maxPtJet2 = 0.;   // jet 2 need for back to back inclusive
 706   Int_t    index2 = -1;
 707   Double_t maxPtJet3 = 0.;   // jet 3 need for back to back exclusive
 708   Int_t    index3 = -1;
 709   TVector3 jetVect[3];
 710
 711   jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
 712   jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
 713   jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
 714
 715   Int_t nJets = 0;
 716   //TClonesArray* fArrayJets;
 717   TObjArray* arrayJets;
 718   // 1) JETS FROM AOD BRANCH (standard, non-standard or delta)
 719   if (!chargedJets && fAnaType != 4 ) {
 720         AliInfo(" ==== Read AODs  !");
 721         AliInfo(Form(" ====  Reading Branch: %s  ", aodBranch.Data()));
 722         arrayJets = (TObjArray*)fkAOD->GetList()->FindObject(aodBranch.Data());
 723         if (!arrayJets){
 724                AliFatal(" No jet-array! ");
 725                return *jetVect;
 726                }
 727
 728         // Find Leading Jets 1,2,3
 729         // (could be skipped if Jets are sort by Pt...)
 730         nJets=arrayJets->GetEntries();
 731         for( Int_t i=0; i<nJets; ++i ) {
 732                 AliAODJet* jet = (AliAODJet*)arrayJets->At(i);
 733                 Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!!
 734
 735                 if( jetPt > maxPtJet1 ) {
 736                         maxPtJet3 = maxPtJet2; index3 = index2;
 737                         maxPtJet2 = maxPtJet1; index2 = index1;
 738                         maxPtJet1 = jetPt; index1 = i;
 739                         } else if( jetPt > maxPtJet2 ) {
 740                         maxPtJet3 = maxPtJet2; index3 = index2;
 741                         maxPtJet2 = jetPt; index2 = i;
 742                         } else if( jetPt > maxPtJet3 ) {
 743                         maxPtJet3 = jetPt; index3 = i;
 744                         }
 745                 }
 746
 747         if( index1 != -1 ) {
 748                 AliAODJet *jet =(AliAODJet*) arrayJets->At(index1);
 749                 if(jet)jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 750                 }
 751         if( index2 != -1 ) {
 752                 AliAODJet* jet= (AliAODJet*) arrayJets->At(index2);
 753                 if(jet)jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 754                 }
 755         if( index3 != -1 ) {
 756                 AliAODJet* jet = (AliAODJet*) arrayJets->At(index3);
 757                 if(jet)jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 758                 }
 759
 760   }
 761
 762
 763   // 2) ON-THE-FLY CDF ALGORITHM
 764   if (chargedJets){
 765         arrayJets = FindChargedParticleJets(chJetPtMin);
 766         if( arrayJets ) {
 767                 nJets = arrayJets->GetEntriesFast();
 768                 if( nJets > 0 ) {
 769                         index1 = 0;
 770                         AliAODJet* jet = (AliAODJet*)arrayJets->At(0);
 771                         maxPtJet1 = jet->Pt();
 772                         jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 773                         }
 774                 if( nJets > 1 ) {
 775                         index2 = 1;
 776                         AliAODJet* jet = (AliAODJet*)arrayJets->At(1);
 777                         maxPtJet2 = jet->Pt();
 778                         jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 779                         }
 780                 if( nJets > 2 ) {
 781                         index3 = 2;
 782                         AliAODJet* jet = (AliAODJet*)arrayJets->At(2);
 783                         maxPtJet3 = jet->Pt();
 784                         jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 785                         }
 786
 787                 arrayJets->Delete();
 788                 delete arrayJets;
 789                 }
 790         }
 791
 792
 793   // 3) LEADING PARTICLE
 794   if( fAnaType == 4 ){
 795         TObjArray* tracks = SortChargedParticles();
 796         if( tracks ) {
 797                 nJets = tracks->GetEntriesFast();
 798                 if( nJets > 0 ) {
 799                         index1 = 0;
 800                         AliAODTrack* jet = (AliAODTrack*)tracks->At(0);
 801                         fLtLabel = jet->GetLabel();
 802                         maxPtJet1 = jet->Pt();
 803                         jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 804                         }
 805                 tracks->Clear();
 806                 delete tracks;
 807                 }
 808
 809         }
 810   fHistos->FillHistogram("hNJets",nJets);
 811
 812   return *jetVect;
 813
 814 }
 815
 816
 817 //-------------------------------------------------------------------
 818 void AliAnalyseUE::Initialize(AliAnalysisTaskUE& taskUE){
 819 /*void AliAnalyseUE::Initialize(AliAnalysisTask& task){// when correction task is in trunk
 820    if (task->InheritsFrom("AliAnalysisTaskUE")){
 821         AliAnalysisTaskUE *taskUE = dynamic_cast<AliAnalysisTaskUE*> task;
 822         }
 823    else if (task->InheritsFrom("AliAnalysisTaskCorrectionsUE")){
 824         AliAnalysisTaskCorrectionsUE *taskUE = dynamic_cast<AliAnalysisTaskCorrectionsUE*> task;
 825         }
 826
 827 */
 828   //Get principal settings from current instance of UE analysis-task
 829   fAnaType = taskUE.GetAnaTopology();
 830   fkAOD = taskUE.GetAOD();
 831   fConeRadius = taskUE.GetConeRadius();
 832   fDebug = taskUE.GetDebugLevel();
 833   fFilterBit = taskUE.GetFilterBit();
 834   fJet1EtaCut = taskUE.GetJet1EtaCut();
 835   fJet2DeltaPhiCut = taskUE.GetJet2DeltaPhiCut();
 836   fJet2RatioPtCut = taskUE.GetJet2RatioPtCut();
 837   fJet3PtCut = taskUE.GetJet3PtCut();
 838   fOrdering = taskUE.GetPtSumOrdering() ;
 839   fRegionType = taskUE.GetRegionType();
 840   fSimulateChJetPt = taskUE.GetSimulateChJetPt();
 841   fTrackEtaCut = taskUE.GetTrackEtaCut();
 842   fTrackPtCut = taskUE.GetTrackPtCut();
 843   fHistos = taskUE.fHistosUE;
 844   fUseChargeHadrons = taskUE.GetUseChargeHadrons();
 845   fUseChPartJet = taskUE.GetUseChPartJet();
 846   fUsePositiveCharge = taskUE.GetUseNegativeChargeType();
 847   fUseSingleCharge = taskUE.GetUseSingleCharge();
 848
 849 }
 850
 851 //-------------------------------------------------------------------
 852 void  AliAnalyseUE::Initialize(Int_t anaType, AliAODEvent* aod,Double_t coneRadius, Int_t debug, Int_t filterBit, Double_t jet1EtaCut, Double_t jet2DeltaPhiCut, Double_t jet2RatioPtCut, Double_t jet3PtCut, Int_t ordering, Int_t regionType,Bool_t simulateChJetPt, Double_t trackEtaCut, Double_t trackPtCut, Bool_t useChargeHadrons, Bool_t useChPartJet, Bool_t usePositiveCharge, Bool_t useSingleCharge, AliHistogramsUE* histos){
 853
 854    fAnaType = anaType;
 855    fkAOD = aod;
 856    fConeRadius = coneRadius;
 857    fDebug = debug;
 858    fFilterBit = filterBit;
 859    fJet1EtaCut = jet1EtaCut;
 860    fJet2DeltaPhiCut = jet2DeltaPhiCut;
 861    fJet2RatioPtCut = jet2RatioPtCut;
 862    fJet3PtCut = jet3PtCut;
 863    fOrdering = ordering;
 864    fRegionType = regionType;
 865    fSimulateChJetPt = simulateChJetPt;
 866    fTrackEtaCut = trackEtaCut;
 867    fTrackPtCut = trackPtCut;
 868    fUseChargeHadrons = useChargeHadrons;
 869    fUseChPartJet = useChPartJet;
 870    fUsePositiveCharge = usePositiveCharge;
 871    fUseSingleCharge = useSingleCharge;
 872    fHistos = histos;
 873 }
 874
 875
 876 //-------------------------------------------------------------------
 877 Bool_t  AliAnalyseUE::TriggerSelection(AliInputEventHandler* inputHandler){
 878
 879   //Use AliPhysicsSelection to select good events
 880   if( !inputHandler->InheritsFrom("AliAODInputHandler") ) { // input AOD
 881         if (inputHandler->IsEventSelected()) {
 882                 if (fDebug > 1) AliInfo(" Trigger Selection: event ACCEPTED ... ");
 883         } else {
 884                 if (fDebug > 1) AliInfo(" Trigger Selection: event REJECTED ... ");
 885                 return kFALSE;
 886                 }
 887         }
 888
 889   return kTRUE;
 890
 891 }
 892
 893
 894 //-------------------------------------------------------------------
 895 Bool_t  AliAnalyseUE::VertexSelection(AliAODEvent *aod, Int_t tracks, Double_t zed ){
 896
 897   //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
 898   Int_t nVertex = aod->GetNumberOfVertices();
 899   if( nVertex > 0 ) { // Only one vertex (reject pileup)
 900         AliAODVertex* vertex = (AliAODVertex*)aod->GetPrimaryVertex();
 901         Int_t nTracksPrim = vertex->GetNContributors();
 902         Double_t zVertex = vertex->GetZ();
 903         if (fDebug > 1) AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
 904         // Select a quality vertex by number of tracks?
 905         if( nTracksPrim < tracks || TMath::Abs(zVertex) > zed ) {
 906                 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 907                 return kFALSE;
 908                 }
 909         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
 910         } else {
 911                 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 912                 return kFALSE;
 913                 }
 914
 915   return kTRUE;
 916 }
 917
 918 //-------------------------------------------------------------------
 919 Bool_t  AliAnalyseUE::VertexSelectionOld(AliAODEvent *aod ){
 920
 921   AliKFVertex primVtx(*(aod->GetPrimaryVertex()));
 922   Int_t nTracksPrim=primVtx.GetNContributors();
 923     if (fDebug > 1) AliInfo(Form(" Primary-vertex Selection: %d",nTracksPrim));
 924     if(!nTracksPrim){
 925         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
 926         return kFALSE;
 927         }
 928     if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED ...");
 929
 930   return kTRUE;
 931 }
 932
 933 // PRIVATE METHODS **************************************************
 934
 935 TObjArray*  AliAnalyseUE::FindChargedParticleJets( Double_t chJetPtMin )
 936 {
 937   // Return a TObjArray of "charged particle jets"
 938
 939   // Charged particle jet definition from reference:
 940   // "Charged jet evolution and the underlying event
 941   //  in proton-antiproton collisions at 1.8 TeV"
 942   //  PHYSICAL REVIEW D 65 092002, CDF Collaboration
 943
 944   // We defined "jets" as circular regions in eta-phi space with
 945   // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
 946   // Our jet algorithm is as follows:
 947   //   1- Order all charged particles according to their pT .
 948   //   2- Start with the highest pT particle and include in the jet all
 949   //      particles within the radius R=0.7 considering each particle
 950   //      in the order of decreasing pT and recalculating the centroid
 951   //      of the jet after each new particle is added to the jet .
 952   //   3- Go to the next highest pT particle not already included in
 953   //      a jet and add to the jet all particles not already included in
 954   //      a jet within R=0.7.
 955   //   4- Continue until all particles are in a jet.
 956   // We defined the transverse momentum of the jet to be
 957   // the scalar pT sum of all the particles within the jet, where pT
 958   // is measured with respect to the beam axis
 959
 960   //  1 - Order all charged particles according to their pT .
 961   Int_t nTracks = fkAOD->GetNTracks();
 962   if( !nTracks ) return 0;
 963   TObjArray tracks(nTracks);
 964
 965   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 966         AliAODTrack* part = fkAOD->GetTrack( ipart );
 967         if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
 968         if( !part->Charge() ) continue;
 969         if( part->Pt() < fTrackPtCut ) continue;
 970         tracks.AddLast(part);
 971         }
 972   QSortTracks( tracks, 0, tracks.GetEntriesFast() );
 973
 974   nTracks = tracks.GetEntriesFast();
 975   if( !nTracks ) return 0;
 976
 977   TObjArray *jets = new TObjArray(nTracks);
 978   TIter itrack(&tracks);
 979   while( nTracks ) {
 980         // 2- Start with the highest pT particle ...
 981         Float_t px,py,pz,pt;
 982         AliAODTrack* track = (AliAODTrack*)itrack.Next();
 983         if( !track ) continue;
 984         px = track->Px();
 985         py = track->Py();
 986         pz = track->Pz();
 987         pt = track->Pt(); // Use the energy member to store Pt
 988         jets->AddLast( new TLorentzVector(px, py, pz, pt) );
 989         tracks.Remove( track );
 990         TLorentzVector* jet = (TLorentzVector*)jets->Last();
 991         jet->SetPtEtaPhiE( 1., jet->Eta(), jet->Phi(), pt );
 992         // 3- Go to the next highest pT particle not already included...
 993         AliAODTrack* track1;
 994         Double_t fPt = jet->E();
 995         while ( (track1  = (AliAODTrack*)(itrack.Next())) ) {
 996                 Double_t tphi = track1->Phi(); // here Phi is from 0 <-> 2Pi
 997                 if (tphi > TMath::Pi()) tphi -= 2. * TMath::Pi(); // convert to  -Pi <-> Pi
 998                 Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-tphi);
 999                 Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +dphi*dphi );
1000                 if( r < fConeRadius ) {
1001                         fPt   = jet->E()+track1->Pt();  // Scalar sum of Pt
1002                         // recalculating the centroid
1003                         Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
1004                         Double_t phi = jet->Phi()*jet->E()/fPt + tphi*track1->Pt()/fPt;
1005                         jet->SetPtEtaPhiE( 1., eta, phi, fPt );
1006                         tracks.Remove( track1 );
1007                         }
1008                 }
1009
1010         tracks.Compress();
1011         nTracks = tracks.GetEntries();
1012         //   4- Continue until all particles are in a jet.
1013         itrack.Reset();
1014         } // end while nTracks
1015
1016   // Convert to AODjets....
1017   Int_t njets = jets->GetEntriesFast();
1018   TObjArray* aodjets = new TObjArray(njets);
1019   aodjets->SetOwner(kTRUE);
1020   for(Int_t ijet=0; ijet<njets; ++ijet) {
1021         TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
1022         if (jet->E() < chJetPtMin) continue;
1023         Float_t px, py,pz,en; // convert to 4-vector
1024         px = jet->E() * TMath::Cos(jet->Phi());  // Pt * cos(phi)
1025         py = jet->E() * TMath::Sin(jet->Phi());  // Pt * sin(phi)
1026         pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
1027         en = TMath::Sqrt(px * px + py * py + pz * pz);
1028
1029         aodjets->AddLast( new AliAODJet(px, py, pz, en) );
1030         }
1031   jets->Delete();
1032   delete jets;
1033
1034   // Order jets according to their pT .
1035   QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
1036
1037   // debug
1038   if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
1039
1040   return aodjets;
1041 }
1042
1043
1044 //____________________________________________________________________
1045 void AliAnalyseUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
1046 {
1047
1048   // Fill average particle Pt of control regions
1049
1050   // Max cone
1051   fHistos->FillHistogram("hRegionAvePartPtMaxVsEt", leadingE, ptMax);
1052   // Min cone
1053   fHistos->FillHistogram("hRegionAvePartPtMinVsEt", leadingE, ptMin);
1054   // MAke distributions for UE comparison with MB data
1055   fHistos->FillHistogram("hMinRegAvePt", ptMin);
1056
1057 }
1058
1059 //____________________________________________________________________
1060 void AliAnalyseUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin  )
1061 {
1062
1063   // Fill Nch multiplicity of control regions
1064
1065   // Max cone
1066   fHistos->FillHistogram("hRegionMultMaxVsEt", leadingE, nTrackPtmax);
1067   fHistos->FillHistogram("hRegionMultMax",  nTrackPtmax);
1068
1069   // Min cone
1070   fHistos->FillHistogram("hRegionMultMinVsEt", leadingE, nTrackPtmin );
1071   fHistos->FillHistogram("hRegionMultMin", nTrackPtmin);
1072
1073   // MAke distributions for UE comparison with MB data
1074   fHistos->FillHistogram("hMinRegSumPtvsMult", nTrackPtmin,ptMin);
1075
1076 }
1077
1078 //____________________________________________________________________
1079 void AliAnalyseUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
1080 {
1081   // Fill sumPt of control regions
1082
1083   // Max cone
1084   fHistos->FillHistogram("hRegionSumPtMaxVsEt", leadingE, ptMax);
1085
1086   // Min cone
1087   fHistos->FillHistogram("hRegionSumPtMinVsEt", leadingE, ptMin);
1088
1089   // MAke distributions for UE comparison with MB data
1090   fHistos->FillHistogram("hMinRegSumPt", ptMin);
1091   fHistos->FillHistogram("hMinRegSumPtJetPtBin", leadingE, ptMin);
1092   fHistos->FillHistogram("hMaxRegSumPtJetPtBin", leadingE, ptMax);
1093
1094 }
1095
1096 //-------------------------------------------------------------------
1097 Int_t AliAnalyseUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect, Int_t conePosition)
1098 {
1099   // return de region in delta phi
1100   // -1 negative delta phi
1101   //  1 positive delta phi
1102   //  0 outside region
1103   static const Double_t k60rad  = 60.*TMath::Pi()/180.;
1104   static const Double_t k120rad = 120.*TMath::Pi()/180.;
1105
1106   Int_t region = 0;
1107   if( fRegionType == 1 ) {
1108         if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
1109         // transverse regions
1110         if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1; //left
1111         if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1;    //right
1112         if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) < k60rad ) region = 2;    //forward
1113         if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) > k120rad ) region = -2; //backward
1114
1115         } else if( fRegionType == 2 ) {
1116         // Cone regions
1117         Double_t deltaR = 0.;
1118
1119         TVector3 positVect,negatVect;
1120         if (conePosition==1){
1121                 positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
1122                 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
1123         }else if (conePosition==2){
1124                 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
1125                 positVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()+TMath::PiOver2());
1126                 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()-TMath::PiOver2());
1127         }else if (conePosition==3){
1128                 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
1129                 Double_t weightEta = jetVect[0].Eta() * jetVect[0].Pt()/(jetVect[0].Pt() + jetVect[1].Pt()) +
1130                 jetVect[1].Eta() * jetVect[1].Pt()/(jetVect[0].Pt() + jetVect[1].Pt());
1131                 //Double_t weightEta = jetVect[0].Eta() * jetVect[0].Mag()/(jetVect[0].Mag() + jetVect[1].Mag()) +
1132                 // jetVect[1].Eta() * jetVect[1].Mag()/(jetVect[0].Mag() + jetVect[1].Mag());
1133                 positVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()+TMath::PiOver2());
1134                 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()-TMath::PiOver2());
1135                 }
1136   if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) {
1137         region = 1;
1138         deltaR = positVect.DrEtaPhi(*partVect);
1139   } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) {
1140         region = -1;
1141         deltaR = negatVect.DrEtaPhi(*partVect);
1142         }
1143
1144   if (deltaR > fConeRadius) region = 0;
1145
1146   } else
1147         AliError("Unknow region type");
1148
1149         return region;
1150   }
1151
1152
1153
1154 //-------------------------------------------------------------------
1155 void  AliAnalyseUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
1156 {
1157   // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
1158
1159   static TObject *tmp;
1160   static int i;           // "static" to save stack space
1161   int j;
1162
1163   while (last - first > 1) {
1164     i = first;
1165     j = last;
1166     for (;;) {
1167       while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
1168         ;
1169       while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
1170         ;
1171       if (i >= j)
1172         break;
1173
1174       tmp  = a[i];
1175       a[i] = a[j];
1176       a[j] = tmp;
1177     }
1178     if (j == first) {
1179       ++first;
1180       continue;
1181     }
1182     tmp = a[first];
1183     a[first] = a[j];
1184     a[j] = tmp;
1185     if (j - first < last - (j + 1)) {
1186       QSortTracks(a, first, j);
1187       first = j + 1;   // QSortTracks(j + 1, last);
1188     } else {
1189       QSortTracks(a, j + 1, last);
1190       last = j;        // QSortTracks(first, j);
1191     }
1192   }
1193 }
1194
1195
1196 //-------------------------------------------------------------------
1197 void  AliAnalyseUE::QSortTracksMC(TObjArray &a, Int_t first, Int_t last)
1198 {
1199   // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
1200
1201   static TObject *tmp;
1202   static int i;           // "static" to save stack space
1203   int j;
1204
1205   while (last - first > 1) {
1206     i = first;
1207     j = last;
1208     for (;;) {
1209       while (++i < last && ((AliAODMCParticle*)a[i])->Pt() > ((AliAODMCParticle*)a[first])->Pt() )
1210         ;
1211       while (--j > first && ((AliAODMCParticle*)a[j])->Pt() < ((AliAODMCParticle*)a[first])->Pt() )
1212         ;
1213       if (i >= j)
1214         break;
1215
1216       tmp  = a[i];
1217       a[i] = a[j];
1218       a[j] = tmp;
1219     }
1220     if (j == first) {
1221       ++first;
1222       continue;
1223     }
1224     tmp = a[first];
1225     a[first] = a[j];
1226     a[j] = tmp;
1227     if (j - first < last - (j + 1)) {
1228       QSortTracksMC(a, first, j);
1229       first = j + 1;   // QSortTracks(j + 1, last);
1230     } else {
1231       QSortTracksMC(a, j + 1, last);
1232       last = j;        // QSortTracks(first, j);
1233     }
1234   }
1235 }
1236
1237
1238
1239
1240
1241 //-------------------------------------------------------------------
1242 void AliAnalyseUE::SetRegionArea(TVector3 *jetVect)
1243 {
1244   // Set region area
1245   Double_t areaCorrFactor=0.;
1246   Double_t deltaEta = 0.;
1247   if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.;
1248   else if (fRegionType==2){
1249         deltaEta = 0.9-TMath::Abs(jetVect[0].Eta());
1250         if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius;
1251         else{
1252                 areaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
1253                 (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
1254                 fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-areaCorrFactor;
1255                 }
1256         }else AliWarning("Unknown Region Type");
1257   if (fDebug>10) AliInfo(Form("\n dEta=%5.3f Angle =%5.3f Region Area = %5.3f Corr Factor=%5.4f \n",deltaEta,TMath::ACos(deltaEta/fConeRadius),fAreaReg,areaCorrFactor));
1258 }
1259
1260
1261 //____________________________________________________________________
1262 TObjArray*  AliAnalyseUE::SortChargedParticles()
1263 {
1264   //  return an array with all charged particles ordered according to their pT .
1265   Int_t nTracks = fkAOD->GetNTracks();
1266   if( !nTracks ) return 0;
1267   TObjArray* tracks = new TObjArray(nTracks);
1268
1269   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1270         AliAODTrack* part = fkAOD->GetTrack( ipart );
1271         if( !part->TestFilterBit( fFilterBit ) ) continue; // track cut selection
1272         if( !part->Charge() ) continue;
1273         if( part->Pt() < fTrackPtCut ) continue;
1274         tracks->AddLast( part );
1275         }
1276   QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
1277
1278   nTracks = tracks->GetEntriesFast();
1279   if( !nTracks ) return 0;
1280
1281   return tracks;
1282   }
1283
1284 //____________________________________________________________________
1285 /*TObjArray*  AliAnalyseUE::SortChargedParticlesMC()
1286 {
1287   //  return an array with all charged particles ordered according to their pT .
1288   AliStack *mcStack = fkMC->Stack();
1289   Int_t nTracks = mcStack->GetNtrack();
1290   if( !nTracks ) return 0;
1291   TObjArray* tracks = new TObjArray(nTracks);
1292
1293   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1294         if (!(mcStack->IsPhysicalPrimary(ipart))) continue;
1295         TParticle *part = mcStack->Particle(ipart);
1296
1297         if( !part->GetPDG()->Charge() ) continue;
1298         if( part->Pt() < fTrackPtCut ) continue;
1299         tracks->AddLast( part );
1300         }
1301   QSortTracksMC( *tracks, 0, tracks->GetEntriesFast() );
1302
1303   nTracks = tracks->GetEntriesFast();
1304   if( !nTracks ) return 0;
1305
1306   return tracks;
1307   }
1308 */
1309
1310 //____________________________________________________________________
1311 TObjArray*  AliAnalyseUE::SortChargedParticlesMC()
1312 {
1313   //  return an array with all charged particles ordered according to their pT .
1314   TClonesArray *tca = dynamic_cast<TClonesArray*>(fkAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1315   if(!tca){
1316         Printf("No branch!!!");
1317         return 0;
1318   }
1319   Int_t nTracks = tca->GetEntriesFast();
1320   if( !nTracks ) return 0;
1321   TObjArray* tracks = new TObjArray(nTracks);
1322
1323   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1324         AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(ipart));
1325         if(!part)continue;
1326         if(!part->IsPhysicalPrimary())continue;
1327         if( part->Pt() < fTrackPtCut ) continue;
1328         if(!part->Charge()) continue;
1329         tracks->AddLast( part );
1330         }
1331   QSortTracksMC( *tracks, 0, tracks->GetEntriesFast() );
1332
1333   nTracks = tracks->GetEntriesFast();
1334   if( !nTracks ) return 0;
1335
1336   return tracks;
1337   }
1338
1339
1340 //____________________________________________________________________
1341 Bool_t AliAnalyseUE::TrackMCSelected(Double_t charge, Double_t pT, Double_t eta, Int_t pdgCode)const{
1342
1343   // MC track selection
1344   Double_t const kMyTolerance = 0.0000001;
1345   //if (charge == 0. || charge == -99.) return kFALSE;
1346   if (charge < kMyTolerance || charge + 99 < kMyTolerance) return kFALSE;
1347
1348   if ( fUseSingleCharge ) { // Charge selection
1349         if ( fUsePositiveCharge && charge < 0.) return kFALSE; // keep Positives
1350         if ( !fUsePositiveCharge && charge > 0.) return kFALSE; // keep Negatives
1351         }
1352
1353   //Kinematics cuts on particle
1354   if ((pT < fTrackPtCut) || (TMath::Abs(eta) > fTrackEtaCut )) return kFALSE;
1355
1356   Bool_t isHadron = TMath::Abs(pdgCode)==211 ||
1357         TMath::Abs(pdgCode)==2212 ||
1358         TMath::Abs(pdgCode)==321;
1359
1360   if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1361
1362   return kTRUE;
1363
1364 }
1365
1366 //____________________________________________________________________
1367 Bool_t AliAnalyseUE::TrackSelected(AliAODTrack* part)const{
1368
1369   // Real track selection
1370   if ( !part->TestFilterBit(fFilterBit) ) return kFALSE; // track cut selection
1371   if ( !part->IsPrimaryCandidate()) return kFALSE; // reject whatever is not linked to collision point
1372   // PID Selection: Reject everything but hadrons
1373   Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
1374         part->GetMostProbablePID()==AliAODTrack::kKaon ||
1375         part->GetMostProbablePID()==AliAODTrack::kProton;
1376   if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1377
1378   if ( !part->Charge() ) return kFALSE; //Only charged
1379   if ( fUseSingleCharge ) { // Charge selection
1380         if ( fUsePositiveCharge && part->Charge() < 0.) return kFALSE; // keep Positives
1381         if ( !fUsePositiveCharge && part->Charge() > 0.) return kFALSE; // keep Negatives
1382         }
1383
1384   if ( part->Pt() < fTrackPtCut ) return kFALSE;
1385   if( TMath::Abs(part->Eta()) > fTrackEtaCut ) return kFALSE;
1386
1387   return kTRUE;
1388 }
1389
1390 //____________________________________________________________________
1391 Bool_t AliAnalyseUE::TrackSelectedEfficiency(AliAODTrack* part)const{
1392
1393   if (!fStack) AliWarning("Attention: stack not set from mother task");
1394   // Real track selection
1395   if ( !part->TestFilterBit(fFilterBit) ) return kFALSE; // track cut selection
1396   if ( !part->IsPrimaryCandidate()) return kFALSE; // reject whatever is not linked to collision point
1397   // PID Selection: Reject everything but hadrons
1398   Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
1399         part->GetMostProbablePID()==AliAODTrack::kKaon ||
1400         part->GetMostProbablePID()==AliAODTrack::kProton;
1401   if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1402
1403   if ( !part->Charge() ) return kFALSE; //Only charged
1404   if ( fUseSingleCharge ) { // Charge selection
1405         if ( fUsePositiveCharge && part->Charge() < 0.) return kFALSE; // keep Positives
1406         if ( !fUsePositiveCharge && part->Charge() > 0.) return kFALSE; // keep Negatives
1407         }
1408
1409   /*if ( part->Pt() < fTrackPtCut ) return kFALSE;
1410   if( TMath::Abs(part->Eta()) > fTrackEtaCut ) return kFALSE;*/
1411
1412   //Check if there was a primary fulfilling the required cuts
1413   Double_t charge=-999.;
1414   Double_t pt=-999.;
1415   Double_t eta=-999.;
1416   Int_t pdgcode=-999;
1417
1418   Int_t label = part->GetLabel();
1419   TClonesArray *tca=0;
1420   tca = dynamic_cast<TClonesArray*>(fkAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1421   if(!tca)return kFALSE;
1422   //TParticle *partMC = fStack->ParticleFromTreeK(label);
1423   AliAODMCParticle *partMC=dynamic_cast<AliAODMCParticle*>(tca->At(TMath::Abs(label)));
1424   if(!partMC)return kFALSE;
1425   if (!partMC->IsPhysicalPrimary())return kFALSE;
1426   //charge = ((TParticlePDG*)partMC->GetPDG())->Charge();
1427   charge = partMC->Charge();
1428   pt = partMC->Pt();
1429   eta = partMC->Eta();
1430   pdgcode = partMC->GetPdgCode();
1431
1432   if (!TrackMCSelected(charge, pt, eta, pdgcode)) return kFALSE;
1433   return kTRUE;
1434 }
1435
1436