b42fec1bf3641c60447c096a3c1bf61ae91f1d4a
[u/mrichter/AliRoot.git] / 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 <TProfile.h>
27 #include <TCanvas.h>
28 #include <TVector3.h>
29 #include <TLorentzVector.h>
30 #include <TMath.h>
31 #include <TTree.h>
32 #include <TBranch.h>
33 #include <TRandom.h>
34
35 #include "AliAnalysisTaskUE.h"
36 #include "AliAnalysisManager.h"
37 #include "AliMCEventHandler.h"
38 #include "AliMCEvent.h"
39 #include "AliAODEvent.h"
40 #include "AliAODInputHandler.h"
41 #include "AliAODHandler.h"
42 #include "AliStack.h"
43 #include "AliAODJet.h"
44 #include "AliAODTrack.h"
45 #include "AliAODMCParticle.h"
46 #include "AliKFVertex.h"
47
48 #include "AliGenPythiaEventHeader.h"
49 #include "AliAnalysisHelperJetTasks.h"
50 #include "AliInputEventHandler.h"
51 #include "AliStack.h"
52 #include "AliLog.h"
53
54 //__________________________________________________________________
55 // Analysis class for Underlying Event studies
56 //
57 // Look for correlations on the tranverse regions to 
58 // the leading charged jet
59 //
60 // This class needs as input AOD with track and Jets.
61 // The output is a list of histograms
62 //
63 // AOD can be either connected to the InputEventHandler  
64 // for a chain of AOD files 
65 // or 
66 // to the OutputEventHandler
67 // for a chain of ESD files, so this case class should be 
68 // in the train after the Jet finder
69 //
70 //    Arian.Abrahantes.Quintana@cern.ch 
71 //    Ernesto.Lopez.Torres@cern.ch
72 // 
73
74
75 ClassImp( AliAnalysisTaskUE)
76
77 ////////////////////////////////////////////////////////////////////////
78
79
80 //____________________________________________________________________
81 AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name):
82 AliAnalysisTask(name,""),
83 fTrigger(0),
84 fDebug(0),
85 fDeltaAOD(kFALSE),
86 fDeltaAODBranch(""),
87 fAODBranch("jets"),
88 fArrayJets(0x0),           
89 fAOD(0x0),            
90 fAODjets(0x0),
91 fListOfHistos(0x0),  
92 fBinsPtInHist(30),     
93 fMinJetPtInHist(0.),
94 fMaxJetPtInHist(300.), 
95 fIsNorm2Area(kTRUE),
96 fUseMCParticleBranch(kFALSE),
97 fConstrainDistance(kTRUE),
98 fMinDistance(0.2),
99 fSimulateChJetPt(kFALSE),
100 fUseAliStack(kTRUE),
101 fnTracksVertex(3),  // QA tracks pointing to principal vertex (= 3 default) 
102 fZVertex(5.),
103 fAnaType(1),         
104 fRegionType(1),
105 fConeRadius(0.7),
106 fConePosition(1),
107 fAreaReg(1.5393), // Pi*0.7*0.7
108 fUseChPartJet(kFALSE),
109 fUseChPartMaxPt(kFALSE),
110 fUseChargeHadrons(kFALSE),
111 fUseSingleCharge(kFALSE),
112 fUsePositiveCharge(kTRUE),
113 fOrdering(1),
114 fFilterBit(0xFF),
115 fJetsOnFly(kFALSE),
116 fChJetPtMin(5.0),
117 fJet1EtaCut(0.2),
118 fJet2DeltaPhiCut(2.616),    // 150 degrees
119 fJet2RatioPtCut(0.8),
120 fJet3PtCut(15.),
121 fTrackPtCut(0.),
122 fTrackEtaCut(0.9),
123 fAvgTrials(1),
124 fhNJets(0x0),
125 fhEleadingPt(0x0),
126 fhMinRegPtDist(0x0),
127 fhRegionMultMin(0x0),
128 fhMinRegAvePt(0x0), 
129 fhMinRegSumPt(0x0),            
130 fhMinRegMaxPtPart(0x0),
131 fhMinRegSumPtvsMult(0x0),
132 fhdNdEtaPhiDist(0x0),        
133 fhFullRegPartPtDistVsEt(0x0), 
134 fhTransRegPartPtDistVsEt(0x0),
135 fhRegionSumPtMaxVsEt(0x0),
136 fhRegionMultMax(0x0),         
137 fhRegionMultMaxVsEt(0x0),     
138 fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0),         
139 fhRegionMultMinVsEt(0x0),
140 fhRegionAveSumPtVsEt(0x0),
141 fhRegionDiffSumPtVsEt(0x0),
142 fhRegionAvePartPtMaxVsEt(0x0),
143 fhRegionAvePartPtMinVsEt(0x0),
144 fhRegionMaxPartPtMaxVsEt(0x0),
145 //fhRegForwardSumPtVsEt(0x0),
146 //fhRegForwardMultVsEt(0x0),
147 //fhRegBackwardSumPtVsEt(0x0),
148 //fhRegBackwardMultVsEt(0x0),
149 fhRegForwardMult(0x0),
150 fhRegForwardSumPtvsMult(0x0),
151 fhRegBackwardMult(0x0),
152 fhRegBackwardSumPtvsMult(0x0),
153 fhRegForwardPartPtDistVsEt(0x0),
154 fhRegBackwardPartPtDistVsEt(0x0),
155 fhRegTransMult(0x0),
156 fhRegTransSumPtVsMult(0x0),
157 fhMinRegSumPtJetPtBin(0x0),
158 fhMaxRegSumPtJetPtBin(0x0),
159 fhVertexMult(0x0),
160 fh1Xsec(0x0),
161 fh1Trials(0x0),
162 fSettingsTree(0x0)//,   fhValidRegion(0x0)
163 {
164   // Default constructor
165   // Define input and output slots here
166   // Input slot #0 works with a TChain
167   DefineInput(0, TChain::Class());
168   // Output slot #0 writes into a TList container
169   DefineOutput(0, TList::Class());
170
171 }
172
173 //______________________________________________________________
174 Bool_t AliAnalysisTaskUE::Notify()
175 {
176   //
177   // Implemented Notify() to read the cross sections
178   // and number of trials from pyxsec.root
179   // Copy from AliAnalysisTaskJFSystematics
180   fAvgTrials = 1;
181   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
182   Float_t xsection = 0;
183   Float_t ftrials  = 1;
184   if(tree){
185     TFile *curfile = tree->GetCurrentFile();
186     if (!curfile) {
187       Error("Notify","No current file");
188       return kFALSE;
189     }
190     if(!fh1Xsec||!fh1Trials){
191       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
192       return kFALSE;
193     }
194     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials); 
195     fh1Xsec->Fill("<#sigma>",xsection);
196
197    // construct average trials
198    Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
199    if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
200   }
201   
202   return kTRUE;
203 }
204
205 //____________________________________________________________________
206 void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/)
207 {
208   // Connect the input data  
209   
210   // We need AOD with tracks and jets.
211   // Since AOD can be either connected to the InputEventHandler (input chain fron AOD files)
212   // or to the OutputEventHandler ( AOD is create by a previus task in the train )
213   // we need to check where it is and get the pointer to AODEvent in the right way
214   
215   // Delta AODs are also implemented
216   
217  
218   if (fDebug > 1) AliInfo("ConnectInputData() ");
219   
220   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
221   
222   if( handler && handler->InheritsFrom("AliAODInputHandler") ) { // input AOD
223     fAOD = ((AliAODInputHandler*)handler)->GetEvent();
224     if (fDebug > 1) AliInfo(" ==== Tracks from AliAODInputHandler");
225     // Case when jets are reconstructed on the fly from AOD tracks
226     // (the Jet Finder is using the AliJetAODReader) of InputEventHandler
227     // and put in the OutputEventHandler AOD. Useful whe you want to reconstruct jets with
228     // different parameters to default ones stored in the AOD or to use a different algorithm
229     if( fJetsOnFly ) {
230       handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
231       if( handler && handler->InheritsFrom("AliAODHandler") ) {
232         fAODjets = ((AliAODHandler*)handler)->GetAOD();
233         if (fDebug > 1) AliInfo(" ==== Jets from AliAODHandler (on the fly)");
234       }
235     } else {
236       fAODjets = fAOD;
237       if (fDebug > 1) AliInfo(" ==== Jets from AliAODInputHandler");
238     }
239   } else {  //output AOD
240     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
241     if( handler && handler->InheritsFrom("AliAODHandler") ) {
242       fAOD = ((AliAODHandler*)handler)->GetAOD();
243       fAODjets = fAOD;
244       if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
245     } else {
246       AliFatal("I can't get any AOD Event Handler");
247       return;
248     }
249   }
250 }
251
252 //____________________________________________________________________
253 void  AliAnalysisTaskUE::CreateOutputObjects()
254 {
255   // Create the output container
256   //
257   if (fDebug > 1) AliInfo("CreateOutPutData()");
258   //
259   //  Histograms
260
261   // OpenFile(0);
262   CreateHistos();
263   fListOfHistos->SetOwner(kTRUE);  
264   PostData(0,fListOfHistos);
265
266 }
267
268 //____________________________________________________________________
269 void  AliAnalysisTaskUE::Exec(Option_t */*option*/)
270 {
271   //Trigger selection ************************************************
272   /*AliInputEventHandler* inputHandler = (AliInputEventHandler*)
273          ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
274   if (inputHandler->IsEventSelected()) {
275     if (fDebug > 1) AliInfo(" Trigger Selection: event ACCEPTED ... ");
276   } else {
277     if (fDebug > 1) AliInfo(" Trigger Selection: event REJECTED ... ");
278     return;
279   }
280    */                                     
281   //Event selection (vertex) *****************************************
282
283   Int_t nVertex = fAOD->GetNumberOfVertices();
284   if( nVertex > 0 ) { // Only one vertex (reject pileup??)
285      AliAODVertex* vertex = (AliAODVertex*)fAOD->GetPrimaryVertex();
286      TString tpcvertex("TPCVertex");
287      if (vertex->GetName() == tpcvertex){
288         if (fDebug > 1) AliInfo(Form("Primary vertex selection: %s event REJECTED ...",vertex->GetName()));
289         return;
290         }
291      Int_t nTracksPrim = vertex->GetNContributors();
292      Double_t zVertex = vertex->GetZ();
293      if (fDebug > 1) AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
294      // Select a quality vertex by number of tracks?
295      if( nTracksPrim < fnTracksVertex || TMath::Abs(zVertex) > fZVertex ) {
296         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
297         return;
298      }
299      if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
300   } else {
301      if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
302      return;
303   } 
304   /*AliKFVertex primVtx(*(fAOD->GetPrimaryVertex()));
305   Int_t nTracksPrim=primVtx.GetNContributors();
306   if (fDebug > 1) AliInfo(Form(" Primary-vertex Selection: %d",nTracksPrim));
307   if(!nTracksPrim){
308     if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
309     return;
310   }
311   if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED ...");
312   */
313   // Execute analysis for current event
314   //
315   if ( fDebug > 3 ) AliInfo( " Processing event..." );
316   // fetch the pythia header info and get the trials
317   AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
318   Float_t nTrials = 1;
319   if (mcHandler) {
320        AliMCEvent* mcEvent = mcHandler->MCEvent();
321        if (mcEvent) {
322          AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
323          if(pythiaGenHeader){
324            nTrials = pythiaGenHeader->Trials();
325          }
326       }
327     }
328   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
329   AnalyseUE();
330
331   // Post the data
332   PostData(0, fListOfHistos);
333 }
334
335 //____________________________________________________________________
336 void  AliAnalysisTaskUE::AnalyseUE()
337 {
338    Double_t const kMyTolerance = 0.0000001; 
339   //
340   // Look for correlations on the tranverse regions to 
341   // the leading charged jet
342   // 
343   
344   // ------------------------------------------------
345   // Find Leading Jets 1,2,3 
346   // (could be skipped if Jets are sort by Pt...)
347   Double_t maxPtJet1 = 0.; 
348   Int_t    index1 = -1;
349   Double_t maxPtJet2 = 0.;   // jet 2 need for back to back inclusive
350   Int_t    index2 = -1;
351   Double_t maxPtJet3 = 0.;   // jet 3 need for back to back exclusive
352   Int_t    index3 = -1;
353   TVector3 jetVect[3];
354   Int_t nJets = 0;
355   
356   
357   if( !fUseChPartJet && fAnaType != 4 ) {
358     
359     // Use AOD Jets
360     if(fDeltaAOD){
361       if (fDebug > 1) AliInfo(" ==== Jets From  Delta-AODs !");
362       if (fDebug > 1) AliInfo(Form(" ====  Reading Branch: %s  ",fDeltaAODBranch.Data()));
363            fArrayJets = (TClonesArray*)fAODjets->GetList()->FindObject(fDeltaAODBranch.Data());
364     } else {
365       if (fDebug > 1) AliInfo(" ==== Read Standard-AODs  !");
366       if (fDebug > 1) AliInfo(Form(" ====  Reading Branch: %s  ",fAODBranch.Data()));
367       fArrayJets = (TClonesArray*)fAODjets->GetList()->FindObject(fAODBranch.Data());
368     }
369     if (!fArrayJets){
370        AliFatal(" No jet-array! ");
371        return;
372     }
373     nJets=fArrayJets->GetEntries();
374     //printf("AOD %d jets \n", nJets);
375
376     for( Int_t i=0; i<nJets; ++i ) {
377       AliAODJet* jet = (AliAODJet*)fArrayJets->At(i);
378       Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!!
379  
380       if( jetPt > maxPtJet1 ) {
381              maxPtJet3 = maxPtJet2; index3 = index2;
382              maxPtJet2 = maxPtJet1; index2 = index1;
383              maxPtJet1 = jetPt; index1 = i;
384       } else if( jetPt > maxPtJet2 ) {
385              maxPtJet3 = maxPtJet2; index3 = index2;
386              maxPtJet2 = jetPt; index2 = i;
387       } else if( jetPt > maxPtJet3 ) {
388              maxPtJet3 = jetPt; index3 = i;
389       }
390     }
391
392     if( index1 != -1 ) {
393       AliAODJet *jet =(AliAODJet*) fArrayJets->At(index1);
394       if(jet)jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
395     }
396     if( index2 != -1 ) {
397       AliAODJet* jet= (AliAODJet*) fArrayJets->At(index2);
398       if(jet)jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
399     }
400     if( index3 != -1 ) {
401        AliAODJet* jet = (AliAODJet*) fArrayJets->At(index3);
402       if(jet)jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
403     }
404     
405   }
406
407   if( fUseChPartJet )  {
408
409     // Use "Charged Particle Jets"
410     TObjArray* jets = FindChargedParticleJets();
411     if( jets ) {
412       nJets = jets->GetEntriesFast();
413       if( nJets > 0 ) {
414              index1 = 0;
415              AliAODJet* jet = (AliAODJet*)jets->At(0);
416              maxPtJet1 = jet->Pt();
417              jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
418       }
419       if( nJets > 1 ) {
420              index2 = 1;
421              AliAODJet* jet = (AliAODJet*)jets->At(1);
422         maxPtJet2 = jet->Pt();
423              jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
424       }
425       if( nJets > 2 ) {
426              index3 = 2;
427              AliAODJet* jet = (AliAODJet*)jets->At(2);
428         maxPtJet3 = jet->Pt();
429              jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
430       }
431       
432       jets->Delete();
433       delete jets;
434       if (fDebug > 1) AliInfo(" ==== Jets Created in OUR class  !: CDF algorithm ");
435     }
436   }
437
438
439   if( fAnaType == 4 )  {
440
441     // Use "Max Pt Charge Particle"
442     TObjArray* tracks = SortChargedParticles();
443     if( tracks ) {
444       nJets = tracks->GetEntriesFast();
445       if( nJets > 0 ) {
446         index1 = 0;
447         AliAODTrack* jet = (AliAODTrack*)tracks->At(0);
448         maxPtJet1 = jet->Pt();
449         jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
450       }
451       tracks->Clear();
452       delete tracks; 
453     }
454   }
455
456
457   fhNJets->Fill(nJets);
458   
459   if( fDebug > 1 ) {
460     if( index1 < 0 ) {
461       AliInfo("\n   Skipping Event, not jet found...");
462       return;
463     } else
464       AliInfo(Form("\n   Pt Leading Jet = %6.1f eta=%5.3f ",  maxPtJet1, jetVect[0].Eta() ));
465   }
466   
467   
468   // ----------------------------------------------
469   // Cut events by jets topology
470   // fAnaType:
471   //     1 = inclusive,
472   //         - Jet1 |eta| < fJet1EtaCut
473   //     2 = back to back inclusive
474   //         - fulfill case 1
475   //         - |Jet1.Phi - Jet2.Phi| > fJet2DeltaPhiCut
476   //         - Jet2.Pt/Jet1Pt > fJet2RatioPtCut
477   //     3 = back to back exclusive
478   //         - fulfill case 2
479   //         - Jet3.Pt < fJet3PtCut
480   
481   if( index1 < 0 || TMath::Abs(jetVect[0].Eta()) > fJet1EtaCut) {
482     if( fDebug > 1 ) AliInfo("\n   Skipping Event...Jet1 |eta| > fJet1EtaCut");
483     return;
484   }
485   // back to back inclusive
486   if( fAnaType > 1 && fAnaType < 4 && index2 == -1 ) {
487     if( fDebug > 1 ) AliInfo("\n   Skipping Event... no second Jet found");
488     return;
489   }
490   if( fAnaType > 1 && fAnaType < 4 && index2 > -1 ) {
491     if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
492        maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) {
493       if( fDebug > 1 ) AliInfo("\n   Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
494       return;
495     }
496   }
497   // back to back exclusive
498   if( fAnaType > 2 && fAnaType < 4 && index3 > -1 ) {
499     if( maxPtJet3 > fJet3PtCut ) {
500       if( fDebug > 1 ) AliInfo("\n   Skipping Event... Jet3.Pt > fJet3PtCut ");
501       return;
502     }
503   }
504   
505   //fhEleadingPt->Fill( maxPtJet1 );
506
507   // ----------------------------------------------
508   // Find max and min regions
509   Double_t sumPtRegionPosit = 0.;
510   Double_t sumPtRegionNegat = 0.;
511   Double_t sumPtRegionForward = 0;
512   Double_t sumPtRegionBackward = 0;
513   Double_t maxPartPtRegion  = 0.;
514   Int_t    nTrackRegionPosit = 0;
515   Int_t    nTrackRegionNegat = 0;
516   Int_t    nTrackRegionForward = 0;
517   Int_t    nTrackRegionBackward = 0;
518   static Double_t const  kPI     = TMath::Pi();
519   static Double_t const  kTWOPI  = 2.*kPI;
520   static Double_t const  k270rad = 270.*kPI/180.;
521   static Double_t const  k120rad = 120.*kPI/180.;
522   
523   //Area for Normalization Purpose at Display histos
524   // Forward and backward
525   Double_t normArea = 1.;
526   // Transverse
527   if (fIsNorm2Area) {
528     SetRegionArea(jetVect);
529     normArea =  2.*fTrackEtaCut*k120rad ;
530   } else fAreaReg = 1.;
531   
532   if (!fUseMCParticleBranch){
533     AliAODVertex* vertex = (AliAODVertex*)fAOD->GetPrimaryVertex();
534     fhVertexMult->Fill( vertex->GetNContributors() );
535     
536     fhEleadingPt->Fill( maxPtJet1 );
537     Int_t nTracks = fAOD->GetNTracks();
538     if (fDebug > 1) AliInfo(Form(" ==== AOD tracks = %d \n ",nTracks));
539     
540     for (Int_t ipart=0; ipart<nTracks; ++ipart) {
541       
542       AliAODTrack* part = fAOD->GetTrack( ipart );
543       if (fDebug > 1) AliInfo(Form(" ==== AOD track = %d pt = %f charge = %d \n ",ipart,part->Pt(),part->Charge()));
544       if ( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
545       if (!part->IsPrimaryCandidate()) continue; // reject whatever is not linked to collision point
546       // PID Selection: Reject everything but hadrons
547       Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion || 
548                         part->GetMostProbablePID()==AliAODTrack::kKaon || 
549                         part->GetMostProbablePID()==AliAODTrack::kProton;
550       if ( fUseChargeHadrons && !isHadron ) continue;
551       
552       if ( !part->Charge() ) continue; //Only charged
553       if ( fUseSingleCharge ) { // Charge selection
554         if ( fUsePositiveCharge && part->Charge() < 0.) continue; // keep Positives
555         if ( !fUsePositiveCharge && part->Charge() > 0.) continue; // keep Negatives
556       }
557       
558       if ( part->Pt() < fTrackPtCut ) continue;
559       if( TMath::Abs(part->Eta()) > fTrackEtaCut ) continue;
560       
561       TVector3 partVect(part->Px(), part->Py(), part->Pz());
562       Bool_t isFlagPart = kTRUE;
563       Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
564       if( deltaPhi > kTWOPI )  deltaPhi-= kTWOPI;
565       if (fAnaType != 4 ) fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
566       else if (TMath::Abs(deltaPhi-k270rad) >= kMyTolerance && TMath::Abs(jetVect[0].Eta()-partVect.Eta()) >= kMyTolerance){
567          fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
568          isFlagPart = kFALSE;
569       }
570       
571       fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
572       
573       Int_t region = IsTrackInsideRegion( jetVect, &partVect );  
574       
575       if (region == 1) {
576         if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
577         sumPtRegionPosit += part->Pt();
578         nTrackRegionPosit++;
579         fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
580       }
581       if (region == -1) {
582         if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
583         sumPtRegionNegat += part->Pt();
584         nTrackRegionNegat++;
585         fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
586       }
587       if (region == 2){ //forward
588         sumPtRegionForward += part->Pt();
589         nTrackRegionForward++;
590         fhRegForwardPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
591       }
592       if (region == -2){ //backward
593         sumPtRegionBackward += part->Pt();
594         nTrackRegionBackward++;
595         fhRegBackwardPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
596       }
597     }//end loop AOD tracks
598   }
599   else {
600     
601     // this is the part we only use when we have MC information
602     // More than a test for values of it also resumes the reconstruction efficiency of jets
603     // As commented bellow if available for the data, we try to pair reconstructed jets with simulated ones
604     // afterwards we kept angular variables of MC jet to perform UE analysis over MC particles
605     // TODO: Handle Multiple jet environment. 06/2009 just suited for inclusive jet condition ( fAnaType = 1 ) 
606     
607       AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
608       if (!mcHandler) {
609         Printf("ERROR: Could not retrieve MC event handler");
610         return;
611       }
612       
613       AliMCEvent* mcEvent = mcHandler->MCEvent();
614       if (!mcEvent) {
615         Printf("ERROR: Could not retrieve MC event");
616         return;
617       }
618     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
619     if(!pythiaGenHeader){
620       return;
621     }
622     
623     //Get Jets from MC header
624     Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
625     AliAODJet pythiaGenJets[4];
626     TVector3 jetVectnew[4];
627     Int_t iCount = 0;
628     for(int ip = 0;ip < nPythiaGenJets;++ip){
629       if (iCount>3) break;
630       Float_t p[4];
631       pythiaGenHeader->TriggerJet(ip,p);
632       TVector3 tempVect(p[0],p[1],p[2]);
633       if ( TMath::Abs(tempVect.Eta())>fJet1EtaCut ) continue;
634       pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
635       jetVectnew[iCount].SetXYZ(pythiaGenJets[iCount].Px(), pythiaGenJets[iCount].Py(), pythiaGenJets[iCount].Pz());
636       iCount++;
637     }
638     
639     if (!iCount) return;// no jet in eta acceptance
640     
641     //Search the index of the nearest MC jet to the leading jet reconstructed from the input data
642     Int_t index = 0;
643     if (fConstrainDistance){
644       Float_t deltaR = 0.;
645       Float_t dRTemp = 0.;
646       for (Int_t i=0; i<iCount; i++){
647          if (!i) {
648          dRTemp = jetVectnew[i].DeltaR(jetVect[0]);
649          index = i;
650          }
651          deltaR = jetVectnew[i].DeltaR(jetVect[0]);
652          if (deltaR < dRTemp){
653          index = i;
654          dRTemp = deltaR;
655          }
656       }
657    
658       if (jetVectnew[index].DeltaR(jetVect[0]) > fMinDistance) return;
659     }
660     //Let's add some taste to jet and simulate pt of charged alone 
661     //eta and phi are kept as original
662     //Play a Normal Distribution
663     Float_t random = 1.;  
664     if (fSimulateChJetPt){
665       while(1){
666         random = gRandom->Gaus(0.6,0.25);
667         if (random > 0. && random < 1. && 
668             (random * jetVectnew[index].Pt()>6.)) break;
669       }
670     }
671     
672     //Set new Pt & Fill histogram accordingly
673     maxPtJet1 = random * jetVectnew[index].Pt();
674     
675     
676     fhEleadingPt->Fill( maxPtJet1 );
677
678     if (fUseAliStack){//Try Stack Information to perform UE analysis
679     
680       AliStack* mcStack = mcEvent->Stack();//Load Stack
681       Int_t nTracksMC = mcStack->GetNtrack();
682       for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
683         //Cuts
684         if(!(mcStack->IsPhysicalPrimary(iTracks))) continue;
685         
686         TParticle* mctrk = mcStack->Particle(iTracks);
687         
688         Double_t charge = mctrk->GetPDG()->Charge();
689         if (charge == 0) continue;
690         
691         if ( fUseSingleCharge ) { // Charge selection
692           if ( fUsePositiveCharge && charge < 0.) continue; // keep Positives
693           if ( !fUsePositiveCharge && charge > 0.) continue; // keep Negatives
694         }
695         
696         //Kinematics cuts on particle
697         if ((mctrk->Pt() < fTrackPtCut) || (TMath::Abs(mctrk->Eta()) > fTrackEtaCut )) continue;
698         
699         Bool_t isHadron = TMath::Abs(mctrk->GetPdgCode())==211 ||
700                           TMath::Abs(mctrk->GetPdgCode())==2212 ||
701                           TMath::Abs(mctrk->GetPdgCode())==321;
702         
703         if ( fUseChargeHadrons && !isHadron ) continue;
704         
705         TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
706
707         Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
708         if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
709         fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
710         
711         fhFullRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
712         
713         //We are not interested on stack organization but don't loose track of info
714         TVector3 tempVector =  jetVectnew[0];
715         jetVectnew[0] = jetVectnew[index];
716         jetVectnew[index] = tempVector;
717         
718         Int_t region = IsTrackInsideRegion( jetVectnew, &partVect );  
719         
720         if (region == 1) {
721           if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
722           sumPtRegionPosit += mctrk->Pt();
723           nTrackRegionPosit++;
724           fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
725         }
726         if (region == -1) {
727           if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
728           sumPtRegionNegat += mctrk->Pt();
729           nTrackRegionNegat++;
730           fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
731         }
732         if (region == 2){ //forward
733           sumPtRegionForward += mctrk->Pt();
734           nTrackRegionForward++;
735           fhRegForwardPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
736         }
737         if (region == -2){ //backward
738           sumPtRegionBackward += mctrk->Pt();
739           nTrackRegionBackward++;
740           fhRegBackwardPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
741         }
742         
743       }  // end loop stack Particles
744       
745     }else{//Try mc Particle
746
747       TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
748        
749       Int_t ntrks = farray->GetEntries();
750       if (fDebug>1) AliInfo(Form("In UE MC analysis tracks %d \n",ntrks));
751       for(Int_t i =0 ; i < ntrks; i++){   
752         AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
753         //Cuts
754         if (!(mctrk->IsPhysicalPrimary())) continue;
755         //if (!(mctrk->IsPrimary())) continue;
756         
757         if (mctrk->Charge() == 0 || mctrk->Charge()==-99) continue;
758         
759         if (mctrk->Pt() < fTrackPtCut ) continue;
760         if( TMath::Abs(mctrk->Eta()) > fTrackEtaCut ) continue;
761         
762         Bool_t isHadron = TMath::Abs(mctrk->GetPdgCode())==211 ||
763                           TMath::Abs(mctrk->GetPdgCode())==2212 ||
764                           TMath::Abs(mctrk->GetPdgCode())==321;
765         
766         if ( fUseChargeHadrons && !isHadron ) continue;
767         
768         TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
769
770         Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
771         if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
772         fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
773
774         fhFullRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
775         
776         //We are not interested on stack organization but don't loose track of info
777         TVector3 tempVector =  jetVectnew[0];
778         jetVectnew[0] = jetVectnew[index];
779         jetVectnew[index] = tempVector;
780         
781         Int_t region = IsTrackInsideRegion( jetVectnew, &partVect );  
782         
783         if (region == 1) { //right
784           if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
785           sumPtRegionPosit += mctrk->Pt();
786           nTrackRegionPosit++;
787           fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
788         }
789         if (region == -1) { //left
790           if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
791           sumPtRegionNegat += mctrk->Pt();
792           nTrackRegionNegat++;
793           fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
794         }
795         if (region == 2){ //forward
796           sumPtRegionForward += mctrk->Pt();
797           nTrackRegionForward++;
798           fhRegForwardPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
799         }
800         if (region == -2){ //backward
801           sumPtRegionBackward += mctrk->Pt();
802           nTrackRegionBackward++;
803           fhRegBackwardPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
804         }
805         
806       }//end loop AliAODMCParticle tracks
807     }
808   }
809   
810   Double_t avePosRegion = (nTrackRegionPosit) ? sumPtRegionPosit/nTrackRegionPosit : 0.;
811   Double_t aveNegRegion = (nTrackRegionNegat) ? sumPtRegionNegat/nTrackRegionNegat : 0.;
812   if( avePosRegion > aveNegRegion ) {
813      FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg );
814   } else {
815      FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg );
816   }
817   
818   //How quantities will be sorted before Fill Min and Max Histogram
819   //  1=Plots will be CDF-like
820   //  2=Plots will be Marchesini-like
821   //  3=Minimum zone is selected as the one having lowest pt per track 
822   if( fOrdering == 1 ) {
823     if( sumPtRegionPosit > sumPtRegionNegat ) {
824       FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
825     } else {
826       FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
827     }
828     if (nTrackRegionPosit > nTrackRegionNegat ) {
829       FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
830     } else {
831       FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
832     }
833   } else if( fOrdering == 2 ) {
834     if (sumPtRegionPosit > sumPtRegionNegat) {
835       FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
836       FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
837     } else {
838       FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
839       FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
840     }
841   } else if( fOrdering == 3 ){
842      if (avePosRegion > aveNegRegion) {
843         FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
844         FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
845      }else{
846         FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
847         FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
848      }
849   }
850   fhRegionMaxPartPtMaxVsEt->Fill(maxPtJet1, maxPartPtRegion );
851   
852   // Compute pedestal like magnitudes
853   fhRegionDiffSumPtVsEt->Fill(maxPtJet1, TMath::Abs(sumPtRegionPosit-sumPtRegionNegat)/(2.0*fAreaReg));
854   fhRegionAveSumPtVsEt->Fill(maxPtJet1, (sumPtRegionPosit+sumPtRegionNegat)/(2.0*fAreaReg));
855   // Transverse as a whole
856   fhRegTransMult->Fill( maxPtJet1, nTrackRegionPosit + nTrackRegionNegat, (nTrackRegionPosit + nTrackRegionNegat)/(2.0*fAreaReg));
857   fhRegTransSumPtVsMult->Fill(maxPtJet1, nTrackRegionPosit + nTrackRegionNegat , (sumPtRegionNegat + sumPtRegionPosit)/(2.0*fAreaReg) );
858   
859   // Fill Histograms for Forward and away side w.r.t. leading jet direction
860   // Pt dependence
861   //fhRegForwardSumPtVsEt->Fill( maxPtJet1, sumPtRegionForward/normArea );
862   //fhRegForwardMultVsEt->Fill( maxPtJet1, nTrackRegionForward/normArea );
863   //fhRegBackwardSumPtVsEt->Fill( maxPtJet1, sumPtRegionBackward/normArea );
864   //fhRegBackwardMultVsEt->Fill( maxPtJet1, nTrackRegionBackward/normArea );
865   // Multiplicity dependence
866   fhRegForwardMult->Fill(maxPtJet1, nTrackRegionForward, nTrackRegionForward/normArea);
867   fhRegForwardSumPtvsMult->Fill(maxPtJet1, nTrackRegionForward,sumPtRegionForward/normArea);
868   fhRegBackwardMult->Fill(maxPtJet1, nTrackRegionBackward, nTrackRegionBackward/normArea );
869   fhRegBackwardSumPtvsMult->Fill(maxPtJet1, nTrackRegionBackward,sumPtRegionBackward/normArea);
870   
871 }
872
873 //____________________________________________________________________
874 void AliAnalysisTaskUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
875 {
876   // Fill sumPt of control regions
877   
878   // Max cone
879   fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax );
880   // Min cone
881   fhRegionSumPtMinVsEt->Fill( leadingE, ptMin );
882   // MAke distributions for UE comparison with MB data
883   fhMinRegSumPt->Fill(ptMin);
884   fhMinRegSumPtJetPtBin->Fill(leadingE, ptMin);
885   fhMaxRegSumPtJetPtBin->Fill(leadingE, ptMax);
886 }
887
888 //____________________________________________________________________
889 void AliAnalysisTaskUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
890 {
891   // Fill average particle Pt of control regions
892   
893   // Max cone
894   fhRegionAvePartPtMaxVsEt->Fill( leadingE, ptMax );
895   // Min cone
896   fhRegionAvePartPtMinVsEt->Fill( leadingE, ptMin );
897   // MAke distributions for UE comparison with MB data
898   fhMinRegAvePt->Fill(ptMin);
899 }
900
901 //____________________________________________________________________
902 void AliAnalysisTaskUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin  )
903 {
904   // Fill Nch multiplicity of control regions
905   
906   // Max cone
907   fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax );
908   fhRegionMultMax->Fill( nTrackPtmax );
909   // Min cone
910   fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin );
911   fhRegionMultMin->Fill( nTrackPtmin );
912   // MAke distributions for UE comparison with MB data
913   fhMinRegSumPtvsMult->Fill(nTrackPtmin,ptMin);
914 }
915
916 //____________________________________________________________________
917 Int_t AliAnalysisTaskUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect) 
918 {  
919   // return de region in delta phi
920   // -1 negative delta phi 
921   //  1 positive delta phi
922   //  0 outside region
923   static const Double_t k60rad  = 60.*TMath::Pi()/180.;
924   static const Double_t k120rad = 120.*TMath::Pi()/180.;
925   
926   Int_t region = 0;
927   if( fRegionType == 1 ) {
928     if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
929     // transverse regions
930     if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1; //left
931     if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1;    //right
932     if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) < k60rad ) region = 2;    //forward
933     if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) > k120rad ) region = -2; //backward
934     
935   } else if( fRegionType == 2 ) {
936     // Cone regions
937     Double_t deltaR = 0.;
938     
939     TVector3 positVect,negatVect;
940     if (fConePosition==1){
941       positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
942       negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
943     }else if (fConePosition==2){
944        if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
945        positVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()+TMath::PiOver2());
946        negatVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()-TMath::PiOver2());
947     }else if (fConePosition==3){
948        if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
949        Double_t weightEta = jetVect[0].Eta() * jetVect[0].Pt()/(jetVect[0].Pt() + jetVect[1].Pt()) + 
950                             jetVect[1].Eta() * jetVect[1].Pt()/(jetVect[0].Pt() + jetVect[1].Pt());
951        //Double_t weightEta = jetVect[0].Eta() * jetVect[0].Mag()/(jetVect[0].Mag() + jetVect[1].Mag()) + 
952        //                     jetVect[1].Eta() * jetVect[1].Mag()/(jetVect[0].Mag() + jetVect[1].Mag());
953        positVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()+TMath::PiOver2());
954        negatVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()-TMath::PiOver2());
955     }
956     if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) { 
957       region = 1;  
958       deltaR = positVect.DrEtaPhi(*partVect);
959     } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) { 
960       region = -1;  
961       deltaR = negatVect.DrEtaPhi(*partVect);
962     }
963     
964     if (deltaR > fConeRadius) region = 0;
965     
966   } else 
967     AliError("Unknow region type");
968   
969   // For debug (to be removed)
970   if( fDebug > 5 && region != 0 ) AliInfo(Form("Delta Phi = %3.2f region = %d \n", jetVect[0].DeltaPhi(*partVect), region)); //fhValidRegion->Fill( partVect->Eta()-jetVect[0].Eta(), jetVect[0].DeltaPhi(*partVect) );
971   
972   return region;
973 }
974
975 //____________________________________________________________________
976 TObjArray*  AliAnalysisTaskUE::SortChargedParticles()
977 {
978   //  return an array with all charged particles ordered according to their pT .
979   Int_t nTracks = fAOD->GetNTracks();
980   if( !nTracks ) return 0;
981   TObjArray* tracks = new TObjArray(nTracks);
982
983   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
984     AliAODTrack* part = fAOD->GetTrack( ipart );
985     if( !part->TestFilterBit( fFilterBit ) ) continue; // track cut selection
986     if( !part->Charge() ) continue;
987     if( part->Pt() < fTrackPtCut ) continue;
988     tracks->AddLast( part );
989   }
990   QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
991
992   nTracks = tracks->GetEntriesFast();
993   if( !nTracks ) return 0;
994
995   return tracks;
996 }
997
998 //____________________________________________________________________
999 TObjArray*  AliAnalysisTaskUE::FindChargedParticleJets()
1000 {
1001   // Return a TObjArray of "charged particle jets"
1002   //
1003   // Charged particle jet definition from reference:
1004   //
1005   // "Charged jet evolution and the underlying event
1006   //  in proton-antiproton collisions at 1.8 TeV"
1007   //  PHYSICAL REVIEW D 65 092002, CDF Collaboration
1008   //
1009   // We defined "jets" as circular regions in eta-phi space with
1010   // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
1011   // Our jet algorithm is as follows:
1012   //   1- Order all charged particles according to their pT .
1013   //   2- Start with the highest pT particle and include in the jet all
1014   //      particles within the radius R=0.7 considering each particle
1015   //      in the order of decreasing pT and recalculating the centroid
1016   //      of the jet after each new particle is added to the jet .
1017   //   3- Go to the next highest pT particle not already included in
1018   //      a jet and add to the jet all particles not already included in
1019   //      a jet within R=0.7.
1020   //   4- Continue until all particles are in a jet.
1021   // We defined the transverse momentum of the jet to be
1022   // the scalar pT sum of all the particles within the jet, where pT
1023   // is measured with respect to the beam axis
1024   
1025   //  1 - Order all charged particles according to their pT .
1026   Int_t nTracks = fAOD->GetNTracks();
1027   if( !nTracks ) return 0;
1028   TObjArray tracks(nTracks);
1029   
1030   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1031     AliAODTrack* part = fAOD->GetTrack( ipart );
1032     if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
1033     if( !part->Charge() ) continue;
1034     if( part->Pt() < fTrackPtCut ) continue;
1035     tracks.AddLast(part);
1036   }
1037   QSortTracks( tracks, 0, tracks.GetEntriesFast() );
1038   
1039   nTracks = tracks.GetEntriesFast();
1040   if( !nTracks ) return 0;
1041
1042   TObjArray *jets = new TObjArray(nTracks);
1043   TIter itrack(&tracks);
1044   while( nTracks ) {
1045     // 2- Start with the highest pT particle ...
1046     Float_t px,py,pz,pt; 
1047     AliAODTrack* track = (AliAODTrack*)itrack.Next();
1048     if( !track ) continue;
1049     px = track->Px();
1050     py = track->Py();
1051     pz = track->Pz();
1052     pt = track->Pt(); // Use the energy member to store Pt
1053     jets->AddLast( new TLorentzVector(px, py, pz, pt) );
1054     tracks.Remove( track );
1055     TLorentzVector* jet = (TLorentzVector*)jets->Last();
1056     jet->SetPtEtaPhiE( 1., jet->Eta(), jet->Phi(), pt );
1057     // 3- Go to the next highest pT particle not already included...
1058     AliAODTrack* track1;
1059     Double_t fPt = jet->E();
1060     while ( (track1  = (AliAODTrack*)(itrack.Next())) ) {
1061       Double_t tphi = track1->Phi(); // here Phi is from 0 <-> 2Pi
1062       if (tphi > TMath::Pi()) tphi -= 2. * TMath::Pi(); // convert to  -Pi <-> Pi
1063       Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-tphi);
1064       Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +
1065                                dphi*dphi );
1066       if( r < fConeRadius ) {
1067         fPt   = jet->E()+track1->Pt();  // Scalar sum of Pt
1068         // recalculating the centroid
1069         Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
1070         Double_t phi = jet->Phi()*jet->E()/fPt + tphi*track1->Pt()/fPt;
1071         jet->SetPtEtaPhiE( 1., eta, phi, fPt );
1072         tracks.Remove( track1 );
1073       }
1074     }
1075     
1076     tracks.Compress();
1077     nTracks = tracks.GetEntries();
1078     //   4- Continue until all particles are in a jet.
1079     itrack.Reset();
1080   } // end while nTracks
1081   
1082   // Convert to AODjets....
1083   Int_t njets = jets->GetEntriesFast();
1084   TObjArray* aodjets = new TObjArray(njets);
1085   aodjets->SetOwner(kTRUE);
1086   for(Int_t ijet=0; ijet<njets; ++ijet) {
1087     TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
1088     if (jet->E() < fChJetPtMin) continue;
1089     Float_t px, py,pz,en; // convert to 4-vector
1090     px = jet->E() * TMath::Cos(jet->Phi());  // Pt * cos(phi)
1091     py = jet->E() * TMath::Sin(jet->Phi());  // Pt * sin(phi)
1092     pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
1093     en = TMath::Sqrt(px * px + py * py + pz * pz);
1094
1095     aodjets->AddLast( new AliAODJet(px, py, pz, en) );
1096   }
1097   jets->Delete();
1098   delete jets;
1099   
1100   // Order jets according to their pT .
1101   QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
1102   
1103   // debug
1104   if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
1105   
1106   return aodjets;
1107 }
1108
1109 //____________________________________________________________________
1110 void  AliAnalysisTaskUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
1111 {
1112   // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
1113   
1114   static TObject *tmp;
1115   static int i;           // "static" to save stack space
1116   int j;
1117   
1118   while (last - first > 1) {
1119     i = first;
1120     j = last;
1121     for (;;) {
1122       while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
1123         ;
1124       while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
1125         ;
1126       if (i >= j)
1127         break;
1128       
1129       tmp  = a[i];
1130       a[i] = a[j];
1131       a[j] = tmp;
1132     }
1133     if (j == first) {
1134       ++first;
1135       continue;
1136     }
1137     tmp = a[first];
1138     a[first] = a[j];
1139     a[j] = tmp;
1140     if (j - first < last - (j + 1)) {
1141       QSortTracks(a, first, j);
1142       first = j + 1;   // QSortTracks(j + 1, last);
1143     } else {
1144       QSortTracks(a, j + 1, last);
1145       last = j;        // QSortTracks(first, j);
1146     }
1147   }
1148 }
1149
1150 //____________________________________________________________________
1151 void AliAnalysisTaskUE::SetRegionArea(TVector3 *jetVect)
1152 {
1153   Double_t fAreaCorrFactor=0.;
1154   Double_t deltaEta = 0.;
1155   if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.;
1156   else if (fRegionType==2){ 
1157     deltaEta = 0.9-TMath::Abs(jetVect[0].Eta());
1158     if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius;
1159     else{
1160       fAreaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
1161       (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
1162       fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-fAreaCorrFactor;
1163     }
1164   }else AliWarning("Unknown Rgion Type");
1165   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,fAreaCorrFactor));
1166 }
1167
1168 //____________________________________________________________________
1169 void  AliAnalysisTaskUE::CreateHistos()
1170 {
1171   fListOfHistos = new TList();
1172   
1173   
1174   fhNJets = new TH1F("fhNJets", "n Jet",  20, 0, 20);
1175   fhNJets->SetXTitle("# of jets");
1176   fhNJets->Sumw2();
1177   fListOfHistos->Add( fhNJets );                 // At(0) 
1178   
1179   fhEleadingPt  = new TH1F("hEleadingPt",   "Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1180   fhEleadingPt->SetXTitle("P_{T} (GeV/c)");
1181   fhEleadingPt->SetYTitle("dN/dP_{T}");
1182   fhEleadingPt->Sumw2();
1183   fListOfHistos->Add( fhEleadingPt );            // At(1)
1184   
1185   fhMinRegPtDist = new TH1F("hMinRegPtDist",   "P_{T} distribution in Min zone",  50,0.,20.);
1186   fhMinRegPtDist->SetXTitle("P_{T} (GeV/c)");
1187   fhMinRegPtDist->SetYTitle("dN/dP_{T}");
1188   fhMinRegPtDist->Sumw2();
1189   fListOfHistos->Add( fhMinRegPtDist );          // At(2)
1190   
1191   fhRegionMultMin = new TH1F("hRegionMultMin",      "N_{ch}^{90, min}",  21, -0.5,   20.5);
1192   fhRegionMultMin->SetXTitle("N_{ch tracks}");
1193   fhRegionMultMin->Sumw2();
1194   fListOfHistos->Add( fhRegionMultMin );         // At(3)            
1195   
1196   fhMinRegAvePt = new TH1F("hMinRegAvePt", "#LTp_{T}#GT",  50, 0.,   20.);
1197   fhMinRegAvePt->SetXTitle("P_{T} (GeV/c)");
1198   fhMinRegAvePt->Sumw2();
1199   fListOfHistos->Add( fhMinRegAvePt );           // At(4)
1200   
1201   fhMinRegSumPt = new TH1F("hMinRegSumPt", "#Sigma p_{T} ",  50, 0.,   20.);
1202   fhMinRegSumPt->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");  
1203   fhMinRegSumPt->SetXTitle("#Sigma p_{T} (GeV/c)");
1204   fhMinRegSumPt->Sumw2();
1205   fListOfHistos->Add( fhMinRegSumPt );           // At(5)
1206   
1207   fhMinRegMaxPtPart = new TH1F("hMinRegMaxPtPart", "max(p_{T})|_{event} ",  50, 0.,   20.);
1208   fhMinRegMaxPtPart->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");  
1209   fhMinRegMaxPtPart->SetXTitle("p_{T} (GeV/c)");
1210   fhMinRegMaxPtPart->Sumw2();
1211   fListOfHistos->Add( fhMinRegMaxPtPart );       // At(6)
1212   
1213   fhMinRegSumPtvsMult = new TH1F("hMinRegSumPtvsMult", "#Sigma p_{T} vs. Multiplicity ",  21, -0.5,   20.5);
1214   fhMinRegSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
1215   fhMinRegSumPtvsMult->SetXTitle("N_{charge}");
1216   fhMinRegSumPtvsMult->Sumw2();
1217   fListOfHistos->Add( fhMinRegSumPtvsMult );     // At(7);
1218   
1219   fhdNdEtaPhiDist  = new TH2F("hdNdEtaPhiDist", Form("Charge particle density |#eta|<%3.1f vs #Delta#phi", fTrackEtaCut),  
1220                                62, 0.,   2.*TMath::Pi(), fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1221   fhdNdEtaPhiDist->SetXTitle("#Delta#phi");
1222   fhdNdEtaPhiDist->SetYTitle("Leading Jet P_{T}");
1223   fhdNdEtaPhiDist->Sumw2();
1224   fListOfHistos->Add( fhdNdEtaPhiDist );        // At(8)
1225   
1226   // Can be use to get part pt distribution for differente Jet Pt bins
1227   fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading Jet P_{T}", fTrackEtaCut),
1228                                      200,0.,100., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1229   fhFullRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
1230   fhFullRegPartPtDistVsEt->SetXTitle("p_{T}");
1231   fhFullRegPartPtDistVsEt->Sumw2();
1232   fListOfHistos->Add( fhFullRegPartPtDistVsEt );  // At(9) 
1233   
1234   // Can be use to get part pt distribution for differente Jet Pt bins
1235   fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", Form( "dN/dP_{T} in tranvese regions |#eta|<%3.1f vs Leading Jet P_{T}", fTrackEtaCut),  
1236                                      200,0.,100., fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1237   fhTransRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
1238   fhTransRegPartPtDistVsEt->SetXTitle("p_{T}");
1239   fhTransRegPartPtDistVsEt->Sumw2();
1240   fListOfHistos->Add( fhTransRegPartPtDistVsEt );  // At(10)
1241   
1242   
1243   fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt",  "P_{T}^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1244   fhRegionSumPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
1245   fhRegionSumPtMaxVsEt->Sumw2();
1246   fListOfHistos->Add( fhRegionSumPtMaxVsEt );      // At(11)
1247   
1248   fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt",   "P_{T}^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1249   fhRegionSumPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
1250   fhRegionSumPtMinVsEt->Sumw2();
1251   fListOfHistos->Add( fhRegionSumPtMinVsEt );      // At(12)
1252   
1253   fhRegionMultMax = new TH1I("hRegionMultMax",      "N_{ch}^{90, max}",  21, -0.5,   20.5);
1254   fhRegionMultMax->SetXTitle("N_{ch tracks}");
1255   fhRegionMultMax->Sumw2();
1256   fListOfHistos->Add( fhRegionMultMax );           // At(13)
1257   
1258   fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt",  "N_{ch}^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1259   fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)");
1260   fhRegionMultMaxVsEt->Sumw2();
1261   fListOfHistos->Add( fhRegionMultMaxVsEt );       // At(14)
1262   
1263   fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt",  "N_{ch}^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1264   fhRegionMultMinVsEt->SetXTitle("E (GeV/c)");
1265   fhRegionMultMinVsEt->Sumw2();
1266   fListOfHistos->Add( fhRegionMultMinVsEt );      // At(15)
1267   
1268   fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1269   fhRegionAveSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
1270   fhRegionAveSumPtVsEt->Sumw2();
1271   fListOfHistos->Add( fhRegionAveSumPtVsEt );     // At(16)
1272   
1273   fhRegionDiffSumPtVsEt= new TH1F("hRegionPtDiffVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1274   fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
1275   fhRegionDiffSumPtVsEt->Sumw2();
1276   fListOfHistos->Add( fhRegionDiffSumPtVsEt );    // At(17)
1277   
1278   fhRegionAvePartPtMaxVsEt = new TH1F("hRegionAvePartPtMaxVsEt", "#LTp_{T}#GT^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1279   fhRegionAvePartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
1280   fhRegionAvePartPtMaxVsEt->Sumw2();
1281   fListOfHistos->Add( fhRegionAvePartPtMaxVsEt );  // At(18)
1282   
1283   fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1284   fhRegionAvePartPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
1285   fhRegionAvePartPtMinVsEt->Sumw2();
1286   fListOfHistos->Add( fhRegionAvePartPtMinVsEt );   // At(19)
1287   
1288   fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1289   fhRegionMaxPartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
1290   fhRegionMaxPartPtMaxVsEt->Sumw2();
1291   fListOfHistos->Add( fhRegionMaxPartPtMaxVsEt );    // At(20)
1292   /*
1293   fhRegForwardSumPtVsEt = new TH1F("hRegForwardSumPtVsEt", "Forward #sum{p_{T}} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1294   fhRegForwardSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
1295   fhRegForwardSumPtVsEt->Sumw2();
1296   fListOfHistos->Add( fhRegForwardSumPtVsEt );    // At(21)
1297   
1298   fhRegForwardMultVsEt = new TH1F("hRegForwardMultVsEt", "Forward #sum{N_{ch}} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1299   fhRegForwardMultVsEt->SetXTitle("P_{T} (GeV/c)");
1300   fhRegForwardMultVsEt->Sumw2();
1301   fListOfHistos->Add( fhRegForwardMultVsEt );    // At(22)
1302   
1303   fhRegBackwardSumPtVsEt = new TH1F("hRegBackwardSumPtVsEt", "Backward #sum{p_{T}} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1304   fhRegBackwardSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
1305   fhRegBackwardSumPtVsEt->Sumw2();
1306   fListOfHistos->Add( fhRegBackwardSumPtVsEt );    // At(23)
1307   
1308   fhRegBackwardMultVsEt = new TH1F("hRegBackwardMultVsEt", "Backward #sum{N_{ch}} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1309   fhRegBackwardMultVsEt->SetXTitle("P_{T} (GeV/c)");
1310   fhRegBackwardMultVsEt->Sumw2();
1311   fListOfHistos->Add( fhRegBackwardMultVsEt );    // At(24)
1312                                   */
1313   fhRegForwardMult = new TH2F("hRegForwardMult",      "N_{ch}^{forward}",  
1314                               fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5,   20.5);
1315   fhRegForwardMult->SetXTitle("N_{ch tracks}");
1316   fhRegForwardMult->Sumw2();
1317   fListOfHistos->Add( fhRegForwardMult );           // At(25)
1318   
1319   fhRegForwardSumPtvsMult = new TH2F("hRegForwardSumPtvsMult", "Forward #Sigma p_{T} vs. Multiplicity ",  
1320                                      fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5,   20.5);
1321   fhRegForwardSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
1322   fhRegForwardSumPtvsMult->SetXTitle("N_{charge}");
1323   fhRegForwardSumPtvsMult->Sumw2();
1324   fListOfHistos->Add( fhRegForwardSumPtvsMult );     // At(26);
1325   
1326   fhRegBackwardMult = new TH2F("hRegBackwardMult",      "N_{ch}^{backward}",  
1327                                fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5,   20.5);
1328   fhRegBackwardMult->SetXTitle("N_{ch tracks}");
1329   fhRegBackwardMult->Sumw2();
1330   fListOfHistos->Add( fhRegBackwardMult );           // At(27)
1331   
1332   fhRegBackwardSumPtvsMult = new TH2F("hRegBackwardSumPtvsMult", "Backward #Sigma p_{T} vs. Multiplicity ",  
1333                                       fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5,   20.5);
1334   fhRegBackwardSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
1335   fhRegBackwardSumPtvsMult->SetXTitle("N_{charge}");
1336   fhRegBackwardSumPtvsMult->Sumw2();
1337   fListOfHistos->Add( fhRegBackwardSumPtvsMult );     // At(28);
1338   
1339   fhRegForwardPartPtDistVsEt = new TH2F("hRegForwardPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading Jet P_{T}", fTrackEtaCut),
1340                                        200,0.,100., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1341   fhRegForwardPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
1342   fhRegForwardPartPtDistVsEt->SetXTitle("p_{T}");
1343   fhRegForwardPartPtDistVsEt->Sumw2();
1344   fListOfHistos->Add( fhRegForwardPartPtDistVsEt );  // At(29) 
1345   
1346   fhRegBackwardPartPtDistVsEt = new TH2F("hRegBackwardPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading Jet P_{T}", fTrackEtaCut),
1347                                          200,0.,100., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
1348   fhRegBackwardPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
1349   fhRegBackwardPartPtDistVsEt->SetXTitle("p_{T}");
1350   fhRegBackwardPartPtDistVsEt->Sumw2();
1351   fListOfHistos->Add( fhRegBackwardPartPtDistVsEt );  // At(30) 
1352   
1353   fhRegTransMult  = new TH2F("hRegTransMult",      "N_{ch}^{transv}",  
1354                              fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5,   20.5);
1355   fhRegTransMult->SetXTitle("N_{ch tracks}");
1356   fhRegTransMult->Sumw2();
1357   fListOfHistos->Add( fhRegTransMult );           // At(31)
1358   
1359   fhRegTransSumPtVsMult = new TH2F("hRegTransSumPtVsMult", "Transverse #Sigma p_{T} vs. Multiplicity ",  
1360                                    fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5,   20.5);
1361   fhRegTransSumPtVsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
1362   fhRegTransSumPtVsMult->SetXTitle("N_{charge}");
1363   fhRegTransSumPtVsMult->Sumw2();
1364   fListOfHistos->Add( fhRegTransSumPtVsMult );     // At(32);
1365   
1366   fhMinRegSumPtJetPtBin = new TH2F("hMinRegSumPtJetPtBin",      "Transverse Min Reg #Sigma p_{T} per jet Pt Bin",  
1367                                    fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 50, 0.,   20.);
1368   fhMinRegSumPtJetPtBin->SetXTitle("Leading Jet P_{T}");
1369   fhMinRegSumPtJetPtBin->Sumw2();
1370   fListOfHistos->Add( fhMinRegSumPtJetPtBin );           // At(33)
1371   
1372   fhMaxRegSumPtJetPtBin = new TH2F("hMaxRegSumPtJetPtBin",      "Transverse Max Reg #Sigma p_{T} per jet Pt Bin",  
1373                                    fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 50, 0.,   20.);
1374   fhMaxRegSumPtJetPtBin->SetXTitle("Leading Jet P_{T}");
1375   fhMaxRegSumPtJetPtBin->Sumw2();
1376   fListOfHistos->Add( fhMaxRegSumPtJetPtBin );           // At(34)
1377   
1378   fhVertexMult = new TH1F("hVertexMult",      "Multiplicity in Main Vertex", 81, -0.5 , 80.5);
1379   fhVertexMult->SetXTitle("Main Vertex Multiplicity");
1380   fhVertexMult->Sumw2();
1381   fListOfHistos->Add( fhVertexMult ); //At(35)
1382   
1383   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1); 
1384   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
1385   fh1Xsec->Sumw2();
1386   fListOfHistos->Add( fh1Xsec );            //At(36)
1387   
1388   fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
1389   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
1390   fh1Trials->Sumw2();
1391   fListOfHistos->Add( fh1Trials ); //At(37)
1392   
1393   fSettingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
1394   fSettingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
1395   fSettingsTree->Branch("fTrigger", &fTrigger,"TriggerFlag/I");
1396   fSettingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D");
1397   fSettingsTree->Branch("fJet1EtaCut", &fJet1EtaCut, "LeadJetEtaCut/D");
1398   fSettingsTree->Branch("fJet2DeltaPhiCut", &fJet2DeltaPhiCut, "DeltaPhi/D");
1399   fSettingsTree->Branch("fJet2RatioPtCut", &fJet2RatioPtCut, "Jet2Ratio/D");
1400   fSettingsTree->Branch("fJet3PtCut", &fJet3PtCut, "Jet3PtCut/D");
1401   fSettingsTree->Branch("fTrackPtCut", &fTrackPtCut, "TrackPtCut/D");
1402   fSettingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
1403   fSettingsTree->Branch("fAnaType", &fAnaType, "Ana/I");        
1404   fSettingsTree->Branch("fRegionType", &fRegionType,"Reg/I");
1405   fSettingsTree->Branch("fOrdering", &fOrdering,"OrderMeth/I");
1406   fSettingsTree->Branch("fUseChPartJet", &fUseChPartJet,"UseChPart/O");
1407   fSettingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
1408   fSettingsTree->Branch("fUseSingleCharge", &fUseSingleCharge,"UseSingleCh/O");
1409   fSettingsTree->Branch("fUsePositiveCharge", &fUsePositiveCharge,"UsePositiveCh/O");
1410   fSettingsTree->Fill();
1411
1412   
1413   fListOfHistos->Add( fSettingsTree );    // At(38)
1414   
1415   /*   
1416    // For debug region selection
1417    fhValidRegion = new TH2F("hValidRegion", "dN/d#eta/d#phi",      
1418    20, -2.,2., 62, -TMath::Pi(),   TMath::Pi());
1419    fhValidRegion->SetYTitle("#Delta#phi");
1420    fhValidRegion->Sumw2();
1421    fListOfHistos->Add( fhValidRegion );  // At(15)
1422    */
1423 }
1424
1425
1426
1427 //____________________________________________________________________
1428 void  AliAnalysisTaskUE::Terminate(Option_t */*option*/)
1429 {
1430   // Terminate analysis
1431   //
1432   
1433   //Save Analysis Settings
1434   //  WriteSettings();
1435   
1436   // Normalize histos to region area TODO: 
1437   // Normalization done at Analysis, taking into account 
1438   // area variations on per-event basis (cone case)
1439    
1440    //HIGH WARNING!!!!!: DO NOT SCALE ANY OF THE ORIGINAL HISTOGRAMS
1441    //MAKE A COPY, DRAW IT, And later sacale that copy. CAF Issue!!!!!
1442   
1443   // Draw some Test plot to the screen
1444   if (!gROOT->IsBatch()) {
1445     
1446     // Update pointers reading them from the output slot
1447     fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
1448     if( !fListOfHistos ) {
1449       AliError("Histogram List is not available");
1450       return;
1451     }
1452     fhNJets              = (TH1F*)fListOfHistos->At(0);
1453     fhEleadingPt         = (TH1F*)fListOfHistos->At(1);
1454     fhdNdEtaPhiDist      = (TH2F*)fListOfHistos->At(8);
1455     fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11);
1456     fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12);
1457     fhRegionMultMaxVsEt  = (TH1F*)fListOfHistos->At(14);
1458     fhRegionMultMinVsEt  = (TH1F*)fListOfHistos->At(15);
1459     fhRegionAveSumPtVsEt = (TH1F*)fListOfHistos->At(16);
1460     //fhRegForwardSumPtVsEt = (TH1F*)fListOfHistos->At(21);
1461     //fhRegBackwardSumPtVsEt = (TH1F*)fListOfHistos->At(23);
1462     
1463     //fhValidRegion  = (TH2F*)fListOfHistos->At(21);
1464     
1465     // Canvas
1466     TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1100,700);
1467     c1->Divide(2,2);
1468     c1->cd(1);
1469     TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
1470     h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1);
1471     //h1r->Scale( areafactor );
1472     h1r->SetMarkerStyle(20);
1473     h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1474     h1r->SetYTitle("P_{T}^{90, max}");
1475     h1r->DrawCopy("p");
1476     
1477     c1->cd(2);
1478     
1479     TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
1480     h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1);
1481     //h2r->Scale( areafactor );
1482     h2r->SetMarkerStyle(20);
1483     h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1484     h2r->SetYTitle("P_{T}^{90, min}");
1485     h2r->DrawCopy("p");
1486     
1487     c1->cd(3);
1488     TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
1489     //TH1F *h41r = new TH1F("hRegForwvsDiffPt" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
1490     //TH1F *h42r = new TH1F("hRegBackvsDiffPt" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
1491     //h41r->Divide(fhRegForwardSumPtVsEt,fhEleadingPt,1,1);
1492     //h42r->Divide(fhRegBackwardSumPtVsEt,fhEleadingPt,1,1);
1493     h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1);
1494     //h4r->Scale(2.); // make average
1495     //h4r->Scale( areafactor );
1496     h4r->SetYTitle("#DeltaP_{T}^{90}");
1497     h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1498     h4r->SetMarkerStyle(20);
1499     h4r->DrawCopy("p");
1500     /*h41r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1501     h41r->SetMarkerStyle(22);
1502     h41r->DrawCopy("p same");
1503     h42r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1504     h42r->SetMarkerStyle(23);
1505     h42r->SetMarkerColor(kRed);
1506     h42r->DrawCopy("p same");
1507     */
1508     c1->cd(4);
1509     TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading",   "",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1510     TH1F *h6r = new TH1F("hRegionMultMinVsEtleading",   "",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
1511     
1512     h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1);
1513     h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1);
1514     //h5r->Scale( areafactor );
1515     //h6r->Scale( areafactor );
1516     h5r->SetYTitle("N_{Tracks}^{90}");
1517     h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1518     h5r->SetMarkerStyle(20);
1519     h5r->DrawCopy("p");
1520     h6r->SetMarkerStyle(21);
1521     h6r->SetMarkerColor(2);
1522     h6r->SetYTitle("N_{Tracks}^{90}");
1523     h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
1524     h6r->DrawCopy("p same");
1525     c1->Update();
1526     
1527     //Get Normalization
1528     fh1Xsec                = (TProfile*)fListOfHistos->At(21);
1529     fh1Trials              = (TH1F*)fListOfHistos->At(22);
1530     
1531     Double_t xsec = fh1Xsec->GetBinContent(1);
1532     Double_t ntrials = fh1Trials->GetBinContent(1);
1533     Double_t normFactor = xsec/ntrials;
1534     if(fDebug > 1)Printf("xSec %f nTrials %f Norm %f \n",xsec,ntrials,normFactor);
1535     
1536     
1537     
1538     TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800);
1539     TH1 * copy = 0;
1540     c2->Divide(2,2);
1541     c2->cd(1);
1542     fhEleadingPt->SetMarkerStyle(20);
1543     fhEleadingPt->SetMarkerColor(2);
1544     //if( normFactor > 0.) fhEleadingPt->Scale(normFactor);
1545     //fhEleadingPt->Draw("p");
1546     copy = fhEleadingPt->DrawCopy("p");
1547     if( normFactor > 0.) copy->Scale(normFactor);
1548     gPad->SetLogy();
1549     
1550     c2->cd(2);
1551     Int_t xbin1 = fhdNdEtaPhiDist->GetYaxis()->FindFixBin(fMinJetPtInHist);
1552     Int_t xbin2 = fhdNdEtaPhiDist->GetYaxis()->FindFixBin(fMaxJetPtInHist);
1553     TH1D* dNdEtaPhiDistAllJets = fhdNdEtaPhiDist->ProjectionX("dNdEtaPhiDistAllJets",xbin1,xbin2);
1554     dNdEtaPhiDistAllJets->SetMarkerStyle(20);
1555     dNdEtaPhiDistAllJets->SetMarkerColor(2);
1556     dNdEtaPhiDistAllJets->DrawCopy("p");
1557     gPad->SetLogy();
1558     
1559     c2->cd(3);      
1560     fhNJets->DrawCopy();
1561     
1562     //c2->cd(4);      
1563     //fhValidRegion->DrawCopy("p");
1564     
1565     //fhTransRegPartPtDist  = (TH1F*)fListOfHistos->At(2);
1566     fhRegionMultMin          = (TH1F*)fListOfHistos->At(3);
1567     fhMinRegAvePt     = (TH1F*)fListOfHistos->At(4);
1568     fhMinRegSumPt     = (TH1F*)fListOfHistos->At(5);   
1569     //fhMinRegMaxPtPart     = (TH1F*)fListOfHistos->At(6);
1570     fhMinRegSumPtvsMult   = (TH1F*)fListOfHistos->At(7);
1571     
1572     
1573     // Canvas
1574     TCanvas* c3 = new TCanvas("c3"," pT dist",160,160,1200,800);
1575     c3->Divide(2,2);
1576     c3->cd(1);
1577     /*fhTransRegPartPtDist->SetMarkerStyle(20);
1578      fhTransRegPartPtDist->SetMarkerColor(2); 
1579      fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries());
1580      fhTransRegPartPtDist->DrawCopy("p");
1581      gPad->SetLogy();
1582      */
1583     c3->cd(2); 
1584     fhMinRegSumPt->SetMarkerStyle(20);
1585     fhMinRegSumPt->SetMarkerColor(2);  
1586     //fhMinRegSumPt->Scale(areafactor);
1587     fhMinRegSumPt->DrawCopy("p");
1588     gPad->SetLogy();
1589     
1590     c3->cd(3);
1591     fhMinRegAvePt->SetMarkerStyle(20);
1592     fhMinRegAvePt->SetMarkerColor(2);  
1593     //fhMinRegAvePt->Scale(areafactor);
1594     fhMinRegAvePt->DrawCopy("p");
1595     gPad->SetLogy();
1596     
1597     c3->cd(4);
1598     
1599     TH1F *h7r = new TH1F("hRegionMultMinVsMult",   "",  21, -0.5,   20.5);
1600     
1601     h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1);
1602     
1603     h7r->SetMarkerStyle(20);
1604     h7r->SetMarkerColor(2);   
1605     h7r->DrawCopy("p");
1606     
1607     c3->Update();
1608     
1609     
1610     /*    c2->cd(3);      
1611      fhValidRegion->SetMarkerStyle(20);
1612      fhValidRegion->SetMarkerColor(2);
1613      fhValidRegion->DrawCopy("p");
1614      */    
1615     c2->Update();
1616   } else {
1617     AliInfo(" Batch mode, not histograms will be shown...");
1618   }
1619   
1620   if( fDebug > 1 ) 
1621     AliInfo("End analysis");
1622   
1623 }
1624
1625 void  AliAnalysisTaskUE::WriteSettings()
1626
1627   if (fDebug>5){
1628     AliInfo(" All Analysis Settings in Saved Tree");
1629     fSettingsTree->Scan();
1630   }
1631 }