Fixed compiler warnings
[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 "AliAODHandler.h"
41 #include "AliAODInputHandler.h"
42 #include "AliAODJet.h"
43 #include "AliAODMCParticle.h"
44 #include "AliAODTrack.h"
45 #include "AliKFVertex.h"
46 #include "AliMCEvent.h"
47 #include "AliMCEventHandler.h"
48 #include "AliStack.h"
49
50 #include "AliAnalysisHelperJetTasks.h"
51 #include "AliGenPythiaEventHeader.h"
52 #include "AliInputEventHandler.h"
53 #include "AliLog.h"
54 #include "AliStack.h"
55
56 ////////////////////////////////////////////////
57 //--------------------------------------------- 
58 // Class for transverse regions analysis
59 //---------------------------------------------
60 ////////////////////////////////////////////////
61
62
63 using namespace std;
64
65 ClassImp(AliAnalyseUE)
66
67 //-------------------------------------------------------------------
68 AliAnalyseUE::AliAnalyseUE() :
69   TObject(),
70   //fTaskUE(0),
71   fkAOD(0x0),            
72   fDebug(0),
73   fSimulateChJetPt(kFALSE),
74   fAnaType(1),         
75   fAreaReg(1.5393), // Pi*0.7*0.7
76   fConeRadius(0.7),
77   fFilterBit(0xFF),
78   fRegionType(1),
79   fUseChargeHadrons(kFALSE),
80   fUseChPartJet(kFALSE),
81   fUsePositiveCharge(kTRUE),
82   fUseSingleCharge(kFALSE),
83   fOrdering(1),
84   fJet1EtaCut(0.2),
85   fJet2DeltaPhiCut(2.616),    // 150 degrees
86   fJet2RatioPtCut(0.8),
87   fJet3PtCut(15.),
88   fTrackEtaCut(0.9),
89   fTrackPtCut(0.),
90   fHistos(0x0),
91   fSumPtRegionPosit(0.),
92   fSumPtRegionNegat(0.),
93   fSumPtRegionForward(0.),
94   fSumPtRegionBackward(0.),
95   fMaxPartPtRegion(0.),
96   fNTrackRegionPosit(0),
97   fNTrackRegionNegat(0),
98   fNTrackRegionForward(0),
99   fNTrackRegionBackward(0),
100   fSettingsTree(0x0)
101
102 {
103   // constructor
104 }
105
106
107 //-------------------------------------------------------------------
108 AliAnalyseUE::AliAnalyseUE(const AliAnalyseUE & original) :
109   TObject(original),
110   //fTaskUE(original.fTaskUE),
111   fkAOD(original.fkAOD),            
112   fDebug(original.fDebug),
113   fSimulateChJetPt(original.fSimulateChJetPt),
114   fAnaType(original.fAnaType),
115   fAreaReg(original.fAreaReg),
116   fConeRadius(original.fConeRadius),
117   fFilterBit(original.fFilterBit),
118   fRegionType(original.fRegionType),
119   fUseChargeHadrons(original.fUseChargeHadrons),
120   fUseChPartJet(original.fUseChPartJet),
121   fUsePositiveCharge(original.fUsePositiveCharge),
122   fUseSingleCharge(original.fUseSingleCharge),
123   fOrdering(original.fOrdering),
124   fJet1EtaCut(original.fJet1EtaCut),
125   fJet2DeltaPhiCut(original.fJet2DeltaPhiCut),
126   fJet2RatioPtCut(original.fJet2RatioPtCut),
127   fJet3PtCut(original.fJet3PtCut),
128   fTrackEtaCut(original.fTrackEtaCut),
129   fTrackPtCut(original.fTrackPtCut),
130   fHistos(original.fHistos),
131   fSumPtRegionPosit(original.fSumPtRegionPosit),
132   fSumPtRegionNegat(original.fSumPtRegionNegat),
133   fSumPtRegionForward(original.fSumPtRegionForward),
134   fSumPtRegionBackward(original.fSumPtRegionBackward),
135   fMaxPartPtRegion(original.fMaxPartPtRegion),
136   fNTrackRegionPosit(original.fNTrackRegionPosit),
137   fNTrackRegionNegat(original.fNTrackRegionNegat),
138   fNTrackRegionForward(original.fNTrackRegionForward),
139   fNTrackRegionBackward(original.fNTrackRegionBackward),
140   fSettingsTree(original.fSettingsTree)
141 {
142   //copy constructor    
143 }
144
145 //-------------------------------------------------------------------
146 AliAnalyseUE & AliAnalyseUE::operator = (const AliAnalyseUE & /*source*/)
147 {
148   // assignment operator
149   return *this;
150 }
151
152
153 //-------------------------------------------------------------------
154 AliAnalyseUE::~AliAnalyseUE(){
155
156   //clear memory
157   delete[] fkAOD;
158   fkAOD = NULL;
159
160   
161   
162 }
163
164
165 //-------------------------------------------------------------------
166 void AliAnalyseUE::AnalyseMC(TVector3 *jetVect,AliMCEvent *mcEvent, AliGenPythiaEventHeader *pythiaGenHeader,Int_t conePosition, Bool_t useAliStack, Bool_t constrainDistance, Double_t minDistance){
167
168   // Execute the analysis in case of MC input
169   fSumPtRegionPosit = 0.;
170   fSumPtRegionNegat = 0.;
171   fSumPtRegionForward = 0.;
172   fSumPtRegionBackward = 0.;
173   fMaxPartPtRegion = 0.;
174   fNTrackRegionPosit = 0;
175   fNTrackRegionNegat = 0;
176   fNTrackRegionForward = 0;
177   fNTrackRegionBackward = 0;
178
179   static Double_t const  kPI     = TMath::Pi();
180   static Double_t const  k270rad = 270.*kPI/180.;
181
182   //Get Jets from MC header
183   Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
184   AliAODJet pythiaGenJets[4];
185   TVector3 jetVectnew[4];
186   Int_t iCount = 0;
187   for(int ip = 0;ip < nPythiaGenJets;++ip){
188         if (iCount>3) break;
189         Float_t p[4];
190         pythiaGenHeader->TriggerJet(ip,p);
191         TVector3 tempVect(p[0],p[1],p[2]);
192         if ( TMath::Abs(tempVect.Eta())>fJet1EtaCut ) continue;
193         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
194         jetVectnew[iCount].SetXYZ(pythiaGenJets[iCount].Px(), pythiaGenJets[iCount].Py(), pythiaGenJets[iCount].Pz());
195         iCount++;
196         }
197
198   if (!iCount) return;// no jet in eta acceptance
199     
200   //Search the index of the nearest MC jet to the leading jet reconstructed from the input data
201   Int_t index = 0;
202   if (constrainDistance){
203         Float_t deltaR = 0.;
204         Float_t dRTemp = 0.;
205         for (Int_t i=0; i<iCount; i++){
206                 if (!i) {
207                         dRTemp = jetVectnew[i].DeltaR(jetVect[0]);
208                         index = i;
209                         }
210
211                 deltaR = jetVectnew[i].DeltaR(jetVect[0]);
212                 if (deltaR < dRTemp){
213                         index = i;
214                         dRTemp = deltaR;
215                         }
216                 }
217    
218         if (jetVectnew[index].DeltaR(jetVect[0]) > minDistance) return;
219         }
220
221   //Let's add some taste to jet and simulate pt of charged alone 
222   //eta and phi are kept as original
223   //Play a Normal Distribution
224   Float_t random = 1.;  
225   if (fSimulateChJetPt){
226         while(1){
227                 random = gRandom->Gaus(0.6,0.25);
228                 if (random > 0. && random < 1. && 
229                 (random * jetVectnew[index].Pt()>6.)) break;
230                 }
231         }
232     
233   //Set new Pt & Fill histogram accordingly
234   Double_t maxPtJet1 = random * jetVectnew[index].Pt();  
235     
236   fHistos->FillHistogram("hEleadingPt", maxPtJet1 );    
237     
238   if (useAliStack){//Try Stack Information to perform UE analysis
239     
240         AliStack* mcStack = mcEvent->Stack();//Load Stack
241         Int_t nTracksMC = mcStack->GetNtrack();
242         for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
243                 //Cuts
244                 if(!(mcStack->IsPhysicalPrimary(iTracks))) continue;
245         
246                 TParticle* mctrk = mcStack->Particle(iTracks);
247         
248                 Double_t charge = mctrk->GetPDG()->Charge();
249                 Double_t pT = mctrk->Pt();
250                 Double_t eta = mctrk->Eta();
251                 Int_t pdgCode = mctrk->GetPdgCode();
252
253                 if (!TrackMCSelected(charge, pT, eta, pdgCode))continue;
254
255                 TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
256                 Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
257                 if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
258                 fHistos->FillHistogram("hdNdEtaPhiDist",deltaPhi, maxPtJet1 ); 
259         
260                 fHistos->FillHistogram("hFullRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 ); 
261         
262                 //We are not interested on stack organization but don't loose track of info
263         
264                 TVector3 tempVector =  jetVectnew[0];
265                 jetVectnew[0] = jetVectnew[index];
266                 jetVectnew[index] = tempVector;
267         
268                 Int_t region = IsTrackInsideRegion( jetVectnew, &partVect, conePosition );  
269         
270                 if (region == 1) {
271                         if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
272                         fSumPtRegionPosit += mctrk->Pt();
273                         fNTrackRegionPosit++;
274                         fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
275                         }
276                 if (region == -1) {
277                         if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
278                         fSumPtRegionNegat += mctrk->Pt();
279                         fNTrackRegionNegat++;
280                         fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
281                         }
282                 if (region == 2){ //forward
283                         fSumPtRegionForward += mctrk->Pt();
284                         fNTrackRegionForward++;
285                         fHistos->FillHistogram("hRegForwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
286                         }
287                 if (region == -2){ //backward
288                         fSumPtRegionBackward += mctrk->Pt();
289                         fNTrackRegionBackward++;
290                         fHistos->FillHistogram("hRegBackwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
291                         }
292                 } //end loop on stack particles     
293     }else{//Try mc Particle
294
295       TClonesArray* farray = (TClonesArray*)fkAOD->FindListObject("mcparticles");
296        
297       Int_t ntrks = farray->GetEntries();
298       if (fDebug>1) AliInfo(Form("In UE MC analysis tracks %d \n",ntrks));
299       for (Int_t i =0 ; i < ntrks; i++){   
300         AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
301         //Cuts
302         if (!(mctrk->IsPhysicalPrimary())) continue;
303         //if (!(mctrk->IsPrimary())) continue;
304         
305         Double_t charge = mctrk->Charge(); 
306         Double_t pT = mctrk->Pt();
307         Double_t eta = mctrk->Eta();
308         Int_t pdgCode = mctrk->GetPdgCode();
309
310         if (!TrackMCSelected(charge, pT, eta, pdgCode))continue;
311         
312         TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
313
314         Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
315         if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
316         fHistos->FillHistogram("hdNdEtaPhiDist", deltaPhi, maxPtJet1 );
317
318         fHistos->FillHistogram("hFullRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
319         
320         //We are not interested on stack organization but don't loose track of info
321         TVector3 tempVector =  jetVectnew[0];
322         jetVectnew[0] = jetVectnew[index];
323         jetVectnew[index] = tempVector;
324         
325         Int_t region = IsTrackInsideRegion( jetVectnew, &partVect, conePosition );  
326         
327         if (region == 1) { //right
328           if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
329           fSumPtRegionPosit += mctrk->Pt();
330           fNTrackRegionPosit++;
331           fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
332         }
333         if (region == -1) { //left
334           if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
335           fSumPtRegionNegat += mctrk->Pt();
336           fNTrackRegionNegat++;
337           fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
338         }
339         if (region == 2){ //forward
340           fSumPtRegionForward += mctrk->Pt();
341           fNTrackRegionForward++;
342           fHistos->FillHistogram("hRegForwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
343         }
344         if (region == -2){ //backward
345           fSumPtRegionBackward += mctrk->Pt();
346           fNTrackRegionBackward++;
347           fHistos->FillHistogram("hRegBackwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
348         }
349         
350       }//end loop AliAODMCParticle tracks
351    }  
352 }
353
354
355
356 //-------------------------------------------------------------------
357 Bool_t AliAnalyseUE::AnaTypeSelection(TVector3 *jetVect ){
358
359   // Cut events by jets topology
360   // anaType:
361   //     1 = inclusive,
362   //         - Jet1 |eta| < jet1EtaCut
363   //     2 = back to back inclusive
364   //         - fulfill case 1
365   //         - |Jet1.Phi - Jet2.Phi| > jet2DeltaPhiCut
366   //         - Jet2.Pt/Jet1Pt > jet2RatioPtCut
367   //     3 = back to back exclusive
368   //         - fulfill case 2
369   //         - Jet3.Pt < jet3PtCut
370
371   Double_t eta=jetVect[0].Eta();
372   if( TMath::Abs(eta) > fJet1EtaCut) {
373         if( fDebug > 1 ) AliInfo("\n   Skipping Event...Jet1 |eta| > fJet1EtaCut");
374         return kFALSE;
375         }
376   // back to back inclusive
377   if( fAnaType > 1 && fAnaType < 4 && jetVect[1].Pt() < 0. ) {
378         if( fDebug > 1 ) AliInfo("\n   Skipping Event... no second Jet found");
379         return kFALSE;
380         }
381   if( fAnaType > 1 && fAnaType < 4 && jetVect[1].Pt() > 0. ) {
382         if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
383         jetVect[1].Pt()/jetVect[0].Pt() < fJet2RatioPtCut ) {
384                 if( fDebug > 1 ) AliInfo("\n   Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
385                 return kFALSE;
386                 }
387         }
388   // back to back exclusive
389   if( fAnaType > 2 && fAnaType < 4 && jetVect[2].Pt() > 0. ) {
390         if( jetVect[2].Pt() > fJet3PtCut ) {
391                 if( fDebug > 1 ) AliInfo("\n   Skipping Event... Jet3.Pt > fJet3PtCut ");
392                 return kFALSE;
393                 }
394         }
395   return kTRUE;  
396 }
397
398
399 //-------------------------------------------------------------------
400 void AliAnalyseUE::FillRegions(Bool_t isNorm2Area,  TVector3 *jetVect){
401   
402   // Fill the different topological regions
403   Double_t maxPtJet1 = jetVect[0].Pt();
404   static Double_t const  kPI     = TMath::Pi();
405   static Double_t const  k120rad = 120.*kPI/180.;
406   Double_t const kMyTolerance = 0.0000001;
407
408   //Area for Normalization 
409   // Forward and backward
410   Double_t normArea = 1.;
411   // Transverse
412   if (isNorm2Area) {
413         SetRegionArea(jetVect);
414         normArea =  2.*fTrackEtaCut*k120rad ;
415         } else fAreaReg = 1.;
416   
417   Double_t avePosRegion = (fNTrackRegionPosit) ? fSumPtRegionPosit/fNTrackRegionPosit : 0.;
418   Double_t aveNegRegion = (fNTrackRegionNegat) ? fSumPtRegionNegat/fNTrackRegionNegat : 0.;
419   if( avePosRegion > aveNegRegion ) {
420      FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg );
421   } else {
422      FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg );
423   }
424   
425   //How quantities will be sorted before Fill Min and Max Histogram
426   //  1=Plots will be CDF-like
427   //  2=Plots will be Marchesini-like
428   //  3=Minimum zone is selected as the one having lowest pt per track 
429   if( fOrdering == 1 ) {
430     if( fSumPtRegionPosit > fSumPtRegionNegat ) {
431       FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
432     } else {
433       FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
434     }
435     if (fNTrackRegionPosit > fNTrackRegionNegat ) {
436       FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
437     } else {
438       FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
439     }
440   } else if( fOrdering == 2 ) {
441     if (fSumPtRegionPosit > fSumPtRegionNegat) {
442       FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
443       FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
444     } else {
445       FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
446       FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
447     }
448   } else if( fOrdering == 3 ){
449      if (avePosRegion > aveNegRegion) {
450         FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
451         FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
452      }else{
453         FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
454         FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
455      }
456   }
457   fHistos->FillHistogram("hRegionMaxPartPtMaxVsEt",maxPtJet1, fMaxPartPtRegion);  
458   
459   // Compute pedestal like magnitudes
460   fHistos->FillHistogram("hRegionDiffSumPtVsEt",maxPtJet1, (TMath::Abs(fSumPtRegionPosit-fSumPtRegionNegat)/(2.0*fAreaReg))+kMyTolerance);
461   fHistos->FillHistogram("hRegionAveSumPtVsEt", maxPtJet1, (fSumPtRegionPosit+fSumPtRegionNegat)/(2.0*fAreaReg));
462
463   // Transverse as a whole
464   fHistos->FillHistogram("hRegTransMult", maxPtJet1, fNTrackRegionPosit + fNTrackRegionNegat, (fNTrackRegionPosit + fNTrackRegionNegat)/(2.0*fAreaReg));
465  fHistos->FillHistogram("hRegTransSumPtVsMult",maxPtJet1, fNTrackRegionPosit + fNTrackRegionNegat , (fSumPtRegionNegat + fSumPtRegionPosit)/(2.0 *fAreaReg));
466
467   // Fill Histograms for Forward and away side w.r.t. leading jet direction
468   // Pt dependence
469   //fHistos->FillHistogram("hRegForwardSumPtVsEt",maxPtJet1, fSumPtRegionForward/normArea );
470   //fHistos->FillHistogram("hRegForwardMultVsEt",maxPtJet1, fNTrackRegionForward/normArea );
471   //fHistos->FillHistogram("hRegBackwardSumPtVsEt",maxPtJet1, fSumPtRegionBackward/normArea );
472   //fHistos->FillHistogram("hRegBackwardMultVsEt",maxPtJet1, fNTrackRegionBackward/normArea);
473   
474   // Multiplicity dependence
475   fHistos->FillHistogram("hRegForwardMult", maxPtJet1, fNTrackRegionForward, fNTrackRegionForward/normArea);
476   fHistos->FillHistogram("hRegForwardSumPtvsMult", maxPtJet1, fNTrackRegionForward,fSumPtRegionForward/normArea);
477   fHistos->FillHistogram("hRegBackwardMult", maxPtJet1, fNTrackRegionBackward, fNTrackRegionBackward/normArea );
478   fHistos->FillHistogram("hRegBackwardSumPtvsMult", maxPtJet1, fNTrackRegionBackward,fSumPtRegionBackward/normArea);
479 }
480
481
482 //-------------------------------------------------------------------
483 void AliAnalyseUE::FindMaxMinRegions(TVector3 *jetVect, Int_t conePosition){
484   
485   // Identify the different topological zones
486   fSumPtRegionPosit = 0.;
487   fSumPtRegionNegat = 0.;
488   fSumPtRegionForward = 0.;
489   fSumPtRegionBackward = 0.;
490   fMaxPartPtRegion = 0.;
491   fNTrackRegionPosit = 0;
492   fNTrackRegionNegat = 0;
493   fNTrackRegionForward = 0;
494   fNTrackRegionBackward = 0;
495   static Double_t const  kPI     = TMath::Pi();
496   static Double_t const  kTWOPI  = 2.*kPI;
497   static Double_t const  k270rad = 270.*kPI/180.;
498   Double_t const kMyTolerance = 0.0000001;
499
500     Int_t nTracks = fkAOD->GetNTracks();
501     if (fDebug > 1) AliInfo(Form(" ==== AOD tracks = %d \n ",nTracks));
502     
503     for (Int_t ipart=0; ipart<nTracks; ++ipart) {
504       
505     AliAODTrack* part = fkAOD->GetTrack( ipart );
506     if (fDebug > 1) AliInfo(Form(" ==== AOD track = %d pt = %f charge = %d \n ",ipart,part->Pt(),part->Charge()));
507     // track selection
508     if (! TrackSelected(part)) continue;
509       
510       
511     TVector3 partVect(part->Px(), part->Py(), part->Pz());
512     Bool_t isFlagPart = kTRUE;
513     Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
514     if( deltaPhi > kTWOPI )  deltaPhi-= kTWOPI;
515     if (fAnaType != 4 ) fHistos->FillHistogram("hdNdEtaPhiDist",deltaPhi, jetVect[0].Pt());
516     else if (TMath::Abs(deltaPhi-k270rad) >= kMyTolerance && TMath::Abs(jetVect[0].Eta()-partVect.Eta()) >= kMyTolerance){
517         fHistos->FillHistogram("hdNdEtaPhiDist",deltaPhi, jetVect[0].Pt());
518         isFlagPart = kFALSE;
519         }
520       
521     fHistos->FillHistogram("hFullRegPartPtDistVsEt", part->Pt(), jetVect[0].Pt());  
522     
523     Int_t region = IsTrackInsideRegion( jetVect, &partVect, conePosition );  
524     if (region == 1) {
525         if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
526         fSumPtRegionPosit += part->Pt();
527         fNTrackRegionPosit++;
528         fHistos->FillHistogram("hTransRegPartPtDistVsEt",part->Pt(), jetVect[0].Pt());
529         }
530     if (region == -1) {
531         if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
532         fSumPtRegionNegat += part->Pt();
533         fNTrackRegionNegat++;
534         fHistos->FillHistogram("hTransRegPartPtDistVsEt",part->Pt(), jetVect[0].Pt());
535         }
536     if (region == 2){ //forward
537         fSumPtRegionForward += part->Pt();
538         fNTrackRegionForward++;
539         fHistos->FillHistogram("hRegForwardPartPtDistVsEt",part->Pt(), jetVect[0].Pt());
540         }
541     if (region == -2){ //backward
542         fSumPtRegionBackward += part->Pt();
543         fNTrackRegionBackward++;
544         fHistos->FillHistogram("hRegBackwardPartPtDistVsEt",part->Pt(), jetVect[0].Pt());
545         }
546     }//end loop AOD tracks
547
548 }
549
550 //-------------------------------------------------------------------
551 TVector3 AliAnalyseUE::GetOrderedClusters(TString aodBranch, Bool_t chargedJets, Double_t chJetPtMin){ 
552
553   // jets from AOD, on-the-fly or leading particle
554   Double_t maxPtJet1 = 0.; 
555   Int_t    index1 = -1;
556   Double_t maxPtJet2 = 0.;   // jet 2 need for back to back inclusive
557   Int_t    index2 = -1;
558   Double_t maxPtJet3 = 0.;   // jet 3 need for back to back exclusive
559   Int_t    index3 = -1;
560   TVector3 jetVect[3];
561   
562   jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
563   jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
564   jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
565   
566   Int_t nJets = 0;
567   //TClonesArray* fArrayJets;
568   TObjArray* arrayJets;
569   // 1) JETS FROM AOD BRANCH (standard, non-standard or delta)
570   if (!chargedJets && fAnaType != 4 ) { 
571         AliInfo(" ==== Read AODs  !");
572         AliInfo(Form(" ====  Reading Branch: %s  ", aodBranch.Data()));
573         arrayJets = (TObjArray*)fkAOD->GetList()->FindObject(aodBranch.Data());
574         if (!arrayJets){
575                AliFatal(" No jet-array! ");
576                return *jetVect;
577                }
578
579         // Find Leading Jets 1,2,3 
580         // (could be skipped if Jets are sort by Pt...)
581         nJets=arrayJets->GetEntries();
582         for( Int_t i=0; i<nJets; ++i ) {
583                 AliAODJet* jet = (AliAODJet*)arrayJets->At(i);
584                 Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!!
585  
586                 if( jetPt > maxPtJet1 ) {
587                         maxPtJet3 = maxPtJet2; index3 = index2;
588                         maxPtJet2 = maxPtJet1; index2 = index1;
589                         maxPtJet1 = jetPt; index1 = i;
590                         } else if( jetPt > maxPtJet2 ) {
591                         maxPtJet3 = maxPtJet2; index3 = index2;
592                         maxPtJet2 = jetPt; index2 = i;
593                         } else if( jetPt > maxPtJet3 ) {
594                         maxPtJet3 = jetPt; index3 = i;
595                         }
596                 }
597
598         if( index1 != -1 ) {
599                 AliAODJet *jet =(AliAODJet*) arrayJets->At(index1);
600                 if(jet)jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
601                 }
602         if( index2 != -1 ) {
603                 AliAODJet* jet= (AliAODJet*) arrayJets->At(index2);
604                 if(jet)jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
605                 }
606         if( index3 != -1 ) {
607                 AliAODJet* jet = (AliAODJet*) arrayJets->At(index3);
608                 if(jet)jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
609                 }
610     
611   }
612
613
614   // 2) ON-THE-FLY CDF ALGORITHM
615   if (chargedJets){ 
616         arrayJets = FindChargedParticleJets(chJetPtMin);
617         if( arrayJets ) {
618                 nJets = arrayJets->GetEntriesFast();
619                 if( nJets > 0 ) {
620                         index1 = 0;
621                         AliAODJet* jet = (AliAODJet*)arrayJets->At(0);
622                         maxPtJet1 = jet->Pt();
623                         jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
624                         }
625                 if( nJets > 1 ) {
626                         index2 = 1;
627                         AliAODJet* jet = (AliAODJet*)arrayJets->At(1);
628                         maxPtJet2 = jet->Pt();
629                         jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
630                         }
631                 if( nJets > 2 ) {
632                         index3 = 2;
633                         AliAODJet* jet = (AliAODJet*)arrayJets->At(2);
634                         maxPtJet3 = jet->Pt();
635                         jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
636                         }
637       
638                 arrayJets->Delete();
639                 delete arrayJets;
640                 }
641         }
642   
643
644   // 3) LEADING PARTICLE
645   if( fAnaType == 4 ){
646         TObjArray* tracks = SortChargedParticles();
647         if( tracks ) {
648                 nJets = tracks->GetEntriesFast();
649                 if( nJets > 0 ) {
650                         index1 = 0;
651                         AliAODTrack* jet = (AliAODTrack*)tracks->At(0);
652                         maxPtJet1 = jet->Pt();
653                         jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
654                         }
655                 tracks->Clear();
656                 delete tracks; 
657                 }
658
659         }
660   fHistos->FillHistogram("hNJets",nJets);
661    
662   return *jetVect;
663
664 }
665
666
667 //-------------------------------------------------------------------
668 void AliAnalyseUE::Initialize(AliAnalysisTaskUE& taskUE){
669    
670   //Get principal settings from current instance of UE analysis-task
671   fAnaType = taskUE.GetAnaTopology();         
672   fkAOD = taskUE.GetAOD();           
673   fConeRadius = taskUE.GetConeRadius();
674   fDebug = taskUE.GetDebugLevel();
675   fFilterBit = taskUE.GetFilterBit();
676   fJet1EtaCut = taskUE.GetJet1EtaCut();
677   fJet2DeltaPhiCut = taskUE.GetJet2DeltaPhiCut();
678   fJet2RatioPtCut = taskUE.GetJet2RatioPtCut();
679   fJet3PtCut = taskUE.GetJet3PtCut();
680   fOrdering = taskUE.GetPtSumOrdering() ;
681   fRegionType = taskUE.GetRegionType();
682   fSimulateChJetPt = taskUE.GetSimulateChJetPt();
683   fTrackEtaCut = taskUE.GetTrackEtaCut(); 
684   fTrackPtCut = taskUE.GetTrackPtCut();
685   fHistos = taskUE.fHistosUE;
686   fUseChargeHadrons = taskUE.GetUseChargeHadrons();
687   fUseChPartJet = taskUE.GetUseChPartJet();
688   fUsePositiveCharge = taskUE.GetUseNegativeChargeType();
689   fUseSingleCharge = taskUE.GetUseSingleCharge();
690   
691 }
692
693 //-------------------------------------------------------------------
694 Bool_t  AliAnalyseUE::VertexSelection(AliAODEvent *aod, Int_t tracks, Double_t zed ){
695
696   //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
697   Int_t nVertex = aod->GetNumberOfVertices();
698   if( nVertex > 0 ) { // Only one vertex (reject pileup)
699         AliAODVertex* vertex = (AliAODVertex*)aod->GetPrimaryVertex();
700         Int_t nTracksPrim = vertex->GetNContributors();
701         Double_t zVertex = vertex->GetZ();
702         if (fDebug > 1) AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
703         // Select a quality vertex by number of tracks?
704         if( nTracksPrim < tracks || TMath::Abs(zVertex) > zed ) {
705                 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
706                 return kFALSE;
707                 }
708         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
709         } else {
710                 if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
711                 return kFALSE;
712                 }
713
714   return kTRUE;
715 }
716
717 //-------------------------------------------------------------------
718 Bool_t  AliAnalyseUE::VertexSelectionOld(AliAODEvent *aod ){
719
720   AliKFVertex primVtx(*(aod->GetPrimaryVertex()));
721   Int_t nTracksPrim=primVtx.GetNContributors();
722     if (fDebug > 1) AliInfo(Form(" Primary-vertex Selection: %d",nTracksPrim));
723     if(!nTracksPrim){
724         if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
725         return kFALSE;
726         }
727     if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED ...");
728
729   return kTRUE;
730 }
731
732 // PRIVATE METHODS **************************************************
733
734 TObjArray*  AliAnalyseUE::FindChargedParticleJets( Double_t chJetPtMin )
735 {
736   // Return a TObjArray of "charged particle jets"
737   
738   // Charged particle jet definition from reference:
739   // "Charged jet evolution and the underlying event
740   //  in proton-antiproton collisions at 1.8 TeV"
741   //  PHYSICAL REVIEW D 65 092002, CDF Collaboration
742   
743   // We defined "jets" as circular regions in eta-phi space with
744   // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
745   // Our jet algorithm is as follows:
746   //   1- Order all charged particles according to their pT .
747   //   2- Start with the highest pT particle and include in the jet all
748   //      particles within the radius R=0.7 considering each particle
749   //      in the order of decreasing pT and recalculating the centroid
750   //      of the jet after each new particle is added to the jet .
751   //   3- Go to the next highest pT particle not already included in
752   //      a jet and add to the jet all particles not already included in
753   //      a jet within R=0.7.
754   //   4- Continue until all particles are in a jet.
755   // We defined the transverse momentum of the jet to be
756   // the scalar pT sum of all the particles within the jet, where pT
757   // is measured with respect to the beam axis
758   
759   //  1 - Order all charged particles according to their pT .
760   Int_t nTracks = fkAOD->GetNTracks();
761   if( !nTracks ) return 0;
762   TObjArray tracks(nTracks);
763   
764   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
765         AliAODTrack* part = fkAOD->GetTrack( ipart );
766         if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
767         if( !part->Charge() ) continue;
768         if( part->Pt() < fTrackPtCut ) continue;
769         tracks.AddLast(part);
770         }
771   QSortTracks( tracks, 0, tracks.GetEntriesFast() );
772   
773   nTracks = tracks.GetEntriesFast();
774   if( !nTracks ) return 0;
775
776   TObjArray *jets = new TObjArray(nTracks);
777   TIter itrack(&tracks);
778   while( nTracks ) {
779         // 2- Start with the highest pT particle ...
780         Float_t px,py,pz,pt; 
781         AliAODTrack* track = (AliAODTrack*)itrack.Next();
782         if( !track ) continue;
783         px = track->Px();
784         py = track->Py();
785         pz = track->Pz();
786         pt = track->Pt(); // Use the energy member to store Pt
787         jets->AddLast( new TLorentzVector(px, py, pz, pt) );
788         tracks.Remove( track );
789         TLorentzVector* jet = (TLorentzVector*)jets->Last();
790         jet->SetPtEtaPhiE( 1., jet->Eta(), jet->Phi(), pt );
791         // 3- Go to the next highest pT particle not already included...
792         AliAODTrack* track1;
793         Double_t fPt = jet->E();
794         while ( (track1  = (AliAODTrack*)(itrack.Next())) ) {
795                 Double_t tphi = track1->Phi(); // here Phi is from 0 <-> 2Pi
796                 if (tphi > TMath::Pi()) tphi -= 2. * TMath::Pi(); // convert to  -Pi <-> Pi
797                 Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-tphi);
798                 Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +dphi*dphi );
799                 if( r < fConeRadius ) {
800                         fPt   = jet->E()+track1->Pt();  // Scalar sum of Pt
801                         // recalculating the centroid
802                         Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
803                         Double_t phi = jet->Phi()*jet->E()/fPt + tphi*track1->Pt()/fPt;
804                         jet->SetPtEtaPhiE( 1., eta, phi, fPt );
805                         tracks.Remove( track1 );
806                         }
807                 }
808     
809         tracks.Compress();
810         nTracks = tracks.GetEntries();
811         //   4- Continue until all particles are in a jet.
812         itrack.Reset();
813         } // end while nTracks
814   
815   // Convert to AODjets....
816   Int_t njets = jets->GetEntriesFast();
817   TObjArray* aodjets = new TObjArray(njets);
818   aodjets->SetOwner(kTRUE);
819   for(Int_t ijet=0; ijet<njets; ++ijet) {
820         TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
821         if (jet->E() < chJetPtMin) continue;
822         Float_t px, py,pz,en; // convert to 4-vector
823         px = jet->E() * TMath::Cos(jet->Phi());  // Pt * cos(phi)
824         py = jet->E() * TMath::Sin(jet->Phi());  // Pt * sin(phi)
825         pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
826         en = TMath::Sqrt(px * px + py * py + pz * pz);
827
828         aodjets->AddLast( new AliAODJet(px, py, pz, en) );
829         }
830   jets->Delete();
831   delete jets;
832   
833   // Order jets according to their pT .
834   QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
835   
836   // debug
837   if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
838   
839   return aodjets;
840 }
841
842
843 //____________________________________________________________________
844 void AliAnalyseUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
845 {
846
847   // Fill average particle Pt of control regions
848   
849   // Max cone
850   fHistos->FillHistogram("hRegionAvePartPtMaxVsEt", leadingE, ptMax);
851   // Min cone
852   fHistos->FillHistogram("hRegionAvePartPtMinVsEt", leadingE, ptMin);
853   // MAke distributions for UE comparison with MB data
854   fHistos->FillHistogram("hMinRegAvePt", ptMin);
855
856 }
857
858 //____________________________________________________________________
859 void AliAnalyseUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin  )
860 {
861
862   // Fill Nch multiplicity of control regions
863   
864   // Max cone
865   fHistos->FillHistogram("hRegionMultMaxVsEt", leadingE, nTrackPtmax);
866   fHistos->FillHistogram("hRegionMultMax",  nTrackPtmax);
867   
868   // Min cone
869   fHistos->FillHistogram("hRegionMultMinVsEt", leadingE, nTrackPtmin );
870   fHistos->FillHistogram("hRegionMultMin", nTrackPtmin);
871   
872   // MAke distributions for UE comparison with MB data
873   fHistos->FillHistogram("hMinRegSumPtvsMult", nTrackPtmin,ptMin);
874
875 }
876
877 //____________________________________________________________________
878 void AliAnalyseUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
879 {
880   // Fill sumPt of control regions
881   
882   // Max cone
883   fHistos->FillHistogram("hRegionSumPtMaxVsEt", leadingE, ptMax);
884   
885   // Min cone
886   fHistos->FillHistogram("hRegionSumPtMinVsEt", leadingE, ptMin);
887   
888   // MAke distributions for UE comparison with MB data
889   fHistos->FillHistogram("hMinRegSumPt", ptMin);
890   fHistos->FillHistogram("hMinRegSumPtJetPtBin", leadingE, ptMin);
891   fHistos->FillHistogram("hMaxRegSumPtJetPtBin", leadingE, ptMax);
892
893 }
894
895 //-------------------------------------------------------------------
896 Int_t AliAnalyseUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect, Int_t conePosition) 
897 {  
898   // return de region in delta phi
899   // -1 negative delta phi 
900   //  1 positive delta phi
901   //  0 outside region
902   static const Double_t k60rad  = 60.*TMath::Pi()/180.;
903   static const Double_t k120rad = 120.*TMath::Pi()/180.;
904  
905   Int_t region = 0;
906   if( fRegionType == 1 ) {
907         if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
908         // transverse regions
909         if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1; //left
910         if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1;    //right
911         if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) < k60rad ) region = 2;    //forward
912         if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) > k120rad ) region = -2; //backward
913     
914         } else if( fRegionType == 2 ) {
915         // Cone regions
916         Double_t deltaR = 0.;
917     
918         TVector3 positVect,negatVect;
919         if (conePosition==1){
920                 positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
921                 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
922         }else if (conePosition==2){
923                 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
924                 positVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()+TMath::PiOver2());
925                 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()-TMath::PiOver2());
926         }else if (conePosition==3){
927                 if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
928                 Double_t weightEta = jetVect[0].Eta() * jetVect[0].Pt()/(jetVect[0].Pt() + jetVect[1].Pt()) + 
929                 jetVect[1].Eta() * jetVect[1].Pt()/(jetVect[0].Pt() + jetVect[1].Pt());
930                 //Double_t weightEta = jetVect[0].Eta() * jetVect[0].Mag()/(jetVect[0].Mag() + jetVect[1].Mag()) + 
931                 // jetVect[1].Eta() * jetVect[1].Mag()/(jetVect[0].Mag() + jetVect[1].Mag());
932                 positVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()+TMath::PiOver2());
933                 negatVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()-TMath::PiOver2());
934                 }
935   if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) { 
936         region = 1;  
937         deltaR = positVect.DrEtaPhi(*partVect);
938   } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) { 
939         region = -1;  
940         deltaR = negatVect.DrEtaPhi(*partVect);
941         }
942     
943   if (deltaR > fConeRadius) region = 0;
944     
945   } else
946         AliError("Unknow region type");
947         
948         return region;
949   }
950
951
952
953 //-------------------------------------------------------------------
954 void  AliAnalyseUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
955 {
956   // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
957   
958   static TObject *tmp;
959   static int i;           // "static" to save stack space
960   int j;
961   
962   while (last - first > 1) {
963     i = first;
964     j = last;
965     for (;;) {
966       while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
967         ;
968       while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
969         ;
970       if (i >= j)
971         break;
972       
973       tmp  = a[i];
974       a[i] = a[j];
975       a[j] = tmp;
976     }
977     if (j == first) {
978       ++first;
979       continue;
980     }
981     tmp = a[first];
982     a[first] = a[j];
983     a[j] = tmp;
984     if (j - first < last - (j + 1)) {
985       QSortTracks(a, first, j);
986       first = j + 1;   // QSortTracks(j + 1, last);
987     } else {
988       QSortTracks(a, j + 1, last);
989       last = j;        // QSortTracks(first, j);
990     }
991   }
992 }
993
994
995
996
997 //-------------------------------------------------------------------
998 void AliAnalyseUE::SetRegionArea(TVector3 *jetVect)
999 {
1000   // Set region area
1001   Double_t areaCorrFactor=0.;
1002   Double_t deltaEta = 0.;
1003   if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.;
1004   else if (fRegionType==2){ 
1005         deltaEta = 0.9-TMath::Abs(jetVect[0].Eta());
1006         if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius;
1007         else{
1008                 areaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
1009                 (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
1010                 fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-areaCorrFactor;
1011                 }
1012         }else AliWarning("Unknown Region Type");
1013   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));
1014 }
1015
1016
1017 //____________________________________________________________________
1018 TObjArray*  AliAnalyseUE::SortChargedParticles()
1019 {
1020   //  return an array with all charged particles ordered according to their pT .
1021   Int_t nTracks = fkAOD->GetNTracks();
1022   if( !nTracks ) return 0;
1023   TObjArray* tracks = new TObjArray(nTracks);
1024
1025   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
1026         AliAODTrack* part = fkAOD->GetTrack( ipart );
1027         if( !part->TestFilterBit( fFilterBit ) ) continue; // track cut selection
1028         if( !part->Charge() ) continue;
1029         if( part->Pt() < fTrackPtCut ) continue;
1030         tracks->AddLast( part );
1031         }
1032   QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
1033
1034   nTracks = tracks->GetEntriesFast();
1035   if( !nTracks ) return 0;
1036
1037   return tracks;
1038   }
1039
1040
1041 //____________________________________________________________________
1042 Bool_t AliAnalyseUE::TrackMCSelected(Double_t charge, Double_t pT, Double_t eta, Int_t pdgCode)const{
1043   
1044   // MC track selection 
1045   Double_t const kMyTolerance = 0.0000001;
1046   //if (charge == 0. || charge == -99.) return kFALSE;
1047   if (charge < kMyTolerance || charge + 99 < kMyTolerance) return kFALSE;
1048         
1049   if ( fUseSingleCharge ) { // Charge selection
1050         if ( fUsePositiveCharge && charge < 0.) return kFALSE; // keep Positives
1051         if ( !fUsePositiveCharge && charge > 0.) return kFALSE; // keep Negatives
1052         }
1053         
1054   //Kinematics cuts on particle
1055   if ((pT < fTrackPtCut) || (TMath::Abs(eta) > fTrackEtaCut )) return kFALSE;
1056         
1057   Bool_t isHadron = TMath::Abs(pdgCode)==211 ||
1058         TMath::Abs(pdgCode)==2212 ||
1059         TMath::Abs(pdgCode)==321;
1060         
1061   if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1062         
1063   return kTRUE;
1064
1065 }
1066
1067 //____________________________________________________________________
1068 Bool_t AliAnalyseUE::TrackSelected(AliAODTrack* part)const{
1069
1070   // Real track selection
1071   if ( !part->TestFilterBit(fFilterBit) ) return kFALSE; // track cut selection
1072   if ( !part->IsPrimaryCandidate()) return kFALSE; // reject whatever is not linked to collision point
1073   // PID Selection: Reject everything but hadrons
1074   Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion || 
1075         part->GetMostProbablePID()==AliAODTrack::kKaon || 
1076         part->GetMostProbablePID()==AliAODTrack::kProton;
1077   if ( fUseChargeHadrons && !isHadron ) return kFALSE;
1078       
1079   if ( !part->Charge() ) return kFALSE; //Only charged
1080   if ( fUseSingleCharge ) { // Charge selection
1081         if ( fUsePositiveCharge && part->Charge() < 0.) return kFALSE; // keep Positives
1082         if ( !fUsePositiveCharge && part->Charge() > 0.) return kFALSE; // keep Negatives
1083         } 
1084     
1085   if ( part->Pt() < fTrackPtCut ) return kFALSE;
1086   if( TMath::Abs(part->Eta()) > fTrackEtaCut ) return kFALSE;
1087
1088   return kTRUE;
1089 }
1090
1091