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