PWG4/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
  31 #include "AliAnalysisTaskUE.h"
  32 #include "AliAnalysisManager.h"
  33 #include "AliMCEventHandler.h"
  34 #include "AliAODEvent.h"
  35 #include "AliAODInputHandler.h"
  36 #include "AliAODHandler.h"
  37 #include "AliStack.h"
  38 #include "AliAODJet.h"
  39 #include "AliAODTrack.h"
  40
  41 #include "AliLog.h"
  42
  43 //
  44 // Analysis class for Underlying Event studies
  45 //
  46 // Look for correlations on the tranverse regions to
  47 // the leading charged jet
  48 //
  49 // This class needs as input AOD with track and Jets
  50 // the output is a list of histograms
  51 //
  52 // AOD can be either connected to the InputEventHandler
  53 // for a chain of AOD files
  54 // or
  55 // to the OutputEventHandler
  56 // for a chain of ESD files, so this case class should be
  57 // in the train after the Jet finder
  58 //
  59 //    Arian.Abrahantes.Quintana@cern.ch
  60 //    Ernesto.Lopez.Torres@cern.ch
  61 //
  62
  63
  64 ClassImp( AliAnalysisTaskUE)
  65
  66 ////////////////////////////////////////////////////////////////////////
  67
  68
  69 //____________________________________________________________________
  70 AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name):
  71   AliAnalysisTask(name, ""),
  72   fDebug(kFALSE),
  73   fAOD(0x0),
  74   fListOfHistos(0x0),
  75   fBinsPtInHist(30),
  76   fMinJetPtInHist(0.),
  77   fMaxJetPtInHist(300.),
  78   fAnaType(1),
  79   fRegionType(1),
  80   fConeRadius(0.7),
  81   fOrdering(1),
  82   fJet1EtaCut(0.2),
  83   fJet2DeltaPhiCut(2.616),    // 150 degrees
  84   fJet2RatioPtCut(0.8),
  85   fJet3PtCut(15.),
  86   fTrackPtCut(0.),
  87   fTrackEtaCut(0.9),
  88   fhNJets(0x0),
  89   fhEleadingPt(0x0),
  90   fhMinRegPtDist(0x0),
  91   fhRegionMultMin(0x0),
  92   fhMinRegAvePt(0x0),
  93   fhMinRegSumPt(0x0),
  94   fhMinRegMaxPtPart(0x0),
  95   fhMinRegSumPtvsMult(0x0),
  96   fhdNdEta_PhiDist(0x0),
  97   fhFullRegPartPtDistVsEt(0x0),
  98   fhTransRegPartPtDistVsEt(0x0),
  99   fhRegionSumPtMaxVsEt(0x0),
 100   fhRegionMultMax(0x0),
 101   fhRegionMultMaxVsEt(0x0),
 102   fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0),
 103   fhRegionMultMinVsEt(0x0),
 104   fhRegionAveSumPtVsEt(0x0),
 105   fhRegionDiffSumPtVsEt(0x0),
 106   fhRegionAvePartPtMaxVsEt(0x0),
 107   fhRegionAvePartPtMinVsEt(0x0),
 108   fhRegionMaxPartPtMaxVsEt(0x0)//,   fhValidRegion(0x0)
 109 {
 110   // Default constructor
 111   // Define input and output slots here
 112   // Input slot #0 works with a TChain
 113   DefineInput(0, TChain::Class());
 114   // Output slot #0 writes into a TList container
 115   DefineOutput(0, TList::Class());
 116 }
 117
 118
 119 //____________________________________________________________________
 120 void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/)
 121 {
 122 // Connect the input data
 123
 124 // We need AOD with tracks and jets.
 125 // Since AOD can be either connected to the InputEventHandler or to the OutputEventHandler
 126 // we need to check where it is and get the pointer to AODEvent in the right way
 127
 128   if (fDebug > 1) AliInfo("ConnectInputData() \n");
 129
 130
 131   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 132
 133   if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
 134      fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
 135   } else {
 136       handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 137       if( handler && handler->InheritsFrom("AliAODHandler") ) {
 138          fAOD  = ((AliAODHandler*)handler)->GetAOD();
 139       } else {
 140         AliFatal("I can't get any AOD Event Handler");
 141       }
 142    }
 143 }
 144
 145 //____________________________________________________________________
 146 void  AliAnalysisTaskUE::CreateOutputObjects()
 147 {
 148 // Create the output container
 149 //
 150   if (fDebug > 1) AliInfo("CreateOutPutData()");
 151   //
 152   //  Histograms
 153   CreateHistos();
 154   fListOfHistos->SetOwner(kTRUE);
 155 //  OpenFile(0);
 156 }
 157
 158
 159 //____________________________________________________________________
 160 void  AliAnalysisTaskUE::Exec(Option_t */*option*/)
 161 {
 162 // Execute analysis for current event
 163 //
 164   AnalyseUE();
 165   // Post the data
 166   PostData(0, fListOfHistos);
 167 }
 168
 169
 170 //____________________________________________________________________
 171 void  AliAnalysisTaskUE::Terminate(Option_t */*option*/)
 172 {
 173 // Terminate analysis
 174 //
 175
 176   // Normalize histos to region area TODO:
 177   Double_t areafactor = 1.;
 178   if( fRegionType == 1 ) {
 179     areafactor = 1./ (fTrackEtaCut *2. * TMath::Pi()*2.);
 180   } else {
 181     areafactor = 1./ ( fConeRadius * fConeRadius * TMath::Pi() );
 182   }
 183
 184   // Draw some Test plot to the screen
 185   if (!gROOT->IsBatch()) {
 186
 187     // Update pointers reading them from the output slot
 188     fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
 189     if( !fListOfHistos ) {
 190       AliError("Histogram List is not available");
 191       return;
 192     }
 193     fhEleadingPt         = (TH1F*)fListOfHistos->At(1);
 194     fhdNdEta_PhiDist     = (TH1F*)fListOfHistos->At(8);
 195     fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11);
 196     fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12);
 197     fhRegionMultMaxVsEt  = (TH1F*)fListOfHistos->At(14);
 198     fhRegionMultMinVsEt  = (TH1F*)fListOfHistos->At(15);
 199     fhRegionAveSumPtVsEt = (TH1F*)fListOfHistos->At(16);
 200
 201     fhNJets              = (TH1F*)fListOfHistos->At(0);
 202
 203     // Canvas
 204     TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1200,800);
 205     c1->Divide(2,2);
 206     c1->cd(1);
 207     TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
 208     h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1);
 209     h1r->Scale( areafactor );
 210     h1r->SetMarkerStyle(20);
 211     h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
 212     h1r->SetYTitle("P_{T}^{90, max}");
 213     h1r->DrawCopy("p");
 214
 215     c1->cd(2);
 216
 217     TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
 218     h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1);
 219     h2r->Scale( areafactor );
 220     h2r->SetMarkerStyle(20);
 221     h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
 222     h2r->SetYTitle("P_{T}^{90, min}");
 223     h2r->DrawCopy("p");
 224
 225     c1->cd(3);
 226     TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
 227     h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1);
 228     h4r->Scale(2.); // make average
 229     h4r->Scale( areafactor );
 230     h4r->SetYTitle("#DeltaP_{T}^{90}");
 231     h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
 232     h4r->SetMarkerStyle(20);
 233     h4r->DrawCopy("p");
 234
 235     c1->cd(4);
 236     TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading",   "",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 237     TH1F *h6r = new TH1F("hRegionMultMinVsEtleading",   "",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 238
 239     h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1);
 240     h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1);
 241     h5r->Scale( areafactor );
 242     h6r->Scale( areafactor );
 243     h5r->SetYTitle("N_{Tracks}^{90}");
 244     h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
 245     h5r->SetMarkerStyle(20);
 246     h5r->DrawCopy("p");
 247     h6r->SetMarkerStyle(21);
 248     h6r->SetMarkerColor(2);
 249     h6r->SetYTitle("N_{Tracks}^{90}");
 250     h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
 251     h6r->DrawCopy("p same");
 252     c1->Update();
 253
 254     TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800);
 255     c2->Divide(2,2);
 256     c2->cd(1);
 257     fhEleadingPt->SetMarkerStyle(20);
 258     fhEleadingPt->SetMarkerColor(2);
 259     fhEleadingPt->DrawCopy("p");
 260     gPad->SetLogy();
 261
 262     c2->cd(2);
 263     fhdNdEta_PhiDist->SetMarkerStyle(20);
 264     fhdNdEta_PhiDist->SetMarkerColor(2);
 265     fhdNdEta_PhiDist->DrawCopy("p");
 266     gPad->SetLogy();
 267
 268     c2->cd(3);
 269     fhNJets->DrawCopy();
 270
 271     //fhTransRegPartPtDist  = (TH1F*)fListOfHistos->At(2);
 272     fhRegionMultMin          = (TH1F*)fListOfHistos->At(3);
 273     fhMinRegAvePt     = (TH1F*)fListOfHistos->At(4);
 274     fhMinRegSumPt     = (TH1F*)fListOfHistos->At(5);
 275     //fhMinRegMaxPtPart     = (TH1F*)fListOfHistos->At(6);
 276     fhMinRegSumPtvsMult   = (TH1F*)fListOfHistos->At(7);
 277
 278     // Canvas
 279     TCanvas* c3 = new TCanvas("c3"," p_{T} dist",160,160,1200,800);
 280     c3->Divide(2,2);
 281     c3->cd(1);
 282     /*fhTransRegPartPtDist->SetMarkerStyle(20);
 283     fhTransRegPartPtDist->SetMarkerColor(2);
 284     fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries());
 285     fhTransRegPartPtDist->DrawCopy("p");
 286     gPad->SetLogy();
 287     */
 288     c3->cd(2);
 289     fhMinRegSumPt->SetMarkerStyle(20);
 290     fhMinRegSumPt->SetMarkerColor(2);
 291     fhMinRegSumPt->Scale(areafactor);
 292     fhMinRegSumPt->DrawCopy("p");
 293     gPad->SetLogy();
 294
 295     c3->cd(3);
 296     fhMinRegAvePt->SetMarkerStyle(20);
 297     fhMinRegAvePt->SetMarkerColor(2);
 298     fhMinRegAvePt->Scale(areafactor);
 299     fhMinRegAvePt->DrawCopy("p");
 300     gPad->SetLogy();
 301
 302     c3->cd(4);
 303
 304     TH1F *h7r = new TH1F("hRegionMultMinVsMult",   "",  21, -0.5,   20.5);
 305
 306     h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1);
 307
 308     h7r->SetMarkerStyle(20);
 309     h7r->SetMarkerColor(2);
 310     h7r->DrawCopy("p");
 311
 312     c3->Update();
 313
 314
 315 /*    c2->cd(3);
 316     fhValidRegion->SetMarkerStyle(20);
 317     fhValidRegion->SetMarkerColor(2);
 318     fhValidRegion->DrawCopy("p");
 319 */
 320     c2->Update();
 321   } else {
 322     AliInfo(" Batch mode, not histograms will shown...");
 323   }
 324
 325   if( fDebug > 1 )
 326      AliInfo("End analysis");
 327
 328 }
 329
 330
 331 //____________________________________________________________________
 332 void  AliAnalysisTaskUE::AnalyseUE()
 333 {
 334   //
 335   // Look for correlations on the tranverse regions to
 336   // the leading charged jet
 337   //
 338   // ------------------------------------------------
 339   // Find Leading Jets 1,2,3
 340   // (could be skipped if Jets are sort by Pt...)
 341   Double_t maxPtJet1 = 0.;
 342   Int_t    index1 = -1;
 343   Double_t maxPtJet2 = 0.;   // jet 2 need for back to back inclusive
 344   Int_t    index2 = -1;
 345   Double_t maxPtJet3 = 0.;   // jet 3 need for back to back exclusive
 346   Int_t    index3 = -1;
 347   Int_t nJets = fAOD->GetNJets();
 348   for( Int_t i=0; i<nJets; ++i ) {
 349     AliAODJet* jet = fAOD->GetJet(i);
 350     if( jet->Pt() > maxPtJet1 ) {
 351       maxPtJet3 = maxPtJet2; index3 = index2;
 352       maxPtJet2 = maxPtJet1; index2 = index1;
 353       maxPtJet1 = jet->Pt(); index1 = i;
 354     } else if( jet->Pt() > maxPtJet2 ) {
 355       maxPtJet3 = maxPtJet2; index3 = index2;
 356       maxPtJet2 = jet->Pt(); index2 = i;
 357     } else if( jet->Pt() > maxPtJet3 ) {
 358       maxPtJet3 = jet->Pt(); index3 = i;
 359     }
 360   }
 361
 362         if( index1 < 0 ) {
 363       AliInfo("\nSkipping Event, not jet found...");
 364       return;
 365         } else
 366       AliInfo(Form("\nPt Leading Jet = %6.1f eta=%5.1f ",  maxPtJet1, fAOD->GetJet(index1)->Eta() ));
 367
 368   fhNJets->Fill(nJets);
 369
 370   TVector3 jetVect[2];
 371
 372   // ----------------------------------------------
 373   // Cut events by jets topology
 374   // fAnaType:
 375   //     1 = inclusive,
 376   //         - Jet1 |eta| < fJet1EtaCut
 377   //     2 = back to back inclusive
 378   //         - fulfill case 1
 379   //         - |Jet1.Phi - Jet2.Phi| > fJet2DeltaPhiCut
 380   //         - Jet2.Pt/Jet1Pt > fJet2RatioPtCut
 381   //     3 = back to back exclusive
 382   //         - fulfill case 2
 383   //         - Jet3.Pt < fJet3PtCut
 384
 385   Bool_t isInclusive = kFALSE;
 386
 387   if(  TMath::Abs(fAOD->GetJet(index1)->Eta()) > fJet1EtaCut) {
 388     if( fDebug > 1 ) AliInfo("Skipping Event...Jet1 |eta| > fJet1EtaCut");
 389     return;
 390   }
 391   isInclusive = kTRUE;
 392   jetVect[0].SetXYZ(fAOD->GetJet(index1)->Px(),
 393                     fAOD->GetJet(index1)->Py(),
 394                     fAOD->GetJet(index1)->Pz());
 395   // back to back inclusive
 396   Bool_t isB2Binclusive = kFALSE;
 397   if( fAnaType > 1 && index2 > 0 && isInclusive) {
 398     jetVect[1].SetXYZ(fAOD->GetJet(index2)->Px(),
 399                       fAOD->GetJet(index2)->Py(),
 400                       fAOD->GetJet(index2)->Pz());
 401     if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
 402         maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) {
 403       if( fDebug > 1 ) AliInfo("Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
 404       return;
 405     }
 406     isB2Binclusive = kTRUE;
 407   }
 408   if (isInclusive && !isB2Binclusive && fAnaType>1) return;
 409   // back to back exclusive
 410   Bool_t isB2Bexclusive = kFALSE;
 411   if( fAnaType > 2 && index3 > 0 && isB2Binclusive) {
 412     if( fAOD->GetJet(index3)->Pt() > fJet3PtCut ) {
 413       if( fDebug > 1 ) AliInfo("Skipping Event... Jet3.Pt > fJet3PtCut ");
 414       return;
 415     }
 416     isB2Bexclusive = kTRUE;
 417   }
 418   if (isB2Binclusive && !isB2Bexclusive && fAnaType>2) return;
 419
 420   AliInfo(Form("njets = %d",nJets));
 421   fhEleadingPt->Fill( maxPtJet1 );
 422
 423   // ----------------------------------------------
 424   // Find max and min regions
 425   Double_t sumPtRegionPosit = 0.;
 426   Double_t sumPtRegionNegat = 0.;
 427   Double_t maxPartPtRegion  = 0.;
 428   Int_t    nTrackRegionPosit = 0;
 429   Int_t    nTrackRegionNegat = 0;
 430   static const Double_t k270rad = 270.*TMath::Pi()/180.;
 431
 432   Int_t nTracks = fAOD->GetNTracks();
 433   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
 434     AliAODTrack* part = fAOD->GetTrack( ipart );
 435
 436     if ( !part->Charge() ) continue;
 437     if ( part->Pt() < fTrackPtCut ) continue;
 438
 439     TVector3 partVect(part->Px(), part->Py(), part->Pz());
 440
 441     Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
 442     if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
 443     fhdNdEta_PhiDist->Fill( deltaPhi );
 444     fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
 445
 446     Int_t region = IsTrackInsideRegion( jetVect, &partVect );
 447
 448     if (region > 0) {
 449        if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
 450        sumPtRegionPosit += part->Pt();
 451        nTrackRegionPosit++;
 452        fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
 453     }
 454     if (region < 0) {
 455        if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
 456        sumPtRegionNegat += part->Pt();
 457        nTrackRegionNegat++;
 458        fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
 459     }
 460   }
 461
 462   //How quantities will be sorted before Fill Min and Max Histogram
 463   //  1=Plots will be CDF-like
 464   //  2=Plots will be Marchesini-like
 465   if (fOrdering==1){
 466     if (sumPtRegionPosit > sumPtRegionNegat) {
 467       FillSumPtRegion( maxPtJet1, sumPtRegionPosit, sumPtRegionNegat );
 468     } else {
 469       FillSumPtRegion( maxPtJet1, sumPtRegionNegat, sumPtRegionPosit );
 470     }
 471     if (nTrackRegionPosit > nTrackRegionNegat ) {
 472       FillMultRegion( maxPtJet1, nTrackRegionPosit, nTrackRegionNegat, sumPtRegionNegat );
 473     } else {
 474       FillMultRegion( maxPtJet1, nTrackRegionNegat, nTrackRegionPosit, sumPtRegionPosit );
 475     }
 476   } else if (fOrdering==2) {
 477     if (sumPtRegionPosit > sumPtRegionNegat) {
 478       FillSumPtRegion( maxPtJet1, sumPtRegionPosit, sumPtRegionNegat );
 479       FillMultRegion( maxPtJet1, nTrackRegionPosit, nTrackRegionNegat, sumPtRegionNegat );
 480     } else {
 481       FillSumPtRegion( maxPtJet1, sumPtRegionNegat, sumPtRegionPosit );
 482       FillMultRegion( maxPtJet1, nTrackRegionNegat, nTrackRegionPosit, sumPtRegionPosit );
 483     }
 484   }
 485
 486   Double_t avePosRegion = (nTrackRegionPosit) ? sumPtRegionPosit/nTrackRegionPosit : 0.;
 487   Double_t aveNegRegion = (nTrackRegionNegat) ? sumPtRegionNegat/nTrackRegionNegat : 0.;
 488   if (avePosRegion > aveNegRegion) {
 489     FillAvePartPtRegion( maxPtJet1, avePosRegion, aveNegRegion );
 490   } else {
 491     FillAvePartPtRegion( maxPtJet1, aveNegRegion, avePosRegion );
 492   }
 493
 494   fhRegionMaxPartPtMaxVsEt->Fill(maxPtJet1, maxPartPtRegion );
 495
 496   // Compute pedestal like magnitude
 497   fhRegionDiffSumPtVsEt->Fill(maxPtJet1, TMath::Abs(sumPtRegionPosit-sumPtRegionNegat)/2.0);
 498   fhRegionAveSumPtVsEt->Fill(maxPtJet1, (sumPtRegionPosit+sumPtRegionNegat)/2.0);
 499
 500 }
 501
 502
 503 //____________________________________________________________________
 504 void AliAnalysisTaskUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
 505 {
 506    // Fill sumPt of control regions
 507
 508    // Max cone
 509    fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax );
 510    // Min cone
 511    fhRegionSumPtMinVsEt->Fill( leadingE, ptMin );
 512    // MAke distributions for UE comparison with MB data
 513    fhMinRegSumPt->Fill(ptMin);
 514 }
 515
 516 //____________________________________________________________________
 517 void AliAnalysisTaskUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
 518 {
 519    // Fill average particle Pt of control regions
 520
 521    // Max cone
 522    fhRegionAvePartPtMaxVsEt->Fill( leadingE, ptMax );
 523    // Min cone
 524    fhRegionAvePartPtMinVsEt->Fill( leadingE, ptMin );
 525    // MAke distributions for UE comparison with MB data
 526    fhMinRegAvePt->Fill(ptMin);
 527 }
 528
 529 //____________________________________________________________________
 530 void AliAnalysisTaskUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin  )
 531 {
 532    // Fill Nch multiplicity of control regions
 533
 534    // Max cone
 535    fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax );
 536    fhRegionMultMax->Fill( nTrackPtmax );
 537    // Min cone
 538    fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin );
 539    fhRegionMultMin->Fill( nTrackPtmin );
 540    // MAke distributions for UE comparison with MB data
 541    fhMinRegSumPtvsMult->Fill(nTrackPtmin,ptMin);
 542
 543 }
 544
 545 //____________________________________________________________________
 546 Int_t AliAnalysisTaskUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect)
 547 {
 548   // return de region in delta phi
 549   // -1 negative delta phi
 550   //  1 positive delta phi
 551   //  0 outside region
 552   static const Double_t k60rad  = 60.*TMath::Pi()/180.;
 553   static const Double_t k120rad = 120.*TMath::Pi()/180.;
 554
 555   Int_t region = 0;
 556   if( fRegionType == 1 ) {
 557     if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
 558     // transverse regions
 559     if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1;
 560     if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1;
 561
 562   } else if( fRegionType == 2 ) {
 563     // Cone regions
 564     Double_t deltaR = 0.;
 565
 566     TVector3 positVect,negatVect;
 567     positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
 568     negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
 569
 570     if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) {
 571        region = 1;
 572        deltaR = positVect.DrEtaPhi(*partVect);
 573     } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) {
 574        region = -1;
 575        deltaR = negatVect.DrEtaPhi(*partVect);
 576     }
 577
 578     if (deltaR > fConeRadius) region = 0;
 579
 580   } else
 581     AliError("Unknow region type");
 582
 583   // For debug (to be removed)
 584   //if( region != 0 )  fhValidRegion->Fill( partVect->Eta()-jetVect[0].Eta(), jetVect[0].DeltaPhi(*partVect) );
 585
 586   return region;
 587 }
 588
 589 //____________________________________________________________________
 590 void  AliAnalysisTaskUE::CreateHistos()
 591 {
 592   fListOfHistos = new TList();
 593
 594
 595   fhNJets = new TH1F("fhNJets", "n Jet",  10, 0, 10);
 596   fhNJets->SetXTitle("# of jets");
 597   fhNJets->Sumw2();
 598   fhNJets->Sumw2();
 599   fListOfHistos->Add( fhNJets );                 // At(0)
 600
 601   fhEleadingPt  = new TH1F("hEleadingPt",   "Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 602   fhEleadingPt->SetXTitle("P_{T} (GeV/c)");
 603   fhEleadingPt->SetYTitle("dN/dP_{T}");
 604   fhEleadingPt->Sumw2();
 605   fListOfHistos->Add( fhEleadingPt );            // At(1)
 606
 607   fhMinRegPtDist = new TH1F("hMinRegPtDist",   "P_{T} distribution in Min zone",  50,0.,20.);
 608   fhMinRegPtDist->SetXTitle("P_{T} (GeV/c)");
 609   fhMinRegPtDist->SetYTitle("dN/dP_{T}");
 610   fhMinRegPtDist->Sumw2();
 611   fListOfHistos->Add( fhMinRegPtDist );          // At(2)
 612
 613   fhRegionMultMin = new TH1F("hRegionMultMin",      "N_{ch}^{90, min}",  21, -0.5,   20.5);
 614   fhRegionMultMin->SetXTitle("N_{ch tracks}");
 615   fhRegionMultMin->Sumw2();
 616   fListOfHistos->Add( fhRegionMultMin );         // At(3)
 617
 618   fhMinRegAvePt = new TH1F("hMinRegAvePt", "#LTp_{T}#GT",  50, 0.,   20.);
 619   fhMinRegAvePt->SetXTitle("P_{T} (GeV/c)");
 620   fhMinRegAvePt->Sumw2();
 621   fListOfHistos->Add( fhMinRegAvePt );           // At(4)
 622
 623   fhMinRegSumPt = new TH1F("hMinRegSumPt", "#Sigma p_{T} ",  50, 0.,   20.);
 624   fhMinRegSumPt->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");
 625   fhMinRegSumPt->SetXTitle("#Sigma p_{T} (GeV/c)");
 626   fhMinRegSumPt->Sumw2();
 627   fListOfHistos->Add( fhMinRegSumPt );           // At(5)
 628
 629   fhMinRegMaxPtPart = new TH1F("hMinRegMaxPtPart", "max(p_{T})|_{event} ",  50, 0.,   20.);
 630   fhMinRegMaxPtPart->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");
 631   fhMinRegMaxPtPart->SetXTitle("p_{T} (GeV/c)");
 632   fhMinRegMaxPtPart->Sumw2();
 633   fListOfHistos->Add( fhMinRegMaxPtPart );       // At(6)
 634
 635   fhMinRegSumPtvsMult = new TH1F("hMinRegSumPtvsMult", "#Sigma p_{T} vs. Multiplicity ",  21, -0.5,   20.5);
 636   fhMinRegSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");
 637   fhMinRegSumPtvsMult->SetXTitle("N_{charge}");
 638   fhMinRegSumPtvsMult->Sumw2();
 639   fListOfHistos->Add( fhMinRegSumPtvsMult );     // At(7);
 640
 641   fhdNdEta_PhiDist  = new TH1F("hdNdEta_PhiDist",   "Charge particle density |#eta|< 1 vs #Delta#phi",  120, 0.,   2.*TMath::Pi());
 642   fhdNdEta_PhiDist->SetXTitle("#Delta#phi");
 643   fhdNdEta_PhiDist->SetYTitle("dN_{ch}/d#etad#phi");
 644   fhdNdEta_PhiDist->Sumw2();
 645   fListOfHistos->Add( fhdNdEta_PhiDist );        // At(8)
 646
 647   // Can be use to get part pt distribution for differente Jet Pt bins
 648   fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", "dN/dP_{T} |#eta|<1 vs Leading Jet P_{T}",
 649                50,0.,50., fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 650   fhFullRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
 651   fhFullRegPartPtDistVsEt->SetXTitle("p_{T}");
 652   fhFullRegPartPtDistVsEt->Sumw2();
 653   fListOfHistos->Add( fhFullRegPartPtDistVsEt );  // At(9)
 654
 655    // Can be use to get part pt distribution for differente Jet Pt bins
 656   fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", "dN/dP_{T} in tranvese regions |#eta|<1 vs Leading Jet P_{T}",
 657                50,0.,50., fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 658   fhTransRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
 659   fhTransRegPartPtDistVsEt->SetXTitle("p_{T}");
 660   fhTransRegPartPtDistVsEt->Sumw2();
 661   fListOfHistos->Add( fhTransRegPartPtDistVsEt );  // At(10)
 662
 663
 664   fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt",  "P_{T}^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 665   fhRegionSumPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
 666   fhRegionSumPtMaxVsEt->Sumw2();
 667   fListOfHistos->Add( fhRegionSumPtMaxVsEt );      // At(11)
 668
 669   fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt",   "P_{T}^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 670   fhRegionSumPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
 671   fhRegionSumPtMinVsEt->Sumw2();
 672   fListOfHistos->Add( fhRegionSumPtMinVsEt );      // At(12)
 673
 674   fhRegionMultMax = new TH1I("hRegionMultMax",      "N_{ch}^{90, max}",  21, -0.5,   20.5);
 675   fhRegionMultMax->SetXTitle("N_{ch tracks}");
 676   fhRegionMultMax->Sumw2();
 677   fListOfHistos->Add( fhRegionMultMax );           // At(13)
 678
 679   fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt",  "N_{ch}^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 680   fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)");
 681   fhRegionMultMaxVsEt->Sumw2();
 682   fListOfHistos->Add( fhRegionMultMaxVsEt );       // At(14)
 683
 684   fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt",  "N_{ch}^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 685   fhRegionMultMinVsEt->SetXTitle("E (GeV/c)");
 686   fhRegionMultMinVsEt->Sumw2();
 687   fListOfHistos->Add( fhRegionMultMinVsEt );      // At(15)
 688
 689   fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 690   fhRegionAveSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
 691   fhRegionAveSumPtVsEt->Sumw2();
 692   fListOfHistos->Add( fhRegionAveSumPtVsEt );     // At(16)
 693
 694   fhRegionDiffSumPtVsEt= new TH1F("hRegionPtDiffVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 695   fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
 696   fhRegionDiffSumPtVsEt->Sumw2();
 697   fListOfHistos->Add( fhRegionDiffSumPtVsEt );    // At(17)
 698
 699   fhRegionAvePartPtMaxVsEt = new TH1F("hRegionAvePartPtMaxVsEt", "#LTp_{T}#GT^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 700   fhRegionAvePartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
 701   fhRegionAvePartPtMaxVsEt->Sumw2();
 702   fListOfHistos->Add( fhRegionAvePartPtMaxVsEt );  // At(18)
 703
 704   fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 705   fhRegionAvePartPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
 706   fhRegionAvePartPtMinVsEt->Sumw2();
 707   fListOfHistos->Add( fhRegionAvePartPtMinVsEt );   // At(19)
 708
 709   fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
 710   fhRegionMaxPartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
 711   fhRegionMaxPartPtMaxVsEt->Sumw2();
 712   fListOfHistos->Add( fhRegionMaxPartPtMaxVsEt );    // At(20)
 713
 714 /*
 715   // For debug region selection
 716   fhValidRegion = new TH2F("hValidRegion", "dN/d#eta/d#phi",
 717                20, -2.,2., 62, -TMath::Pi(),   TMath::Pi());
 718   fhValidRegion->SetYTitle("#Delta#phi");
 719   fhValidRegion->Sumw2();
 720   fListOfHistos->Add( fhValidRegion );  // At(15)
 721 */
 722 }
 723
 724
 725