added some control histos, adapted the eta range and removed dependcies on generated...
[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 fhdNdEta_PhiDist(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     fhdNdEta_PhiDist->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 Pt   = jet->E()+track1->Pt();  // Scalar sum of Pt
574         // recalculating the centroid
575         Double_t eta = jet->Eta()*jet->E()/Pt + track1->Eta()/Pt;
576         Double_t phi = jet->Phi()*jet->E()/Pt + track1->Phi()/Pt;
577         jet->SetPtEtaPhiE( 1., eta, phi, Pt );
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   fhdNdEta_PhiDist  = new TH1F("hdNdEta_PhiDist",   "Charge particle density |#eta|< 1 vs #Delta#phi",  120, 0.,   2.*TMath::Pi());
721   fhdNdEta_PhiDist->SetXTitle("#Delta#phi");
722   fhdNdEta_PhiDist->SetYTitle("dN_{ch}/d#etad#phi");
723   fhdNdEta_PhiDist->Sumw2();
724   fListOfHistos->Add( fhdNdEta_PhiDist );        // 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     fhdNdEta_PhiDist     = (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     fhdNdEta_PhiDist->SetMarkerStyle(20);
899     fhdNdEta_PhiDist->SetMarkerColor(2);
900     fhdNdEta_PhiDist->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 }