Added PID task, some fixes for coding conventions
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskUE.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id:$ */
17
18 #include <TROOT.h>
19 #include <TSystem.h>
20 #include <TChain.h>
21 #include <TFile.h>
22 #include <TList.h>
23 #include <TH1I.h>
24 #include <TH1F.h>
25 #include <TH2F.h>
26 #include <TCanvas.h>
27 #include <TVector3.h>
28 #include <TLorentzVector.h>
29 #include <TMath.h>
30 #include <TTree.h>
31 #include <TBranch.h>
32
33 #include "AliAnalysisTaskUE.h"
34 #include "AliAnalysisManager.h"
35 #include "AliMCEventHandler.h"
36 #include "AliAODEvent.h"
37 #include "AliAODInputHandler.h"
38 #include "AliAODHandler.h"
39 #include "AliStack.h"
40 #include "AliAODJet.h"
41 #include "AliAODTrack.h"
42
43 #include "AliLog.h"
44
45 //
46 // Analysis class for Underlying Event studies
47 //
48 // Look for correlations on the tranverse regions to 
49 // the leading charged jet
50 //
51 // This class needs as input AOD with track and Jets
52 // the output is a list of histograms
53 //
54 // AOD can be either connected to the InputEventHandler  
55 // for a chain of AOD files 
56 // or 
57 // to the OutputEventHandler
58 // for a chain of ESD files, so this case class should be 
59 // in the train after the Jet finder
60 //
61 //    Arian.Abrahantes.Quintana@cern.ch 
62 //    Ernesto.Lopez.Torres@cern.ch
63 // 
64
65
66 ClassImp( AliAnalysisTaskUE)
67
68 ////////////////////////////////////////////////////////////////////////
69
70
71 //____________________________________________________________________
72 AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name):
73 AliAnalysisTask(name, ""),
74 fDebug(kFALSE),          
75 fAOD(0x0),            
76 fAODjets(0x0),
77 fListOfHistos(0x0),  
78 fBinsPtInHist(30),     
79 fMinJetPtInHist(0.),
80 fMaxJetPtInHist(300.),  
81 fAnaType(1),         
82 fRegionType(1),
83 fConeRadius(0.7),
84 fAreaReg(1.5393), // Pi*0.7*0.7
85 fUseChPartJet(kFALSE),
86 fUseSingleCharge(kFALSE),
87 fUsePositiveCharge(kTRUE),
88 fOrdering(1),
89 fFilterBit(0xFF),
90 fJetsOnFly(kFALSE),
91 fJet1EtaCut(0.2),
92 fJet2DeltaPhiCut(2.616),    // 150 degrees
93 fJet2RatioPtCut(0.8),
94 fJet3PtCut(15.),
95 fTrackPtCut(0.),
96 fTrackEtaCut(0.9),
97 fhNJets(0x0),
98 fhEleadingPt(0x0),
99 fhMinRegPtDist(0x0),
100 fhRegionMultMin(0x0),
101 fhMinRegAvePt(0x0), 
102 fhMinRegSumPt(0x0),            
103 fhMinRegMaxPtPart(0x0),
104 fhMinRegSumPtvsMult(0x0),
105 fhdNdEtaPhiDist(0x0),        
106 fhFullRegPartPtDistVsEt(0x0), 
107 fhTransRegPartPtDistVsEt(0x0),
108 fhRegionSumPtMaxVsEt(0x0),
109 fhRegionMultMax(0x0),         
110 fhRegionMultMaxVsEt(0x0),     
111 fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0),         
112 fhRegionMultMinVsEt(0x0),     
113 fhRegionAveSumPtVsEt(0x0),    
114 fhRegionDiffSumPtVsEt(0x0),   
115 fhRegionAvePartPtMaxVsEt(0x0),
116 fhRegionAvePartPtMinVsEt(0x0),
117 fhRegionMaxPartPtMaxVsEt(0x0),
118 fSettingsTree(0x0)//,   fhValidRegion(0x0)
119 {
120   // Default constructor
121   // Define input and output slots here
122   // Input slot #0 works with a TChain
123   DefineInput(0, TChain::Class());
124   // Output slot #0 writes into a TList container
125   DefineOutput(0, TList::Class());
126 }
127
128 //____________________________________________________________________
129 void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/)
130 {
131   // Connect the input data  
132   
133   // We need AOD with tracks and jets.
134   // Since AOD can be either connected to the InputEventHandler (input chain fron AOD files)
135   // or to the OutputEventHandler ( AOD is created by a previus task in the train )
136   // we need to check where it is and get the pointer to AODEvent in the right way
137   
138   if (fDebug > 1) AliInfo("ConnectInputData() \n");
139   
140   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
141   
142   if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
143     fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
144     if (fDebug > 1) AliInfo("  ==== Tracks from AliAODInputHandler");
145     // Case when jets are reconstructed on the fly from AOD tracks
146     // (the Jet Finder is using the AliJetAODReader) of InputEventHandler
147     // and put in the OutputEventHandler AOD. Useful whe you want to reconstruct jets with
148     // different parameters to default ones stored in the AOD or to use a different algorithm
149     if( fJetsOnFly ) {
150       handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
151       if( handler && handler->InheritsFrom("AliAODHandler") ) {
152         fAODjets = ((AliAODHandler*)handler)->GetAOD();
153         if (fDebug > 1) AliInfo("  ==== Jets from AliAODHandler");
154       }
155     } else {
156       fAODjets = fAOD;
157       if (fDebug > 1) AliInfo("  ==== Jets from AliAODInputHandler");
158     }
159   } else {
160     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
161     if( handler && handler->InheritsFrom("AliAODHandler") ) {
162       fAOD  = ((AliAODHandler*)handler)->GetAOD();
163       fAODjets = fAOD;
164       if (fDebug > 1) AliInfo("  ==== Tracks and Jets from AliAODHandler");
165     } else {
166       AliFatal("I can't get any AOD Event Handler");
167       return;
168     }
169   }
170 }
171
172 //____________________________________________________________________
173 void  AliAnalysisTaskUE::CreateOutputObjects()
174 {
175   // Create the output container
176   //
177   if (fDebug > 1) AliInfo("CreateOutPutData()");
178   //
179   //  Histograms
180   CreateHistos();
181   fListOfHistos->SetOwner(kTRUE);  
182   //  OpenFile(0);
183 }
184
185
186 //____________________________________________________________________
187 void  AliAnalysisTaskUE::Exec(Option_t */*option*/)
188 {
189   // Execute analysis for current event
190   //
191   if ( fDebug > 3 ) AliInfo( " Processing event..." );
192   AnalyseUE();
193   
194   // Post the data
195   PostData(0, fListOfHistos);
196 }
197
198 //____________________________________________________________________
199 void  AliAnalysisTaskUE::AnalyseUE()
200 {
201   //
202   // Look for correlations on the tranverse regions to 
203   // the leading charged jet
204   // 
205   
206   // ------------------------------------------------
207   // Find Leading Jets 1,2,3 
208   // (could be skipped if Jets are sort by Pt...)
209   Double_t maxPtJet1 = 0.; 
210   Int_t    index1 = -1;
211   Double_t maxPtJet2 = 0.;   // jet 2 need for back to back inclusive
212   Int_t    index2 = -1;
213   Double_t maxPtJet3 = 0.;   // jet 3 need for back to back exclusive
214   Int_t    index3 = -1;
215   TVector3 jetVect[3];
216   Int_t nJets = 0;
217   
218   if( !fUseChPartJet ) {
219     
220     // Use AOD Jets
221     nJets = fAODjets->GetNJets();
222     //  printf("AOD %d jets \n", nJets);
223     for( Int_t i=0; i<nJets; ++i ) {
224       AliAODJet* jet = fAODjets->GetJet(i);
225       Double_t jetPt = jet->Pt();//*1.666;  // FIXME Jet Pt Correction ?????!!!
226       if( jetPt > maxPtJet1 ) {
227         maxPtJet3 = maxPtJet2; index3 = index2;
228         maxPtJet2 = maxPtJet1; index2 = index1;
229         maxPtJet1 = jetPt; index1 = i;
230       } else if( jetPt > maxPtJet2 ) {
231         maxPtJet3 = maxPtJet2; index3 = index2;
232         maxPtJet2 = jetPt; index2 = i;
233       } else if( jetPt > maxPtJet3 ) {
234         maxPtJet3 = jetPt; index3 = i;
235       }
236     }
237     if( index1 != -1 ) {
238       AliAODJet* jet = fAODjets->GetJet(index1);
239       jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
240     }
241     if( index2 != -1 ) {
242       AliAODJet* jet = fAODjets->GetJet(index2);
243       jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
244     }
245     if( index3 != -1 ) {
246       AliAODJet* jet = fAODjets->GetJet(index3);
247       jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
248     }
249     
250   } else {
251     
252     // Use "Charged Particle Jets"
253     TObjArray* jets = FindChargedParticleJets();
254     if( jets ) {
255       nJets = jets->GetEntriesFast();
256       if( nJets > 0 ) {
257         index1 = 0;
258         AliAODJet* jet = (AliAODJet*)jets->At(0);
259         maxPtJet1 = jet->Pt();
260         jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
261       }
262       if( nJets > 1 ) {
263         index2 = 1;
264         AliAODJet* jet = (AliAODJet*)jets->At(1);
265         maxPtJet1 = jet->Pt();
266         jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
267       }
268       if( nJets > 2 ) {
269         index3 = 2;
270         AliAODJet* jet = (AliAODJet*)jets->At(2);
271         maxPtJet1 = jet->Pt();
272         jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
273       }
274       
275       jets->Delete();
276       delete jets;
277     }
278   }
279   
280   
281   fhNJets->Fill(nJets);
282   
283   if( fDebug > 1 ) {
284     if( index1 < 0 ) {
285       AliInfo("\n   Skipping Event, not jet found...");
286       return;
287     } else
288       AliInfo(Form("\n   Pt Leading Jet = %6.1f eta=%5.3f ",  maxPtJet1, jetVect[0].Eta() ));
289   }
290   
291   
292   // ----------------------------------------------
293   // Cut events by jets topology
294   // fAnaType:
295   //     1 = inclusive,
296   //         - Jet1 |eta| < fJet1EtaCut
297   //     2 = back to back inclusive
298   //         - fulfill case 1
299   //         - |Jet1.Phi - Jet2.Phi| > fJet2DeltaPhiCut
300   //         - Jet2.Pt/Jet1Pt > fJet2RatioPtCut
301   //     3 = back to back exclusive
302   //         - fulfill case 2
303   //         - Jet3.Pt < fJet3PtCut
304   
305   if( index1 < 0 || TMath::Abs(jetVect[0].Eta()) > fJet1EtaCut) {
306     if( fDebug > 1 ) AliInfo("\n   Skipping Event...Jet1 |eta| > fJet1EtaCut");
307     return;
308   }
309   // back to back inclusive
310   if( fAnaType > 1 && index2 == -1 ) {
311     if( fDebug > 1 ) AliInfo("\n   Skipping Event... no second Jet found");
312     return;
313   }
314   if( fAnaType > 1 && index2 > -1 ) {
315     if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
316        maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) {
317       if( fDebug > 1 ) AliInfo("\n   Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
318       return;
319     }
320   }
321   // back to back exclusive
322   if( fAnaType > 2 && index3 > -1 ) {
323     if( maxPtJet3 > fJet3PtCut ) {
324       if( fDebug > 1 ) AliInfo("\n   Skipping Event... Jet3.Pt > fJet3PtCut ");
325       return;
326     }
327   }
328   
329   fhEleadingPt->Fill( maxPtJet1 );
330   //Area for Normalization Purpose at Display histos
331   SetRegionArea(jetVect);
332   
333   // ----------------------------------------------
334   // Find max and min regions
335   Double_t sumPtRegionPosit = 0.;
336   Double_t sumPtRegionNegat = 0.;
337   Double_t maxPartPtRegion  = 0.;
338   Int_t    nTrackRegionPosit = 0;
339   Int_t    nTrackRegionNegat = 0;
340   static const Double_t k270rad = 270.*TMath::Pi()/180.;
341   
342   Int_t nTracks = fAOD->GetNTracks();
343   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
344     AliAODTrack* part = fAOD->GetTrack( ipart );
345     if ( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
346     // PID Selection: Reject everything but hadrons
347     Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion || 
348                       part->GetMostProbablePID()==AliAODTrack::kKaon || 
349                       part->GetMostProbablePID()==AliAODTrack::kProton;
350     if ( !isHadron ) continue;
351     
352     if ( !part->Charge() ) continue; //Only charged
353     if ( fUseSingleCharge ) { // Charge selection
354       if ( fUsePositiveCharge && part->Charge() < 0.) continue; // keep Positives
355       if ( !fUsePositiveCharge && part->Charge() > 0.) continue; // keep Negatives
356     }
357     
358     if ( part->Pt() < fTrackPtCut ) continue;
359     
360     TVector3 partVect(part->Px(), part->Py(), part->Pz());
361     
362     Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
363     if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
364     fhdNdEtaPhiDist->Fill( deltaPhi );
365     fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
366     
367     Int_t region = IsTrackInsideRegion( jetVect, &partVect );  
368     
369     if (region > 0) {
370       if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
371       sumPtRegionPosit += part->Pt();
372       nTrackRegionPosit++;
373       fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
374     }
375     if (region < 0) {
376       if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
377       sumPtRegionNegat += part->Pt();
378       nTrackRegionNegat++;
379       fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
380     }
381   }
382   
383   //How quantities will be sorted before Fill Min and Max Histogram
384   //  1=Plots will be CDF-like
385   //  2=Plots will be Marchesini-like
386   if( fOrdering == 1 ) {
387     if( sumPtRegionPosit > sumPtRegionNegat ) {
388       FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
389     } else {
390       FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
391     }
392     if (nTrackRegionPosit > nTrackRegionNegat ) {
393       FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
394     } else {
395       FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
396     }
397   } else if( fOrdering == 2 ) {
398     if (sumPtRegionPosit > sumPtRegionNegat) {
399       FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
400       FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
401     } else {
402       FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
403       FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
404     }
405   }
406   
407   Double_t avePosRegion = (nTrackRegionPosit) ? sumPtRegionPosit/nTrackRegionPosit : 0.;
408   Double_t aveNegRegion = (nTrackRegionNegat) ? sumPtRegionNegat/nTrackRegionNegat : 0.;
409   if( avePosRegion > aveNegRegion ) {
410     FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg );
411   } else {
412     FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg );
413   }
414   
415   fhRegionMaxPartPtMaxVsEt->Fill(maxPtJet1, maxPartPtRegion );
416   
417   // Compute pedestal like magnitudes
418   fhRegionDiffSumPtVsEt->Fill(maxPtJet1, TMath::Abs(sumPtRegionPosit-sumPtRegionNegat)/(2.0*fAreaReg));
419   fhRegionAveSumPtVsEt->Fill(maxPtJet1, (sumPtRegionPosit+sumPtRegionNegat)/(2.0*fAreaReg));
420   
421 }
422
423 //____________________________________________________________________
424 void AliAnalysisTaskUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
425 {
426   // Fill sumPt of control regions
427   
428   // Max cone
429   fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax );
430   // Min cone
431   fhRegionSumPtMinVsEt->Fill( leadingE, ptMin );
432   // MAke distributions for UE comparison with MB data
433   fhMinRegSumPt->Fill(ptMin);
434 }
435
436 //____________________________________________________________________
437 void AliAnalysisTaskUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
438 {
439   // Fill average particle Pt of control regions
440   
441   // Max cone
442   fhRegionAvePartPtMaxVsEt->Fill( leadingE, ptMax );
443   // Min cone
444   fhRegionAvePartPtMinVsEt->Fill( leadingE, ptMin );
445   // MAke distributions for UE comparison with MB data
446   fhMinRegAvePt->Fill(ptMin);
447 }
448
449 //____________________________________________________________________
450 void AliAnalysisTaskUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin  )
451 {
452   // Fill Nch multiplicity of control regions
453   
454   // Max cone
455   fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax );
456   fhRegionMultMax->Fill( nTrackPtmax );
457   // Min cone
458   fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin );
459   fhRegionMultMin->Fill( nTrackPtmin );
460   // MAke distributions for UE comparison with MB data
461   fhMinRegSumPtvsMult->Fill(nTrackPtmin,ptMin);
462 }
463
464 //____________________________________________________________________
465 Int_t AliAnalysisTaskUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect) 
466 {  
467   // return de region in delta phi
468   // -1 negative delta phi 
469   //  1 positive delta phi
470   //  0 outside region
471   static const Double_t k60rad  = 60.*TMath::Pi()/180.;
472   static const Double_t k120rad = 120.*TMath::Pi()/180.;
473   
474   Int_t region = 0;
475   if( fRegionType == 1 ) {
476     if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
477     // transverse regions
478     if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1;
479     if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1;
480     
481   } else if( fRegionType == 2 ) {
482     // Cone regions
483     Double_t deltaR = 0.;
484     
485     TVector3 positVect,negatVect;
486     positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
487     negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
488     
489     if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) { 
490       region = 1;  
491       deltaR = positVect.DrEtaPhi(*partVect);
492     } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) { 
493       region = -1;  
494       deltaR = negatVect.DrEtaPhi(*partVect);
495     }
496     
497     if (deltaR > fConeRadius) region = 0;
498     
499   } else 
500     AliError("Unknow region type");
501   
502   // For debug (to be removed)
503   //if( region != 0 )  fhValidRegion->Fill( partVect->Eta()-jetVect[0].Eta(), jetVect[0].DeltaPhi(*partVect) );
504   
505   return region;
506 }
507
508
509 //____________________________________________________________________
510 TObjArray*  AliAnalysisTaskUE::FindChargedParticleJets()
511 {
512   // Return a TObjArray of "charged particle jets"
513   //
514   // Charged particle jet definition from reference:
515   //
516   // "Charged jet evolution and the underlying event
517   //  in proton-antiproton collisions at 1.8 TeV"
518   //  PHYSICAL REVIEW D 65 092002, CDF Collaboration
519   //
520   // We define "jets" as circular regions in eta-phi space with
521   // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
522   // Our jet algorithm is as follows:
523   //   1- Order all charged particles according to their pT .
524   //   2- Start with the highest pT particle and include in the jet all
525   //      particles within the radius R=0.7 considering each particle
526   //      in the order of decreasing pT and recalculating the centroid
527   //      of the jet after each new particle is added to the jet .
528   //   3- Go to the next highest pT particle not already included in
529   //      a jet and add to the jet all particles not already included in
530   //      a jet within R=0.7.
531   //   4- Continue until all particles are in a jet.
532   // We define the transverse momentum of the jet to be
533   // the scalar pT sum of all the particles within the jet, where pT
534   // is measured with respect to the beam axis
535   
536   //  1 - Order all charged particles according to their pT .
537   Int_t nTracks = fAOD->GetNTracks();
538   if( !nTracks ) return 0;
539   TObjArray tracks(nTracks);
540   
541   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
542     AliAODTrack* part = fAOD->GetTrack( ipart );
543     if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
544     if( !part->Charge() ) continue;
545     if( part->Pt() < fTrackPtCut ) continue;
546     tracks.AddLast(part);
547   }
548   QSortTracks( tracks, 0, tracks.GetEntriesFast() );
549   
550   nTracks = tracks.GetEntriesFast();
551   if( !nTracks ) return 0;
552   TObjArray *jets = new TObjArray(nTracks);
553   TIter itrack(&tracks);
554   while( nTracks ) {
555     // 2- Start with the highest pT particle ...
556     Float_t px,py,pz,pt; 
557     AliAODTrack* track = (AliAODTrack*)itrack.Next();
558     if( !track ) continue;
559     px = track->Px();
560     py = track->Py();
561     pz = track->Pz();
562     pt = track->Pt(); // Use the energy member to store Pt
563     jets->AddLast( new TLorentzVector(px, py, pz, pt) );
564     tracks.Remove( track );
565     TLorentzVector* jet = (TLorentzVector*)jets->Last();
566     // 3- Go to the next highest pT particle not already included...
567     AliAODTrack* track1;
568     while ( (track1  = (AliAODTrack*)(itrack.Next())) ) {
569       Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-track1->Phi());
570       Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +
571                                dphi*dphi );
572       if( r < fConeRadius ) {
573         Double_t fPt   = jet->E()+track1->Pt();  // Scalar sum of Pt
574         // recalculating the centroid
575         Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()/fPt;
576         Double_t phi = jet->Phi()*jet->E()/fPt + track1->Phi()/fPt;
577         jet->SetPtEtaPhiE( 1., eta, phi, fPt );
578         tracks.Remove( track1 );
579       }
580     }
581     
582     tracks.Compress();
583     nTracks = tracks.GetEntries();
584     //   4- Continue until all particles are in a jet.
585     itrack.Reset();
586   } // end while nTracks
587   
588   // Convert to AODjets....
589   Int_t njets = jets->GetEntriesFast();
590   TObjArray* aodjets = new TObjArray(njets);
591   aodjets->SetOwner(kTRUE);
592   for(Int_t ijet=0; ijet<njets; ++ijet) {
593     TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
594     Float_t px, py,pz,en; // convert to 4-vector
595     px = jet->E() * TMath::Cos(jet->Phi());  // Pt * cos(phi)
596     py = jet->E() * TMath::Sin(jet->Phi());  // Pt * sin(phi)
597     pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
598     en = TMath::Sqrt(px * px + py * py + pz * pz);
599     aodjets->AddLast( new AliAODJet(px, py, pz, en) );
600   }
601   // Order jets according to their pT .
602   QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
603   
604   // debug
605   if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
606   
607   return aodjets;
608 }
609
610 //____________________________________________________________________
611 void  AliAnalysisTaskUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
612 {
613   // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
614   
615   static TObject *tmp;
616   static int i;           // "static" to save stack space
617   int j;
618   
619   while (last - first > 1) {
620     i = first;
621     j = last;
622     for (;;) {
623       while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
624         ;
625       while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
626         ;
627       if (i >= j)
628         break;
629       
630       tmp  = a[i];
631       a[i] = a[j];
632       a[j] = tmp;
633     }
634     if (j == first) {
635       ++first;
636       continue;
637     }
638     tmp = a[first];
639     a[first] = a[j];
640     a[j] = tmp;
641     if (j - first < last - (j + 1)) {
642       QSortTracks(a, first, j);
643       first = j + 1;   // QSortTracks(j + 1, last);
644     } else {
645       QSortTracks(a, j + 1, last);
646       last = j;        // QSortTracks(first, j);
647     }
648   }
649 }
650
651 //____________________________________________________________________
652 void AliAnalysisTaskUE::SetRegionArea(TVector3 *jetVect)
653 {
654   Double_t fAreaCorrFactor=0.;
655   Double_t deltaEta = 0.;
656   if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.;
657   else if (fRegionType==2){ 
658     deltaEta = 0.9-TMath::Abs(jetVect[0].Eta());
659     if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius;
660     else{
661       fAreaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
662       (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
663       fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-fAreaCorrFactor;
664     }
665   }else AliWarning("Unknown Rgion Type");
666   if (fDebug>5) AliInfo(Form("\n dEta=%5.3f Angle =%5.3f Region Area = %5.3f Corr Factor=%5.4f \n",deltaEta,TMath::ACos(deltaEta/fConeRadius),fAreaReg,fAreaCorrFactor));
667 }
668
669 //____________________________________________________________________
670 void  AliAnalysisTaskUE::CreateHistos()
671 {
672   fListOfHistos = new TList();
673   
674   
675   fhNJets = new TH1F("fhNJets", "n Jet",  10, 0, 10);
676   fhNJets->SetXTitle("# of jets");
677   fhNJets->Sumw2();
678   fListOfHistos->Add( fhNJets );                 // At(0) 
679   
680   fhEleadingPt  = new TH1F("hEleadingPt",   "Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
681   fhEleadingPt->SetXTitle("P_{T} (GeV/c)");
682   fhEleadingPt->SetYTitle("dN/dP_{T}");
683   fhEleadingPt->Sumw2();
684   fListOfHistos->Add( fhEleadingPt );            // At(1)
685   
686   fhMinRegPtDist = new TH1F("hMinRegPtDist",   "P_{T} distribution in Min zone",  50,0.,20.);
687   fhMinRegPtDist->SetXTitle("P_{T} (GeV/c)");
688   fhMinRegPtDist->SetYTitle("dN/dP_{T}");
689   fhMinRegPtDist->Sumw2();
690   fListOfHistos->Add( fhMinRegPtDist );          // At(2)
691   
692   fhRegionMultMin = new TH1F("hRegionMultMin",      "N_{ch}^{90, min}",  21, -0.5,   20.5);
693   fhRegionMultMin->SetXTitle("N_{ch tracks}");
694   fhRegionMultMin->Sumw2();
695   fListOfHistos->Add( fhRegionMultMin );         // At(3)            
696   
697   fhMinRegAvePt = new TH1F("hMinRegAvePt", "#LTp_{T}#GT",  50, 0.,   20.);
698   fhMinRegAvePt->SetXTitle("P_{T} (GeV/c)");
699   fhMinRegAvePt->Sumw2();
700   fListOfHistos->Add( fhMinRegAvePt );           // At(4)
701   
702   fhMinRegSumPt = new TH1F("hMinRegSumPt", "#Sigma p_{T} ",  50, 0.,   20.);
703   fhMinRegSumPt->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");  
704   fhMinRegSumPt->SetXTitle("#Sigma p_{T} (GeV/c)");
705   fhMinRegSumPt->Sumw2();
706   fListOfHistos->Add( fhMinRegSumPt );           // At(5)
707   
708   fhMinRegMaxPtPart = new TH1F("hMinRegMaxPtPart", "max(p_{T})|_{event} ",  50, 0.,   20.);
709   fhMinRegMaxPtPart->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");  
710   fhMinRegMaxPtPart->SetXTitle("p_{T} (GeV/c)");
711   fhMinRegMaxPtPart->Sumw2();
712   fListOfHistos->Add( fhMinRegMaxPtPart );       // At(6)
713   
714   fhMinRegSumPtvsMult = new TH1F("hMinRegSumPtvsMult", "#Sigma p_{T} vs. Multiplicity ",  21, -0.5,   20.5);
715   fhMinRegSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
716   fhMinRegSumPtvsMult->SetXTitle("N_{charge}");
717   fhMinRegSumPtvsMult->Sumw2();
718   fListOfHistos->Add( fhMinRegSumPtvsMult );     // At(7);
719   
720   fhdNdEtaPhiDist  = new TH1F("hdNdEtaPhiDist",   "Charge particle density |#eta|< 1 vs #Delta#phi",  120, 0.,   2.*TMath::Pi());
721   fhdNdEtaPhiDist->SetXTitle("#Delta#phi");
722   fhdNdEtaPhiDist->SetYTitle("dN_{ch}/d#etad#phi");
723   fhdNdEtaPhiDist->Sumw2();
724   fListOfHistos->Add( fhdNdEtaPhiDist );        // At(8)
725   
726   // Can be use to get part pt distribution for differente Jet Pt bins
727   fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", "dN/dP_{T} |#eta|<1 vs Leading Jet P_{T}",  
728                                      50,0.,50., fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
729   fhFullRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
730   fhFullRegPartPtDistVsEt->SetXTitle("p_{T}");
731   fhFullRegPartPtDistVsEt->Sumw2();
732   fListOfHistos->Add( fhFullRegPartPtDistVsEt );  // At(9) 
733   
734   // Can be use to get part pt distribution for differente Jet Pt bins
735   fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", "dN/dP_{T} in tranvese regions |#eta|<1 vs Leading Jet P_{T}",  
736                                       50,0.,50., fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
737   fhTransRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
738   fhTransRegPartPtDistVsEt->SetXTitle("p_{T}");
739   fhTransRegPartPtDistVsEt->Sumw2();
740   fListOfHistos->Add( fhTransRegPartPtDistVsEt );  // At(10)
741   
742   
743   fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt",  "P_{T}^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
744   fhRegionSumPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
745   fhRegionSumPtMaxVsEt->Sumw2();
746   fListOfHistos->Add( fhRegionSumPtMaxVsEt );      // At(11)
747   
748   fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt",   "P_{T}^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
749   fhRegionSumPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
750   fhRegionSumPtMinVsEt->Sumw2();
751   fListOfHistos->Add( fhRegionSumPtMinVsEt );      // At(12)
752   
753   fhRegionMultMax = new TH1I("hRegionMultMax",      "N_{ch}^{90, max}",  21, -0.5,   20.5);
754   fhRegionMultMax->SetXTitle("N_{ch tracks}");
755   fhRegionMultMax->Sumw2();
756   fListOfHistos->Add( fhRegionMultMax );           // At(13)
757   
758   fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt",  "N_{ch}^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
759   fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)");
760   fhRegionMultMaxVsEt->Sumw2();
761   fListOfHistos->Add( fhRegionMultMaxVsEt );       // At(14)
762   
763   fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt",  "N_{ch}^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
764   fhRegionMultMinVsEt->SetXTitle("E (GeV/c)");
765   fhRegionMultMinVsEt->Sumw2();
766   fListOfHistos->Add( fhRegionMultMinVsEt );      // At(15)
767   
768   fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
769   fhRegionAveSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
770   fhRegionAveSumPtVsEt->Sumw2();
771   fListOfHistos->Add( fhRegionAveSumPtVsEt );     // At(16)
772   
773   fhRegionDiffSumPtVsEt= new TH1F("hRegionPtDiffVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
774   fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
775   fhRegionDiffSumPtVsEt->Sumw2();
776   fListOfHistos->Add( fhRegionDiffSumPtVsEt );    // At(17)
777   
778   fhRegionAvePartPtMaxVsEt = new TH1F("hRegionAvePartPtMaxVsEt", "#LTp_{T}#GT^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
779   fhRegionAvePartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
780   fhRegionAvePartPtMaxVsEt->Sumw2();
781   fListOfHistos->Add( fhRegionAvePartPtMaxVsEt );  // At(18)
782   
783   fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
784   fhRegionAvePartPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
785   fhRegionAvePartPtMinVsEt->Sumw2();
786   fListOfHistos->Add( fhRegionAvePartPtMinVsEt );   // At(19)
787   
788   fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
789   fhRegionMaxPartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
790   fhRegionMaxPartPtMaxVsEt->Sumw2();
791   fListOfHistos->Add( fhRegionMaxPartPtMaxVsEt );    // At(20)
792   
793   fSettingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
794   fListOfHistos->Add( fSettingsTree );    // At(21)
795   
796   /*   
797    // For debug region selection
798    fhValidRegion = new TH2F("hValidRegion", "dN/d#eta/d#phi",      
799    20, -2.,2., 62, -TMath::Pi(),   TMath::Pi());
800    fhValidRegion->SetYTitle("#Delta#phi");
801    fhValidRegion->Sumw2();
802    fListOfHistos->Add( fhValidRegion );  // At(15)
803    */
804 }
805
806 //____________________________________________________________________
807 void  AliAnalysisTaskUE::Terminate(Option_t */*option*/)
808 {
809   // Terminate analysis
810   //
811   
812   //Save Analysis Settings
813   WriteSettings();
814   
815   // Normalize histos to region area TODO: 
816   // Normalization done at Analysis, taking into account 
817   // area variations on per-event basis (cone case)
818   
819   // Draw some Test plot to the screen
820   if (!gROOT->IsBatch()) {
821     
822     // Update pointers reading them from the output slot
823     fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
824     if( !fListOfHistos ) {
825       AliError("Histogram List is not available");
826       return;
827     }
828     fhEleadingPt         = (TH1F*)fListOfHistos->At(1);
829     fhdNdEtaPhiDist     = (TH1F*)fListOfHistos->At(8);
830     fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11);
831     fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12);
832     fhRegionMultMaxVsEt  = (TH1F*)fListOfHistos->At(14);
833     fhRegionMultMinVsEt  = (TH1F*)fListOfHistos->At(15);
834     fhRegionAveSumPtVsEt = (TH1F*)fListOfHistos->At(16);
835     
836     fhNJets              = (TH1F*)fListOfHistos->At(0);
837     
838     // Canvas
839     TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1200,800);
840     c1->Divide(2,2);
841     c1->cd(1);
842     TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
843     h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1);
844     //h1r->Scale( areafactor );
845     h1r->SetMarkerStyle(20);
846     h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
847     h1r->SetYTitle("P_{T}^{90, max}");
848     h1r->DrawCopy("p");
849     
850     c1->cd(2);
851     
852     TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
853     h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1);
854     //h2r->Scale( areafactor );
855     h2r->SetMarkerStyle(20);
856     h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
857     h2r->SetYTitle("P_{T}^{90, min}");
858     h2r->DrawCopy("p");
859     
860     c1->cd(3);
861     TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
862     h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1);
863     h4r->Scale(2.); // make average
864     //h4r->Scale( areafactor );
865     h4r->SetYTitle("#DeltaP_{T}^{90}");
866     h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
867     h4r->SetMarkerStyle(20);
868     h4r->DrawCopy("p");
869     
870     c1->cd(4);
871     TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading",   "",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
872     TH1F *h6r = new TH1F("hRegionMultMinVsEtleading",   "",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
873     
874     h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1);
875     h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1);
876     //h5r->Scale( areafactor );
877     //h6r->Scale( areafactor );
878     h5r->SetYTitle("N_{Tracks}^{90}");
879     h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
880     h5r->SetMarkerStyle(20);
881     h5r->DrawCopy("p");
882     h6r->SetMarkerStyle(21);
883     h6r->SetMarkerColor(2);
884     h6r->SetYTitle("N_{Tracks}^{90}");
885     h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
886     h6r->DrawCopy("p same");
887     c1->Update();
888     
889     TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800);
890     c2->Divide(2,2);
891     c2->cd(1);
892     fhEleadingPt->SetMarkerStyle(20);
893     fhEleadingPt->SetMarkerColor(2);
894     fhEleadingPt->DrawCopy("p");
895     gPad->SetLogy();
896     
897     c2->cd(2);
898     fhdNdEtaPhiDist->SetMarkerStyle(20);
899     fhdNdEtaPhiDist->SetMarkerColor(2);
900     fhdNdEtaPhiDist->DrawCopy("p");
901     gPad->SetLogy();
902     
903     c2->cd(3);      
904     fhNJets->DrawCopy();
905     
906     //fhTransRegPartPtDist  = (TH1F*)fListOfHistos->At(2);
907     fhRegionMultMin          = (TH1F*)fListOfHistos->At(3);
908     fhMinRegAvePt     = (TH1F*)fListOfHistos->At(4);
909     fhMinRegSumPt     = (TH1F*)fListOfHistos->At(5);   
910     //fhMinRegMaxPtPart     = (TH1F*)fListOfHistos->At(6);
911     fhMinRegSumPtvsMult   = (TH1F*)fListOfHistos->At(7);
912     
913     // Canvas
914     TCanvas* c3 = new TCanvas("c3"," p_{T} dist",160,160,1200,800);
915     c3->Divide(2,2);
916     c3->cd(1);
917     /*fhTransRegPartPtDist->SetMarkerStyle(20);
918      fhTransRegPartPtDist->SetMarkerColor(2); 
919      fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries());
920      fhTransRegPartPtDist->DrawCopy("p");
921      gPad->SetLogy();
922      */
923     c3->cd(2); 
924     fhMinRegSumPt->SetMarkerStyle(20);
925     fhMinRegSumPt->SetMarkerColor(2);  
926     //fhMinRegSumPt->Scale(areafactor);
927     fhMinRegSumPt->DrawCopy("p");
928     gPad->SetLogy();
929     
930     c3->cd(3);
931     fhMinRegAvePt->SetMarkerStyle(20);
932     fhMinRegAvePt->SetMarkerColor(2);  
933     //fhMinRegAvePt->Scale(areafactor);
934     fhMinRegAvePt->DrawCopy("p");
935     gPad->SetLogy();
936     
937     c3->cd(4);
938     
939     TH1F *h7r = new TH1F("hRegionMultMinVsMult",   "",  21, -0.5,   20.5);
940     
941     h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1);
942     
943     h7r->SetMarkerStyle(20);
944     h7r->SetMarkerColor(2);   
945     h7r->DrawCopy("p");
946     
947     c3->Update();
948     
949     
950     /*    c2->cd(3);      
951      fhValidRegion->SetMarkerStyle(20);
952      fhValidRegion->SetMarkerColor(2);
953      fhValidRegion->DrawCopy("p");
954      */    
955     c2->Update();
956   } else {
957     AliInfo(" Batch mode, not histograms will be shown...");
958   }
959   
960   if( fDebug > 1 ) 
961     AliInfo("End analysis");
962   
963 }
964
965 void  AliAnalysisTaskUE::WriteSettings()
966 {
967   
968   fSettingsTree->Branch("fAnaType", &fAnaType, "Ana/I");        
969   fSettingsTree->Branch("fRegionType", &fRegionType,"Reg/I");
970   fSettingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D");
971   fSettingsTree->Branch("fUseChPartJet", &fUseChPartJet,"UseChPart/O");
972   fSettingsTree->Branch("fUseSingleCharge", &fUseSingleCharge,"UseSingleCh/O");
973   fSettingsTree->Branch("fUsePositiveCharge", &fUsePositiveCharge,"UsePositiveCh/O");
974   fSettingsTree->Branch("fOrdering", &fOrdering,"OrderMeth/I");
975   fSettingsTree->Branch("fJet1EtaCut", &fJet1EtaCut, "LeadJetEtaCut/D");
976   fSettingsTree->Branch("fJet2DeltaPhiCut", &fJet2DeltaPhiCut, "DeltaPhi/D");
977   fSettingsTree->Branch("fJet2RatioPtCut", &fJet2RatioPtCut, "Jet2Ratio/D");
978   fSettingsTree->Branch("fJet3PtCut", &fJet3PtCut, "Jet3PtCut/D");
979   fSettingsTree->Branch("fTrackPtCut", &fTrackPtCut, "TrackPtCut/D");
980   fSettingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
981   fSettingsTree->Fill();
982   
983   if (fDebug>5){
984     AliInfo(" All Analysis Settings in Saved Tree");
985     fSettingsTree->Scan();
986   }
987 }