PWG4/JetTasks/AliAnalysisTaskUE.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id:$ */
  17
  18 #include <TROOT.h>
  19 #include <TSystem.h>
  20 #include <TChain.h>
  21 #include <TFile.h>
  22 #include <TList.h>
  23 #include <TH1I.h>
  24 #include <TH1F.h>
  25 #include <TH2F.h>
  26 #include <TCanvas.h>
  27 #include <TVector3.h>
  28 #include <TLorentzVector.h>
  29 #include <TMath.h>
  30 #include <TTree.h>
  31 #include <TBranch.h>
  32
  33 #include "AliAnalysisTaskUE.h"
  34 #include "AliAnalysisManager.h"
  35 #include "AliMCEventHandler.h"
  36 #include "AliAODEvent.h"
  37 #include "AliAODInputHandler.h"
  38 #include "AliAODHandler.h"
  39 #include "AliStack.h"
  40 #include "AliAODJet.h"
  41 #include "AliAODTrack.h"
  42
  43 #include "AliLog.h"
  44
  45 //
  46 // Analysis class for Underlying Event studies
  47 //
  48 // Look for correlations on the tranverse regions to
  49 // the leading charged jet
  50 //
  51 // This class needs as input AOD with track and Jets
  52 // the output is a list of histograms
  53 //
  54 // AOD can be either connected to the InputEventHandler
  55 // for a chain of AOD files
  56 // or
  57 // to the OutputEventHandler
  58 // for a chain of ESD files, so this case class should be
  59 // in the train after the Jet finder
  60 //
  61 //    Arian.Abrahantes.Quintana@cern.ch
  62 //    Ernesto.Lopez.Torres@cern.ch
  63 //
  64
  65
  66 ClassImp( AliAnalysisTaskUE)
  67
  68 ////////////////////////////////////////////////////////////////////////
  69
  70
  71 //____________________________________________________________________
  72 AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name):
  73 AliAnalysisTask(name, ""),
  74 fDebug(kFALSE),
  75 fAOD(0x0),
  76 fAODjets(0x0),
  77 fListOfHistos(0x0),
  78 fBinsPtInHist(30),
  79 fMinJetPtInHist(0.),
  80 fMaxJetPtInHist(300.),
  81 fAnaType(1),
  82 fRegionType(1),
  83 fConeRadius(0.7),
  84 fAreaReg(1.5393), // Pi*0.7*0.7
  85 fUseChPartJet(kFALSE),
  86 fUseSingleCharge(kFALSE),
  87 fUsePositiveCharge(kTRUE),
  88 fOrdering(1),
  89 fFilterBit(0xFF),
  90 fJetsOnFly(kFALSE),
  91 fJet1EtaCut(0.2),
  92 fJet2DeltaPhiCut(2.616),    // 150 degrees
  93 fJet2RatioPtCut(0.8),
  94 fJet3PtCut(15.),
  95 fTrackPtCut(0.),
  96 fTrackEtaCut(0.9),
  97 fhNJets(0x0),
  98 fhEleadingPt(0x0),
  99 fhMinRegPtDist(0x0),
 100 fhRegionMultMin(0x0),
 101 fhMinRegAvePt(0x0),
 102 fhMinRegSumPt(0x0),
 103 fhMinRegMaxPtPart(0x0),
 104 fhMinRegSumPtvsMult(0x0),
 105 fhdNdEta_PhiDist(0x0),
 106 fhFullRegPartPtDistVsEt(0x0),
 107 fhTransRegPartPtDistVsEt(0x0),
 108 fhRegionSumPtMaxVsEt(0x0),
 109 fhRegionMultMax(0x0),
 110 fhRegionMultMaxVsEt(0x0),
 111 fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0),
 112 fhRegionMultMinVsEt(0x0),
 113 fhRegionAveSumPtVsEt(0x0),
 114 fhRegionDiffSumPtVsEt(0x0),
 115 fhRegionAvePartPtMaxVsEt(0x0),
 116 fhRegionAvePartPtMinVsEt(0x0),
 117 fhRegionMaxPartPtMaxVsEt(0x0),
 118 fSettingsTree(0x0)//,   fhValidRegion(0x0)
 119 {
 120   // Default constructor
 121   // Define input and output slots here
 122   // Input slot #0 works with a TChain
 123   DefineInput(0, TChain::Class());
 124   // Output slot #0 writes into a TList container
 125   DefineOutput(0, TList::Class());
 126 }
 127
 128 //____________________________________________________________________
 129 void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/)
 130 {
 131   // Connect the input data
 132
 133   // We need AOD with tracks and jets.
 134   // Since AOD can be either connected to the InputEventHandler (input chain fron AOD files)
 135   // or to the OutputEventHandler ( AOD is created by a previus task in the train )
 136   // we need to check where it is and get the pointer to AODEvent in the right way
 137
 138   if (fDebug > 1) AliInfo("ConnectInputData() \n");
 139
 140   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 141
 142   if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
 143     fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
 144     if (fDebug > 1) AliInfo("  ==== Tracks from AliAODInputHandler");
 145     // Case when jets are reconstructed on the fly from AOD tracks
 146     // (the Jet Finder is using the AliJetAODReader) of InputEventHandler
 147     // and put in the OutputEventHandler AOD. Useful whe you want to reconstruct jets with
 148     // different parameters to default ones stored in the AOD or to use a different algorithm
 149     if( fJetsOnFly ) {
 150       handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 151       if( handler && handler->InheritsFrom("AliAODHandler") ) {
 152         fAODjets = ((AliAODHandler*)handler)->GetAOD();
 153         if (fDebug > 1) AliInfo("  ==== Jets from AliAODHandler");
 154       }
 155     } else {
 156       fAODjets = fAOD;
 157       if (fDebug > 1) AliInfo("  ==== Jets from AliAODInputHandler");
 158     }
 159   } else {
 160     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 161     if( handler && handler->InheritsFrom("AliAODHandler") ) {
 162       fAOD  = ((AliAODHandler*)handler)->GetAOD();
 163       fAODjets = fAOD;
 164       if (fDebug > 1) AliInfo("  ==== Tracks and Jets from AliAODHandler");
 165     } else {
 166       AliFatal("I can't get any AOD Event Handler");
 167       return;
 168     }
 169   }
 170 }
 171
 172 //____________________________________________________________________
 173 void  AliAnalysisTaskUE::CreateOutputObjects()
 174 {
 175   // Create the output container
 176   //
 177   if (fDebug > 1) AliInfo("CreateOutPutData()");
 178   //
 179   //  Histograms
 180   CreateHistos();
 181   fListOfHistos->SetOwner(kTRUE);
 182   //  OpenFile(0);
 183 }
 184
 185
 186 //____________________________________________________________________
 187 void  AliAnalysisTaskUE::Exec(Option_t */*option*/)
 188 {
 189   // Execute analysis for current event
 190   //
 191   if ( fDebug > 3 ) AliInfo( " Processing event..." );
 192   AnalyseUE();
 193
 194   // Post the data
 195   PostData(0, fListOfHistos);
 196 }
 197
 198 //____________________________________________________________________
 199 void  AliAnalysisTaskUE::AnalyseUE()
 200 {
 201   //
 202   // Look for correlations on the tranverse regions to
 203   // the leading charged jet
 204   //
 205
 206   // ------------------------------------------------
 207   // Find Leading Jets 1,2,3
 208   // (could be skipped if Jets are sort by Pt...)
 209   Double_t maxPtJet1 = 0.;
 210   Int_t    index1 = -1;
 211   Double_t maxPtJet2 = 0.;   // jet 2 need for back to back inclusive
 212   Int_t    index2 = -1;
 213   Double_t maxPtJet3 = 0.;   // jet 3 need for back to back exclusive
 214   Int_t    index3 = -1;
 215   TVector3 jetVect[3];
 216   Int_t nJets = 0;
 217
 218   if( !fUseChPartJet ) {
 219
 220     // Use AOD Jets
 221     nJets = fAODjets->GetNJets();
 222     //  printf("AOD %d jets \n", nJets);
 223     for( Int_t i=0; i<nJets; ++i ) {
 224       AliAODJet* jet = fAODjets->GetJet(i);
 225       Double_t jetPt = jet->Pt();//*1.666;  // FIXME Jet Pt Correction ?????!!!
 226       if( jetPt > maxPtJet1 ) {
 227         maxPtJet3 = maxPtJet2; index3 = index2;
 228         maxPtJet2 = maxPtJet1; index2 = index1;
 229         maxPtJet1 = jetPt; index1 = i;
 230       } else if( jetPt > maxPtJet2 ) {
 231         maxPtJet3 = maxPtJet2; index3 = index2;
 232         maxPtJet2 = jetPt; index2 = i;
 233       } else if( jetPt > maxPtJet3 ) {
 234         maxPtJet3 = jetPt; index3 = i;
 235       }
 236     }
 237     if( index1 != -1 ) {
 238       AliAODJet* jet = fAODjets->GetJet(index1);
 239       jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 240     }
 241     if( index2 != -1 ) {
 242       AliAODJet* jet = fAODjets->GetJet(index2);
 243       jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 244     }
 245     if( index3 != -1 ) {
 246       AliAODJet* jet = fAODjets->GetJet(index3);
 247       jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 248     }
 249
 250   } else {
 251
 252     // Use "Charged Particle Jets"
 253     TObjArray* jets = FindChargedParticleJets();
 254     if( jets ) {
 255       nJets = jets->GetEntriesFast();
 256       if( nJets > 0 ) {
 257         index1 = 0;
 258         AliAODJet* jet = (AliAODJet*)jets->At(0);
 259         maxPtJet1 = jet->Pt();
 260         jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 261       }
 262       if( nJets > 1 ) {
 263         index2 = 1;
 264         AliAODJet* jet = (AliAODJet*)jets->At(1);
 265         maxPtJet1 = jet->Pt();
 266         jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 267       }
 268       if( nJets > 2 ) {
 269         index3 = 2;
 270         AliAODJet* jet = (AliAODJet*)jets->At(2);
 271         maxPtJet1 = jet->Pt();
 272         jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
 273       }
 274
 275       jets->Delete();
 276       delete jets;
 277     }
 278   }
 279
 280
 281   fhNJets->Fill(nJets);
 282
 283   if( fDebug > 1 ) {
 284     if( index1 < 0 ) {
 285       AliInfo("\n   Skipping Event, not jet found...");
 286       return;
 287     } else
 288       AliInfo(Form("\n   Pt Leading Jet = %6.1f eta=%5.3f ",  maxPtJet1, jetVect[0].Eta() ));
 289   }
 290
 291
 292   // ----------------------------------------------
 293   // Cut events by jets topology
 294   // fAnaType:
 295   //     1 = inclusive,
 296   //         - Jet1 |eta| < fJet1EtaCut
 297   //     2 = back to back inclusive
 298   //         - fulfill case 1
 299   //         - |Jet1.Phi - Jet2.Phi| > fJet2DeltaPhiCut
 300   //         - Jet2.Pt/Jet1Pt > fJet2RatioPtCut
 301   //     3 = back to back exclusive
 302   //         - fulfill case 2
 303   //         - Jet3.Pt < fJet3PtCut
 304
 305   if( index1 < 0 || TMath::Abs(jetVect[0].Eta()) > fJet1EtaCut) {
 306     if( fDebug > 1 ) AliInfo("\n   Skipping Event...Jet1 |eta| > fJet1EtaCut");
 307     return;
 308   }
 309   // back to back inclusive
 310   if( fAnaType > 1 && index2 == -1 ) {
 311     if( fDebug > 1 ) AliInfo("\n   Skipping Event... no second Jet found");
 312     return;
 313   }
 314   if( fAnaType > 1 && index2 > -1 ) {
 315     if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
 316        maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) {
 317       if( fDebug > 1 ) AliInfo("\n   Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
 318       return;
 319     }
 320   }
 321   // back to back exclusive
 322   if( fAnaType > 2 && index3 > -1 ) {
 323     if( maxPtJet3 > fJet3PtCut ) {
 324       if( fDebug > 1 ) AliInfo("\n   Skipping Event... Jet3.Pt > fJet3PtCut ");
 325       return;
 326     }
 327   }
 328
 329   fhEleadingPt->Fill( maxPtJet1 );
 330   //Area for Normalization Purpose at Display histos
 331   SetRegionArea(jetVect);
 332
 333   // ----------------------------------------------
 334   // Find max and min regions
 335   Double_t sumPtRegionPosit = 0.;
 336   Double_t sumPtRegionNegat = 0.;
 337   Double_t maxPartPtRegion  = 0.;
 338   Int_t    nTrackRegionPosit = 0;
 339   Int_t    nTrackRegionNegat = 0;
 340   static const Double_t k270rad = 270.*TMath::Pi()/180.;
 341
 342   Int_t nTracks = fAOD->GetNTracks();
 343   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 344     AliAODTrack* part = fAOD->GetTrack( ipart );
 345     if ( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
 346     // PID Selection: Reject everything but hadrons
 347     Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion ||
 348                       part->GetMostProbablePID()==AliAODTrack::kKaon ||
 349                       part->GetMostProbablePID()==AliAODTrack::kProton;
 350     if ( !isHadron ) continue;
 351
 352     if ( !part->Charge() ) continue; //Only charged
 353     if ( fUseSingleCharge ) { // Charge selection
 354       if ( fUsePositiveCharge && part->Charge() < 0.) continue; // keep Positives
 355       if ( !fUsePositiveCharge && part->Charge() > 0.) continue; // keep Negatives
 356     }
 357
 358     if ( part->Pt() < fTrackPtCut ) continue;
 359
 360     TVector3 partVect(part->Px(), part->Py(), part->Pz());
 361
 362     Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
 363     if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
 364     fhdNdEta_PhiDist->Fill( deltaPhi );
 365     fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
 366
 367     Int_t region = IsTrackInsideRegion( jetVect, &partVect );
 368
 369     if (region > 0) {
 370       if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
 371       sumPtRegionPosit += part->Pt();
 372       nTrackRegionPosit++;
 373       fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
 374     }
 375     if (region < 0) {
 376       if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
 377       sumPtRegionNegat += part->Pt();
 378       nTrackRegionNegat++;
 379       fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
 380     }
 381   }
 382
 383   //How quantities will be sorted before Fill Min and Max Histogram
 384   //  1=Plots will be CDF-like
 385   //  2=Plots will be Marchesini-like
 386   if( fOrdering == 1 ) {
 387     if( sumPtRegionPosit > sumPtRegionNegat ) {
 388       FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
 389     } else {
 390       FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
 391     }
 392     if (nTrackRegionPosit > nTrackRegionNegat ) {
 393       FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
 394     } else {
 395       FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
 396     }
 397   } else if( fOrdering == 2 ) {
 398     if (sumPtRegionPosit > sumPtRegionNegat) {
 399       FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
 400       FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
 401     } else {
 402       FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
 403       FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
 404     }
 405   }
 406
 407   Double_t avePosRegion = (nTrackRegionPosit) ? sumPtRegionPosit/nTrackRegionPosit : 0.;
 408   Double_t aveNegRegion = (nTrackRegionNegat) ? sumPtRegionNegat/nTrackRegionNegat : 0.;
 409   if( avePosRegion > aveNegRegion ) {
 410     FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg );
 411   } else {
 412     FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg );
 413   }
 414
 415   fhRegionMaxPartPtMaxVsEt->Fill(maxPtJet1, maxPartPtRegion );
 416
 417   // Compute pedestal like magnitudes
 418   fhRegionDiffSumPtVsEt->Fill(maxPtJet1, TMath::Abs(sumPtRegionPosit-sumPtRegionNegat)/(2.0*fAreaReg));
 419   fhRegionAveSumPtVsEt->Fill(maxPtJet1, (sumPtRegionPosit+sumPtRegionNegat)/(2.0*fAreaReg));
 420
 421 }
 422
 423 //____________________________________________________________________
 424 void AliAnalysisTaskUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
 425 {
 426   // Fill sumPt of control regions
 427
 428   // Max cone
 429   fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax );
 430   // Min cone
 431   fhRegionSumPtMinVsEt->Fill( leadingE, ptMin );
 432   // MAke distributions for UE comparison with MB data
 433   fhMinRegSumPt->Fill(ptMin);
 434 }
 435
 436 //____________________________________________________________________
 437 void AliAnalysisTaskUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
 438 {
 439   // Fill average particle Pt of control regions
 440
 441   // Max cone
 442   fhRegionAvePartPtMaxVsEt->Fill( leadingE, ptMax );
 443   // Min cone
 444   fhRegionAvePartPtMinVsEt->Fill( leadingE, ptMin );
 445   // MAke distributions for UE comparison with MB data
 446   fhMinRegAvePt->Fill(ptMin);
 447 }
 448
 449 //____________________________________________________________________
 450 void AliAnalysisTaskUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin  )
 451 {
 452   // Fill Nch multiplicity of control regions
 453
 454   // Max cone
 455   fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax );
 456   fhRegionMultMax->Fill( nTrackPtmax );
 457   // Min cone
 458   fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin );
 459   fhRegionMultMin->Fill( nTrackPtmin );
 460   // MAke distributions for UE comparison with MB data
 461   fhMinRegSumPtvsMult->Fill(nTrackPtmin,ptMin);
 462 }
 463
 464 //____________________________________________________________________
 465 Int_t AliAnalysisTaskUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect)
 466 {
 467   // return de region in delta phi
 468   // -1 negative delta phi
 469   //  1 positive delta phi
 470   //  0 outside region
 471   static const Double_t k60rad  = 60.*TMath::Pi()/180.;
 472   static const Double_t k120rad = 120.*TMath::Pi()/180.;
 473
 474   Int_t region = 0;
 475   if( fRegionType == 1 ) {
 476     if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
 477     // transverse regions
 478     if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1;
 479     if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1;
 480
 481   } else if( fRegionType == 2 ) {
 482     // Cone regions
 483     Double_t deltaR = 0.;
 484
 485     TVector3 positVect,negatVect;
 486     positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
 487     negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
 488
 489     if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) {
 490       region = 1;
 491       deltaR = positVect.DrEtaPhi(*partVect);
 492     } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) {
 493       region = -1;
 494       deltaR = negatVect.DrEtaPhi(*partVect);
 495     }
 496
 497     if (deltaR > fConeRadius) region = 0;
 498
 499   } else
 500     AliError("Unknow region type");
 501
 502   // For debug (to be removed)
 503   //if( region != 0 )  fhValidRegion->Fill( partVect->Eta()-jetVect[0].Eta(), jetVect[0].DeltaPhi(*partVect) );
 504
 505   return region;
 506 }
 507
 508
 509 //____________________________________________________________________
 510 TObjArray*  AliAnalysisTaskUE::FindChargedParticleJets()
 511 {
 512   // Return a TObjArray of "charged particle jets"
 513   //
 514   // Charged particle jet deﬁnition from reference:
 515   //
 516   // "Charged jet evolution and the underlying event
 517   //  in proton-antiproton collisions at 1.8 TeV"
 518   //  PHYSICAL REVIEW D 65 092002, CDF Collaboration
 519   //
 520   // We deﬁne "jets" as circular regions in eta-phi space with
 521   // radius deﬁned by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
 522   // Our jet algorithm is as follows:
 523   //   1- Order all charged particles according to their pT .
 524   //   2- Start with the highest pT particle and include in the jet all
 525   //      particles within the radius R=0.7 considering each particle
 526   //      in the order of decreasing pT and recalculating the centroid
 527   //      of the jet after each new particle is added to the jet .
 528   //   3- Go to the next highest pT particle not already included in
 529   //      a jet and add to the jet all particles not already included in
 530   //      a jet within R=0.7.
 531   //   4- Continue until all particles are in a jet.
 532   // We deﬁne the transverse momentum of the jet to be
 533   // the scalar pT sum of all the particles within the jet, where pT
 534   // is measured with respect to the beam axis
 535
 536   //  1 - Order all charged particles according to their pT .
 537   Int_t nTracks = fAOD->GetNTracks();
 538   if( !nTracks ) return 0;
 539   TObjArray tracks(nTracks);
 540
 541   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 542     AliAODTrack* part = fAOD->GetTrack( ipart );
 543     if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
 544     if( !part->Charge() ) continue;
 545     if( part->Pt() < fTrackPtCut ) continue;
 546     tracks.AddLast(part);
 547   }
 548   QSortTracks( tracks, 0, tracks.GetEntriesFast() );
 549
 550   nTracks = tracks.GetEntriesFast();
 551   if( !nTracks ) return 0;
 552   TObjArray *jets = new TObjArray(nTracks);
 553   TIter itrack(&tracks);
 554   while( nTracks ) {
 555     // 2- Start with the highest pT particle ...
 556     Float_t px,py,pz,pt;
 557     AliAODTrack* track = (AliAODTrack*)itrack.Next();
 558     if( !track ) continue;
 559     px = track->Px();
 560     py = track->Py();
 561     pz = track->Pz();
 562     pt = track->Pt(); // Use the energy member to store Pt
 563     jets->AddLast( new TLorentzVector(px, py, pz, pt) );
 564     tracks.Remove( track );
 565     TLorentzVector* jet = (TLorentzVector*)jets->Last();
 566     // 3- Go to the next highest pT particle not already included...
 567     AliAODTrack* track1;
 568     while ( (track1  = (AliAODTrack*)(itrack.Next())) ) {
 569       Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-track1->Phi());
 570       Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +
 571                                dphi*dphi );
 572       if( r < fConeRadius ) {
 573         Double_t Pt   = jet->E()+track1->Pt();  // Scalar sum of Pt
 574         // recalculating the centroid
 575         Double_t eta = jet->Eta()*jet->E()/Pt + track1->Eta()/Pt;
 576         Double_t phi = jet->Phi()*jet->E()/Pt + track1->Phi()/Pt;
 577         jet->SetPtEtaPhiE( 1., eta, phi, Pt );
 578         tracks.Remove( track1 );
 579       }
 580     }
 581
 582     tracks.Compress();
 583     nTracks = tracks.GetEntries();
 584     //   4- Continue until all particles are in a jet.
 585     itrack.Reset();
 586   } // end while nTracks
 587
 588   // Convert to AODjets....
 589   Int_t njets = jets->GetEntriesFast();
 590   TObjArray* aodjets = new TObjArray(njets);
 591   aodjets->SetOwner(kTRUE);
 592   for(Int_t ijet=0; ijet<njets; ++ijet) {
 593     TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
 594     Float_t px, py,pz,en; // convert to 4-vector
 595     px = jet->E() * TMath::Cos(jet->Phi());  // Pt * cos(phi)
 596     py = jet->E() * TMath::Sin(jet->Phi());  // Pt * sin(phi)
 597     pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
 598     en = TMath::Sqrt(px * px + py * py + pz * pz);
 599     aodjets->AddLast( new AliAODJet(px, py, pz, en) );
 600   }
 601   // Order jets according to their pT .
 602   QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
 603
 604   // debug
 605   if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
 606
 607   return aodjets;
 608 }
 609
 610 //____________________________________________________________________
 611 void  AliAnalysisTaskUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
 612 {
 613   // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
 614
 615   static TObject *tmp;
 616   static int i;           // "static" to save stack space
 617   int j;
 618
 619   while (last - first > 1) {
 620     i = first;
 621     j = last;
 622     for (;;) {
 623       while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
 624         ;
 625       while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
 626         ;
 627       if (i >= j)
 628         break;
 629
 630       tmp  = a[i];
 631       a[i] = a[j];
 632       a[j] = tmp;
 633     }
 634     if (j == first) {
 635       ++first;
 636       continue;
 637     }
 638     tmp = a[first];
 639     a[first] = a[j];
 640     a[j] = tmp;
 641     if (j - first < last - (j + 1)) {
 642       QSortTracks(a, first, j);
 643       first = j + 1;   // QSortTracks(j + 1, last);
 644     } else {
 645       QSortTracks(a, j + 1, last);
 646       last = j;        // QSortTracks(first, j);
 647     }
 648   }
 649 }
 650
 651 //____________________________________________________________________
 652 void AliAnalysisTaskUE::SetRegionArea(TVector3 *jetVect)
 653 {
 654   Double_t fAreaCorrFactor=0.;
 655   Double_t deltaEta = 0.;
 656   if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.;
 657   else if (fRegionType==2){
 658     deltaEta = 0.9-TMath::Abs(jetVect[0].Eta());
 659     if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius;
 660     else{
 661       fAreaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
 662       (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
 663       fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-fAreaCorrFactor;
 664     }
 665   }else AliWarning("Unknown Rgion Type");
 666   if (fDebug>5) AliInfo(Form("\n dEta=%5.3f Angle =%5.3f Region Area = %5.3f Corr Factor=%5.4f \n",deltaEta,TMath::ACos(deltaEta/fConeRadius),fAreaReg,fAreaCorrFactor));
 667 }
 668
 669 //____________________________________________________________________
 670 void  AliAnalysisTaskUE::CreateHistos()
 671 {
 672   fListOfHistos = new TList();
 673
 674
 675   fhNJets = new TH1F("fhNJets", "n Jet",  10, 0, 10);
 676   fhNJets->SetXTitle("# of jets");
 677   fhNJets->Sumw2();
 678   fListOfHistos->Add( fhNJets );                 // At(0)
 679
 680   fhEleadingPt  = new TH1F("hEleadingPt",   "Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 681   fhEleadingPt->SetXTitle("P_{T} (GeV/c)");
 682   fhEleadingPt->SetYTitle("dN/dP_{T}");
 683   fhEleadingPt->Sumw2();
 684   fListOfHistos->Add( fhEleadingPt );            // At(1)
 685
 686   fhMinRegPtDist = new TH1F("hMinRegPtDist",   "P_{T} distribution in Min zone",  50,0.,20.);
 687   fhMinRegPtDist->SetXTitle("P_{T} (GeV/c)");
 688   fhMinRegPtDist->SetYTitle("dN/dP_{T}");
 689   fhMinRegPtDist->Sumw2();
 690   fListOfHistos->Add( fhMinRegPtDist );          // At(2)
 691
 692   fhRegionMultMin = new TH1F("hRegionMultMin",      "N_{ch}^{90, min}",  21, -0.5,   20.5);
 693   fhRegionMultMin->SetXTitle("N_{ch tracks}");
 694   fhRegionMultMin->Sumw2();
 695   fListOfHistos->Add( fhRegionMultMin );         // At(3)
 696
 697   fhMinRegAvePt = new TH1F("hMinRegAvePt", "#LTp_{T}#GT",  50, 0.,   20.);
 698   fhMinRegAvePt->SetXTitle("P_{T} (GeV/c)");
 699   fhMinRegAvePt->Sumw2();
 700   fListOfHistos->Add( fhMinRegAvePt );           // At(4)
 701
 702   fhMinRegSumPt = new TH1F("hMinRegSumPt", "#Sigma p_{T} ",  50, 0.,   20.);
 703   fhMinRegSumPt->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");
 704   fhMinRegSumPt->SetXTitle("#Sigma p_{T} (GeV/c)");
 705   fhMinRegSumPt->Sumw2();
 706   fListOfHistos->Add( fhMinRegSumPt );           // At(5)
 707
 708   fhMinRegMaxPtPart = new TH1F("hMinRegMaxPtPart", "max(p_{T})|_{event} ",  50, 0.,   20.);
 709   fhMinRegMaxPtPart->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");
 710   fhMinRegMaxPtPart->SetXTitle("p_{T} (GeV/c)");
 711   fhMinRegMaxPtPart->Sumw2();
 712   fListOfHistos->Add( fhMinRegMaxPtPart );       // At(6)
 713
 714   fhMinRegSumPtvsMult = new TH1F("hMinRegSumPtvsMult", "#Sigma p_{T} vs. Multiplicity ",  21, -0.5,   20.5);
 715   fhMinRegSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");
 716   fhMinRegSumPtvsMult->SetXTitle("N_{charge}");
 717   fhMinRegSumPtvsMult->Sumw2();
 718   fListOfHistos->Add( fhMinRegSumPtvsMult );     // At(7);
 719
 720   fhdNdEta_PhiDist  = new TH1F("hdNdEta_PhiDist",   "Charge particle density |#eta|< 1 vs #Delta#phi",  120, 0.,   2.*TMath::Pi());
 721   fhdNdEta_PhiDist->SetXTitle("#Delta#phi");
 722   fhdNdEta_PhiDist->SetYTitle("dN_{ch}/d#etad#phi");
 723   fhdNdEta_PhiDist->Sumw2();
 724   fListOfHistos->Add( fhdNdEta_PhiDist );        // At(8)
 725
 726   // Can be use to get part pt distribution for differente Jet Pt bins
 727   fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", "dN/dP_{T} |#eta|<1 vs Leading Jet P_{T}",
 728                                      50,0.,50., fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 729   fhFullRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
 730   fhFullRegPartPtDistVsEt->SetXTitle("p_{T}");
 731   fhFullRegPartPtDistVsEt->Sumw2();
 732   fListOfHistos->Add( fhFullRegPartPtDistVsEt );  // At(9)
 733
 734   // Can be use to get part pt distribution for differente Jet Pt bins
 735   fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", "dN/dP_{T} in tranvese regions |#eta|<1 vs Leading Jet P_{T}",
 736                                       50,0.,50., fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 737   fhTransRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
 738   fhTransRegPartPtDistVsEt->SetXTitle("p_{T}");
 739   fhTransRegPartPtDistVsEt->Sumw2();
 740   fListOfHistos->Add( fhTransRegPartPtDistVsEt );  // At(10)
 741
 742
 743   fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt",  "P_{T}^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 744   fhRegionSumPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
 745   fhRegionSumPtMaxVsEt->Sumw2();
 746   fListOfHistos->Add( fhRegionSumPtMaxVsEt );      // At(11)
 747
 748   fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt",   "P_{T}^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 749   fhRegionSumPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
 750   fhRegionSumPtMinVsEt->Sumw2();
 751   fListOfHistos->Add( fhRegionSumPtMinVsEt );      // At(12)
 752
 753   fhRegionMultMax = new TH1I("hRegionMultMax",      "N_{ch}^{90, max}",  21, -0.5,   20.5);
 754   fhRegionMultMax->SetXTitle("N_{ch tracks}");
 755   fhRegionMultMax->Sumw2();
 756   fListOfHistos->Add( fhRegionMultMax );           // At(13)
 757
 758   fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt",  "N_{ch}^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 759   fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)");
 760   fhRegionMultMaxVsEt->Sumw2();
 761   fListOfHistos->Add( fhRegionMultMaxVsEt );       // At(14)
 762
 763   fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt",  "N_{ch}^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 764   fhRegionMultMinVsEt->SetXTitle("E (GeV/c)");
 765   fhRegionMultMinVsEt->Sumw2();
 766   fListOfHistos->Add( fhRegionMultMinVsEt );      // At(15)
 767
 768   fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 769   fhRegionAveSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
 770   fhRegionAveSumPtVsEt->Sumw2();
 771   fListOfHistos->Add( fhRegionAveSumPtVsEt );     // At(16)
 772
 773   fhRegionDiffSumPtVsEt= new TH1F("hRegionPtDiffVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 774   fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
 775   fhRegionDiffSumPtVsEt->Sumw2();
 776   fListOfHistos->Add( fhRegionDiffSumPtVsEt );    // At(17)
 777
 778   fhRegionAvePartPtMaxVsEt = new TH1F("hRegionAvePartPtMaxVsEt", "#LTp_{T}#GT^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 779   fhRegionAvePartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
 780   fhRegionAvePartPtMaxVsEt->Sumw2();
 781   fListOfHistos->Add( fhRegionAvePartPtMaxVsEt );  // At(18)
 782
 783   fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 784   fhRegionAvePartPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
 785   fhRegionAvePartPtMinVsEt->Sumw2();
 786   fListOfHistos->Add( fhRegionAvePartPtMinVsEt );   // At(19)
 787
 788   fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 789   fhRegionMaxPartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
 790   fhRegionMaxPartPtMaxVsEt->Sumw2();
 791   fListOfHistos->Add( fhRegionMaxPartPtMaxVsEt );    // At(20)
 792
 793   fSettingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
 794   fListOfHistos->Add( fSettingsTree );    // At(21)
 795
 796   /*
 797    // For debug region selection
 798    fhValidRegion = new TH2F("hValidRegion", "dN/d#eta/d#phi",
 799    20, -2.,2., 62, -TMath::Pi(),   TMath::Pi());
 800    fhValidRegion->SetYTitle("#Delta#phi");
 801    fhValidRegion->Sumw2();
 802    fListOfHistos->Add( fhValidRegion );  // At(15)
 803    */
 804 }
 805
 806 //____________________________________________________________________
 807 void  AliAnalysisTaskUE::Terminate(Option_t */*option*/)
 808 {
 809   // Terminate analysis
 810   //
 811
 812   //Save Analysis Settings
 813   WriteSettings();
 814
 815   // Normalize histos to region area TODO:
 816   // Normalization done at Analysis, taking into account
 817   // area variations on per-event basis (cone case)
 818
 819   // Draw some Test plot to the screen
 820   if (!gROOT->IsBatch()) {
 821
 822     // Update pointers reading them from the output slot
 823     fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
 824     if( !fListOfHistos ) {
 825       AliError("Histogram List is not available");
 826       return;
 827     }
 828     fhEleadingPt         = (TH1F*)fListOfHistos->At(1);
 829     fhdNdEta_PhiDist     = (TH1F*)fListOfHistos->At(8);
 830     fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11);
 831     fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12);
 832     fhRegionMultMaxVsEt  = (TH1F*)fListOfHistos->At(14);
 833     fhRegionMultMinVsEt  = (TH1F*)fListOfHistos->At(15);
 834     fhRegionAveSumPtVsEt = (TH1F*)fListOfHistos->At(16);
 835
 836     fhNJets              = (TH1F*)fListOfHistos->At(0);
 837
 838     // Canvas
 839     TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1200,800);
 840     c1->Divide(2,2);
 841     c1->cd(1);
 842     TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
 843     h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1);
 844     //h1r->Scale( areafactor );
 845     h1r->SetMarkerStyle(20);
 846     h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
 847     h1r->SetYTitle("P_{T}^{90, max}");
 848     h1r->DrawCopy("p");
 849
 850     c1->cd(2);
 851
 852     TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
 853     h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1);
 854     //h2r->Scale( areafactor );
 855     h2r->SetMarkerStyle(20);
 856     h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
 857     h2r->SetYTitle("P_{T}^{90, min}");
 858     h2r->DrawCopy("p");
 859
 860     c1->cd(3);
 861     TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
 862     h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1);
 863     h4r->Scale(2.); // make average
 864     //h4r->Scale( areafactor );
 865     h4r->SetYTitle("#DeltaP_{T}^{90}");
 866     h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
 867     h4r->SetMarkerStyle(20);
 868     h4r->DrawCopy("p");
 869
 870     c1->cd(4);
 871     TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading",   "",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 872     TH1F *h6r = new TH1F("hRegionMultMinVsEtleading",   "",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 873
 874     h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1);
 875     h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1);
 876     //h5r->Scale( areafactor );
 877     //h6r->Scale( areafactor );
 878     h5r->SetYTitle("N_{Tracks}^{90}");
 879     h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
 880     h5r->SetMarkerStyle(20);
 881     h5r->DrawCopy("p");
 882     h6r->SetMarkerStyle(21);
 883     h6r->SetMarkerColor(2);
 884     h6r->SetYTitle("N_{Tracks}^{90}");
 885     h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
 886     h6r->DrawCopy("p same");
 887     c1->Update();
 888
 889     TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800);
 890     c2->Divide(2,2);
 891     c2->cd(1);
 892     fhEleadingPt->SetMarkerStyle(20);
 893     fhEleadingPt->SetMarkerColor(2);
 894     fhEleadingPt->DrawCopy("p");
 895     gPad->SetLogy();
 896
 897     c2->cd(2);
 898     fhdNdEta_PhiDist->SetMarkerStyle(20);
 899     fhdNdEta_PhiDist->SetMarkerColor(2);
 900     fhdNdEta_PhiDist->DrawCopy("p");
 901     gPad->SetLogy();
 902
 903     c2->cd(3);
 904     fhNJets->DrawCopy();
 905
 906     //fhTransRegPartPtDist  = (TH1F*)fListOfHistos->At(2);
 907     fhRegionMultMin          = (TH1F*)fListOfHistos->At(3);
 908     fhMinRegAvePt     = (TH1F*)fListOfHistos->At(4);
 909     fhMinRegSumPt     = (TH1F*)fListOfHistos->At(5);
 910     //fhMinRegMaxPtPart     = (TH1F*)fListOfHistos->At(6);
 911     fhMinRegSumPtvsMult   = (TH1F*)fListOfHistos->At(7);
 912
 913     // Canvas
 914     TCanvas* c3 = new TCanvas("c3"," p_{T} dist",160,160,1200,800);
 915     c3->Divide(2,2);
 916     c3->cd(1);
 917     /*fhTransRegPartPtDist->SetMarkerStyle(20);
 918      fhTransRegPartPtDist->SetMarkerColor(2);
 919      fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries());
 920      fhTransRegPartPtDist->DrawCopy("p");
 921      gPad->SetLogy();
 922      */
 923     c3->cd(2);
 924     fhMinRegSumPt->SetMarkerStyle(20);
 925     fhMinRegSumPt->SetMarkerColor(2);
 926     //fhMinRegSumPt->Scale(areafactor);
 927     fhMinRegSumPt->DrawCopy("p");
 928     gPad->SetLogy();
 929
 930     c3->cd(3);
 931     fhMinRegAvePt->SetMarkerStyle(20);
 932     fhMinRegAvePt->SetMarkerColor(2);
 933     //fhMinRegAvePt->Scale(areafactor);
 934     fhMinRegAvePt->DrawCopy("p");
 935     gPad->SetLogy();
 936
 937     c3->cd(4);
 938
 939     TH1F *h7r = new TH1F("hRegionMultMinVsMult",   "",  21, -0.5,   20.5);
 940
 941     h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1);
 942
 943     h7r->SetMarkerStyle(20);
 944     h7r->SetMarkerColor(2);
 945     h7r->DrawCopy("p");
 946
 947     c3->Update();
 948
 949
 950     /*    c2->cd(3);
 951      fhValidRegion->SetMarkerStyle(20);
 952      fhValidRegion->SetMarkerColor(2);
 953      fhValidRegion->DrawCopy("p");
 954      */
 955     c2->Update();
 956   } else {
 957     AliInfo(" Batch mode, not histograms will be shown...");
 958   }
 959
 960   if( fDebug > 1 )
 961     AliInfo("End analysis");
 962
 963 }
 964
 965 void  AliAnalysisTaskUE::WriteSettings()
 966 {
 967
 968   fSettingsTree->Branch("fAnaType", &fAnaType, "Ana/I");
 969   fSettingsTree->Branch("fRegionType", &fRegionType,"Reg/I");
 970   fSettingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D");
 971   fSettingsTree->Branch("fUseChPartJet", &fUseChPartJet,"UseChPart/O");
 972   fSettingsTree->Branch("fUseSingleCharge", &fUseSingleCharge,"UseSingleCh/O");
 973   fSettingsTree->Branch("fUsePositiveCharge", &fUsePositiveCharge,"UsePositiveCh/O");
 974   fSettingsTree->Branch("fOrdering", &fOrdering,"OrderMeth/I");
 975   fSettingsTree->Branch("fJet1EtaCut", &fJet1EtaCut, "LeadJetEtaCut/D");
 976   fSettingsTree->Branch("fJet2DeltaPhiCut", &fJet2DeltaPhiCut, "DeltaPhi/D");
 977   fSettingsTree->Branch("fJet2RatioPtCut", &fJet2RatioPtCut, "Jet2Ratio/D");
 978   fSettingsTree->Branch("fJet3PtCut", &fJet3PtCut, "Jet3PtCut/D");
 979   fSettingsTree->Branch("fTrackPtCut", &fTrackPtCut, "TrackPtCut/D");
 980   fSettingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
 981   fSettingsTree->Fill();
 982
 983   if (fDebug>5){
 984     AliInfo(" All Analysis Settings in Saved Tree");
 985     fSettingsTree->Scan();
 986   }
 987 }