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