841a4e742d1fb4eff7f94fc8dd526b68db13311e
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliHistogramsUE.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
16 /* $Id:$ */
17 ////////////////////////////////////////////////
18 //--------------------------------------------- 
19 // Class  to handle histograms for UE analysis
20 //---------------------------------------------
21 ////////////////////////////////////////////////
22 #include <TROOT.h>
23 #include <TBranch.h>
24 #include <TCanvas.h>
25 #include <TChain.h>
26 #include <TFile.h>
27 #include <TH1F.h>
28 #include <TH1I.h>
29 #include <TH2F.h>
30 #include <TLegend.h>
31 #include <TList.h>
32 #include <TLorentzVector.h>
33 #include <TMath.h>
34 #include <TProfile.h>
35 #include <TRandom.h>
36 #include <TStyle.h>
37 #include <TSystem.h>
38 #include <TTree.h>
39 #include <TVector3.h>
40
41 #include "AliHistogramsUE.h"
42 #include "AliAnalysisTaskUE.h"
43 #include "AliAnalyseUE.h"
44 #include "AliLog.h"
45
46 ClassImp( AliHistogramsUE)
47
48 //____________________________________________________________________
49 AliHistogramsUE:: AliHistogramsUE():
50 TObject(),
51 fBinsPtInHist(0),
52 fMinJetPtInHist(0.),
53 fMaxJetPtInHist(0.),
54 fTrackEtaCut(0.),
55 fListOfHistos(0x0),  
56 fhNJets(0x0),
57 fhEleadingPt(0x0),
58 fhMinRegPtDist(0x0),
59 fhRegionMultMin(0x0),
60 fhMinRegAvePt(0x0), 
61 fhMinRegSumPt(0x0),            
62 fhMinRegMaxPtPart(0x0),
63 fhMinRegSumPtvsMult(0x0),
64 fhdNdEtaPhiDist(0x0),        
65 fhFullRegPartPtDistVsEt(0x0), 
66 fhTransRegPartPtDistVsEt(0x0),
67 fhRegionSumPtMaxVsEt(0x0),
68 fhRegionMultMax(0x0),         
69 fhRegionMultMaxVsEt(0x0),     
70 fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0),         
71 fhRegionMultMinVsEt(0x0),
72 fhRegionAveSumPtVsEt(0x0),
73 fhRegionDiffSumPtVsEt(0x0),
74 fhRegionAvePartPtMaxVsEt(0x0),
75 fhRegionAvePartPtMinVsEt(0x0),
76 fhRegionMaxPartPtMaxVsEt(0x0),
77 fhRegForwardMult(0x0),
78 fhRegForwardSumPtvsMult(0x0),
79 fhRegBackwardMult(0x0),
80 fhRegBackwardSumPtvsMult(0x0),
81 fhRegForwardPartPtDistVsEt(0x0),
82 fhRegBackwardPartPtDistVsEt(0x0),
83 fhRegTransMult(0x0),
84 fhRegTransSumPtVsMult(0x0),
85 fhMinRegSumPtJetPtBin(0x0),
86 fhMaxRegSumPtJetPtBin(0x0),
87 fhVertexMult(0x0),
88 fh1Xsec(0x0),
89 fh1Trials(0x0)
90 {
91   // Default constructor
92
93 }
94
95 //____________________________________________________________________
96 AliHistogramsUE:: AliHistogramsUE(TList *list):
97 TObject(),
98 fBinsPtInHist(0),
99 fMinJetPtInHist(0.),
100 fMaxJetPtInHist(0.),
101 fTrackEtaCut(0.),
102 fListOfHistos(0x0),  
103 fhNJets(0x0),
104 fhEleadingPt(0x0),
105 fhMinRegPtDist(0x0),
106 fhRegionMultMin(0x0),
107 fhMinRegAvePt(0x0), 
108 fhMinRegSumPt(0x0),            
109 fhMinRegMaxPtPart(0x0),
110 fhMinRegSumPtvsMult(0x0),
111 fhdNdEtaPhiDist(0x0),        
112 fhFullRegPartPtDistVsEt(0x0), 
113 fhTransRegPartPtDistVsEt(0x0),
114 fhRegionSumPtMaxVsEt(0x0),
115 fhRegionMultMax(0x0),         
116 fhRegionMultMaxVsEt(0x0),     
117 fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0),         
118 fhRegionMultMinVsEt(0x0),
119 fhRegionAveSumPtVsEt(0x0),
120 fhRegionDiffSumPtVsEt(0x0),
121 fhRegionAvePartPtMaxVsEt(0x0),
122 fhRegionAvePartPtMinVsEt(0x0),
123 fhRegionMaxPartPtMaxVsEt(0x0),
124 fhRegForwardMult(0x0),
125 fhRegForwardSumPtvsMult(0x0),
126 fhRegBackwardMult(0x0),
127 fhRegBackwardSumPtvsMult(0x0),
128 fhRegForwardPartPtDistVsEt(0x0),
129 fhRegBackwardPartPtDistVsEt(0x0),
130 fhRegTransMult(0x0),
131 fhRegTransSumPtVsMult(0x0),
132 fhMinRegSumPtJetPtBin(0x0),
133 fhMaxRegSumPtJetPtBin(0x0),
134 fhVertexMult(0x0),
135 fh1Xsec(0x0),
136 fh1Trials(0x0)
137 {
138   // Constructor, initialize members from list
139         fhNJets = (TH1F*)list->FindObject("hNJets"); 
140         fhEleadingPt = (TH1F*)list->FindObject("hEleadingPt");
141         fhMinRegPtDist = (TH1F*)list->FindObject("hMinRegPtDist");
142         fhRegionMultMin = (TH1F*)list->FindObject("hRegionMultMin");
143         fhMinRegAvePt = (TH1F*)list->FindObject("hMinRegAvePt");
144         fhMinRegSumPt = (TH1F*)list->FindObject("hMinRegSumPt");
145         fhMinRegMaxPtPart = (TH1F*)list->FindObject("hMinRegMaxPtPart");
146         fhMinRegSumPtvsMult = (TH1F*)list->FindObject("hMinRegSumPtvsMult");
147         fhdNdEtaPhiDist = (TH2F*)list->FindObject("hdNdEtaPhiDist");  
148         fhFullRegPartPtDistVsEt = (TH2F*)list->FindObject("hFullRegPartPtDistVsEt");
149         fhTransRegPartPtDistVsEt = (TH2F*)list->FindObject("hTransRegPartPtDistVsEt");
150         fhRegionSumPtMaxVsEt = (TH1F*)list->FindObject("hRegionSumPtMaxVsEt"); 
151         fhRegionMultMax = (TH1I*)list->FindObject("hRegionMultMax");
152         fhRegionMultMaxVsEt = (TH1F*)list->FindObject("hRegionMultMaxVsEt");
153         fhRegionSumPtMinVsEt = (TH1F*)list->FindObject("hRegionSumPtMinVsEt");
154         fhRegionMultMinVsEt = (TH1F*)list->FindObject("hRegionMultMinVsEt");
155         fhRegionAveSumPtVsEt = (TH1F*)list->FindObject("hRegionAveSumPtVsEt");
156         fhRegionDiffSumPtVsEt = (TH1F*)list->FindObject("hRegionDiffSumPtVsEt");
157         fhRegionAvePartPtMaxVsEt = (TH1F*)list->FindObject("hRegionAvePartPtMaxVsEt");
158         fhRegionAvePartPtMinVsEt = (TH1F*)list->FindObject("hRegionAvePartPtMinVsEt");
159         fhRegionMaxPartPtMaxVsEt = (TH1F*)list->FindObject("hRegionMaxPartPtMaxVsEt");
160         fhRegForwardMult = (TH2F*)list->FindObject("hRegForwardMult"); 
161         fhRegForwardSumPtvsMult = (TH2F*)list->FindObject("hRegForwardSumPtvsMult");
162         fhRegBackwardMult = (TH2F*)list->FindObject("hRegBackwardMult");
163         fhRegBackwardSumPtvsMult = (TH2F*)list->FindObject("hRegBackwardSumPtvsMult");
164         fhRegForwardPartPtDistVsEt = (TH2F*)list->FindObject("hRegForwardPartPtDistVsEt");
165         fhRegBackwardPartPtDistVsEt = (TH2F*)list->FindObject("hRegBackwardPartPtDistVsEt");
166         fhRegTransMult = (TH2F*)list->FindObject("hRegTransMult");
167         fhRegTransSumPtVsMult = (TH2F*)list->FindObject("hRegTransSumPtVsMult");
168         fhMinRegSumPtJetPtBin = (TH2F*)list->FindObject("hMinRegSumPtJetPtBin");
169         fhMaxRegSumPtJetPtBin = (TH2F*)list->FindObject("hMaxRegSumPtJetPtBin");
170         fhVertexMult = (TH1F*)list->FindObject("hVertexMult");
171         fh1Xsec = (TProfile*)list->FindObject("h1Xsec");
172         fh1Trials = (TH1F*)list->FindObject("h1Trials"); 
173
174 }
175 //____________________________________________________________________
176 AliHistogramsUE:: AliHistogramsUE(const AliHistogramsUE & original):
177 TObject(),
178 fBinsPtInHist(original.fBinsPtInHist),
179 fMinJetPtInHist(original.fMinJetPtInHist),
180 fMaxJetPtInHist(original.fMaxJetPtInHist),
181 fTrackEtaCut(original.fTrackEtaCut),
182 fListOfHistos(original.fListOfHistos),  
183 fhNJets(original.fhNJets),
184 fhEleadingPt(original.fhEleadingPt),
185 fhMinRegPtDist(original.fhMinRegPtDist),
186 fhRegionMultMin(original.fhRegionMultMin),
187 fhMinRegAvePt(original.fhMinRegAvePt), 
188 fhMinRegSumPt(original.fhMinRegSumPt),            
189 fhMinRegMaxPtPart(original.fhMinRegMaxPtPart),
190 fhMinRegSumPtvsMult(original.fhMinRegSumPtvsMult),
191 fhdNdEtaPhiDist(original.fhdNdEtaPhiDist),        
192 fhFullRegPartPtDistVsEt(original.fhFullRegPartPtDistVsEt), 
193 fhTransRegPartPtDistVsEt(original.fhTransRegPartPtDistVsEt),
194 fhRegionSumPtMaxVsEt(original.fhRegionSumPtMaxVsEt),
195 fhRegionMultMax(original.fhRegionMultMax),         
196 fhRegionMultMaxVsEt(original.fhRegionMultMaxVsEt),     
197 fhRegionSumPtMinVsEt(original.fhRegionSumPtMinVsEt),         
198 fhRegionMultMinVsEt(original.fhRegionMultMinVsEt),
199 fhRegionAveSumPtVsEt(original.fhRegionAveSumPtVsEt),
200 fhRegionDiffSumPtVsEt(original.fhRegionDiffSumPtVsEt),
201 fhRegionAvePartPtMaxVsEt(original.fhRegionAvePartPtMaxVsEt),
202 fhRegionAvePartPtMinVsEt(original.fhRegionAvePartPtMinVsEt),
203 fhRegionMaxPartPtMaxVsEt(original.fhRegionMaxPartPtMaxVsEt),
204 fhRegForwardMult(original.fhRegForwardMult),
205 fhRegForwardSumPtvsMult(original.fhRegForwardSumPtvsMult),
206 fhRegBackwardMult(original.fhRegBackwardMult),
207 fhRegBackwardSumPtvsMult(original.fhRegBackwardSumPtvsMult),
208 fhRegForwardPartPtDistVsEt(original.fhRegForwardPartPtDistVsEt),
209 fhRegBackwardPartPtDistVsEt(original.fhRegBackwardPartPtDistVsEt),
210 fhRegTransMult(original.fhRegTransMult),
211 fhRegTransSumPtVsMult(original.fhRegTransSumPtVsMult),
212 fhMinRegSumPtJetPtBin(original.fhMinRegSumPtJetPtBin),
213 fhMaxRegSumPtJetPtBin(original.fhMaxRegSumPtJetPtBin),
214 fhVertexMult(original.fhVertexMult),
215 fh1Xsec(original.fh1Xsec),
216 fh1Trials(original.fh1Trials)
217 {
218
219   // Copy constructor
220
221 }
222
223
224 //______________________________________________________________
225 AliHistogramsUE & AliHistogramsUE::operator = (const AliHistogramsUE & /*source*/)
226 {
227
228   // assignment operator
229   return *this;
230
231 }
232
233 //______________________________________________________________
234 TObjArray* AliHistogramsUE::CreateCanvas(const Int_t ncanv){
235         
236         // Create canvas for plotting 
237         printf("Creating %d canvas ... \n",ncanv);
238         TObjArray *arr=new TObjArray;
239         TString name;
240         for(Int_t i=0;i<ncanv;i++ ){
241         name=Form("Canvas %d",i);
242         TCanvas *c=new TCanvas(name,name);
243         c->SetFillColor(0);
244         gStyle->SetOptStat(0);
245         gStyle->SetOptTitle(0);
246         arr->Add(c);
247         }
248 return arr;
249 }
250
251 //____________________________________________________________________
252 void  AliHistogramsUE::CreateHistograms(TList *list,Int_t bins, Double_t min, Double_t max, Double_t etacut)
253 {
254
255   // Create all histograms necessary for UE analysis 
256   fBinsPtInHist = bins;
257   fMinJetPtInHist = min;
258   fMaxJetPtInHist = max;
259   fTrackEtaCut= etacut;
260
261   //Number of reconstructed clusters  
262   fhNJets = new TH1F("hNJets", "Number of clusters",  20, 0, 20);
263   fhNJets->SetXTitle("Number of reconstructed clusters");
264   fhNJets->SetYTitle("#");
265   fhNJets->Sumw2();
266   list->Add( fhNJets );                 // At(0) 
267   
268   //pT distribution of leading clusters
269   fhEleadingPt  = new TH1F("hEleadingPt",   "Leading cluster p_{T}",  bins, min,   max);
270   fhEleadingPt->SetXTitle("p_{T} of  cluster (GeV/c)");
271   fhEleadingPt->SetYTitle("1/N_{ev} dN/dp_{T} (|#eta|<0.5)");
272   fhEleadingPt->Sumw2();
273   list->Add( fhEleadingPt );            // At(1)
274   
275   //Track pT distribution in MIN zone
276   fhMinRegPtDist = new TH1F("hMinRegPtDist",   "p_{T} distribution in MIN zone",  50,0.,20.);
277   fhMinRegPtDist->SetXTitle("Track p_{T} (GeV/c)");
278   fhMinRegPtDist->SetYTitle("dN/dp_{T}");
279   fhMinRegPtDist->Sumw2();
280   list->Add( fhMinRegPtDist );          // At(2)
281   
282   //Multiplicity in MIN zone
283   fhRegionMultMin = new TH1F("hRegionMultMin",      "N_{ch}^{90, min}",  21, -0.5,   20.5);
284   fhRegionMultMin->SetXTitle("N_{ch tracks}");
285   fhRegionMultMin->Sumw2();
286   list->Add( fhRegionMultMin );         // At(3)            
287   
288   //Average pT in MIN region
289   fhMinRegAvePt = new TH1F("hMinRegAvePt", "#LTp_{T}#GT",  50, 0.,   20.);
290   fhMinRegAvePt->SetXTitle("p_{T} (GeV/c)");
291   fhMinRegAvePt->Sumw2();
292   list->Add( fhMinRegAvePt );           // At(4)
293   
294   //Sum pT in MIN region
295   fhMinRegSumPt = new TH1F("hMinRegSumPt", "#Sigma p_{T} ",  50, 0.,   20.);
296   fhMinRegSumPt->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");  
297   fhMinRegSumPt->SetXTitle("#Sigma p_{T} (GeV/c)");
298   fhMinRegSumPt->Sumw2();
299   list->Add( fhMinRegSumPt );           // At(5)
300   
301   //Track with maximum pT in MIN region
302   fhMinRegMaxPtPart = new TH1F("hMinRegMaxPtPart", "max(p_{T})|_{event} ",  50, 0.,   20.);
303   fhMinRegMaxPtPart->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");  
304   fhMinRegMaxPtPart->SetXTitle("p_{T} (GeV/c)");
305   fhMinRegMaxPtPart->Sumw2();
306   list->Add( fhMinRegMaxPtPart );       // At(6)
307   
308   //Sum pT vs. multiplicity in MIN region
309   fhMinRegSumPtvsMult = new TH1F("hMinRegSumPtvsMult", "#Sigma p_{T} vs. Multiplicity ",  21, -0.5,   20.5);
310   fhMinRegSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
311   fhMinRegSumPtvsMult->SetXTitle("N_{charge}");
312   fhMinRegSumPtvsMult->Sumw2();
313   list->Add( fhMinRegSumPtvsMult );     // At(7);
314   
315   //Phi correlation track-cluster vs. leading cluster pT 
316   fhdNdEtaPhiDist  = new TH2F("hdNdEtaPhiDist", Form("Charge particle density |#eta|<%3.1f vs #Delta#phi", fTrackEtaCut),62, 0.,   2.*TMath::Pi(), bins, min, max);
317   fhdNdEtaPhiDist->SetXTitle("#Delta#phi");
318   fhdNdEtaPhiDist->SetYTitle("Leading cluster p_{T}");
319   fhdNdEtaPhiDist->Sumw2();
320   list->Add( fhdNdEtaPhiDist );        // At(8)
321   
322   //Can be used to get track pT distribution for different cluster pT bins (full region)
323   fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", Form( "dN/dp_{T} |#eta|<%3.1f vs Leading cluster p_{T}", fTrackEtaCut),100,0.,50., bins, min, max);
324   fhFullRegPartPtDistVsEt->SetYTitle("Leading cluster p_{T}");
325   fhFullRegPartPtDistVsEt->SetXTitle("p_{T}");
326   fhFullRegPartPtDistVsEt->Sumw2();
327   list->Add( fhFullRegPartPtDistVsEt );  // At(9) 
328   
329   //Can be used to get part pT distribution for different cluster pT bins (transverse region)
330   fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", Form( "dN/dp_{T} in tranvese regions |#eta|<%3.1f vs Leading cluster p_{T}", fTrackEtaCut),100,0.,50., bins, min,   max);
331   fhTransRegPartPtDistVsEt->SetYTitle("Leading cluster p_{T}");
332   fhTransRegPartPtDistVsEt->SetXTitle("p_{T}");
333   fhTransRegPartPtDistVsEt->Sumw2();
334   list->Add( fhTransRegPartPtDistVsEt );  // At(10)
335   
336   //Sum pT in MAX region vs. leading-cluster pT
337   fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt",  "P_{T}^{90, max} vs Leading cluster p_{T}",  bins, min,   max);
338   fhRegionSumPtMaxVsEt->SetXTitle("p_{T} (GeV/c)");
339   fhRegionSumPtMaxVsEt->Sumw2();
340   list->Add( fhRegionSumPtMaxVsEt );      // At(11)
341   
342
343   //Sum pT in MIN region vs. leading-cluster pT
344   fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt",   "P_{T}^{90, min} vs Leading cluster p_{T}",  bins, min,   max);
345   fhRegionSumPtMinVsEt->SetXTitle("p_{T} (GeV/c)");
346   fhRegionSumPtMinVsEt->Sumw2();
347   list->Add( fhRegionSumPtMinVsEt );      // At(12)
348  
349   //Multiplicity in MAX region
350   fhRegionMultMax = new TH1I("hRegionMultMax",      "N_{ch}^{90, max}",  21, -0.5,   20.5);
351   fhRegionMultMax->SetXTitle("N_{ch tracks}");
352   fhRegionMultMax->Sumw2();
353   list->Add( fhRegionMultMax );           // At(13)
354   
355   //Multiplicity in MAX region vs. leading-cluster pT
356   fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt",  "N_{ch}^{90, max} vs Leading cluster p_{T}",  bins, min,   max);
357   fhRegionMultMaxVsEt->SetXTitle("p_{T} (GeV/c)");
358   fhRegionMultMaxVsEt->Sumw2();
359   list->Add( fhRegionMultMaxVsEt );       // At(14)
360   
361   //Multiplicity in MIN region vs. leading-cluster pT
362   fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt",  "N_{ch}^{90, min} vs Leading cluster p_{T}",  bins, min,   max);
363   fhRegionMultMinVsEt->SetXTitle("p_{T} (GeV/c)");
364   fhRegionMultMinVsEt->Sumw2();
365   list->Add( fhRegionMultMinVsEt );      // At(15)
366  
367   //Average sum pT in TRANSVERSE(MIN+MAX) region vs. leading-cluster pT 
368   fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading cluster p_{T}",  bins, min,   max);
369   fhRegionAveSumPtVsEt->SetXTitle("p_{T} (GeV/c)");
370   fhRegionAveSumPtVsEt->Sumw2();
371   list->Add( fhRegionAveSumPtVsEt );     // At(16)
372   
373   //Difference sum pT (MAX-MIN) vs.  leading-cluster pT
374   fhRegionDiffSumPtVsEt= new TH1F("hRegionDiffSumPtVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading cluster p_{T}",  bins, min,   max);
375   fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
376   fhRegionDiffSumPtVsEt->Sumw2();
377   list->Add( fhRegionDiffSumPtVsEt );    // At(17)
378   
379   //Average track pT in MAX region vs. leading-cluster pT
380   fhRegionAvePartPtMaxVsEt = new TH1F("hRegionAvePartPtMaxVsEt", "#LTp_{T}#GT^{90, max} vs Leading cluster p_{T}",  bins, min,   max);
381   fhRegionAvePartPtMaxVsEt->SetXTitle("p_{T} (GeV/c)");
382   fhRegionAvePartPtMaxVsEt->Sumw2();
383   list->Add( fhRegionAvePartPtMaxVsEt );  // At(18)
384   
385   //Average track pT in MIN region vs. leading-cluster pT
386   fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading cluster p_{T}",  bins, min,   max);
387   fhRegionAvePartPtMinVsEt->SetXTitle("p_{T} (GeV/c)");
388   fhRegionAvePartPtMinVsEt->Sumw2();
389   list->Add( fhRegionAvePartPtMinVsEt );   // At(19)
390   
391   //Maximum track pT in MAX region vs. leading-cluster pT
392   fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading cluster p_{T}",  bins, min,   max);
393   fhRegionMaxPartPtMaxVsEt->SetXTitle("p_{T} (GeV/c)");
394   fhRegionMaxPartPtMaxVsEt->Sumw2();
395   list->Add( fhRegionMaxPartPtMaxVsEt );    // At(20)
396   
397   //Multiplicity in FORWARD region
398   fhRegForwardMult = new TH2F("hRegForwardMult", "N_{ch}^{forward}", bins, min, max, 21, -0.5,   20.5);
399   fhRegForwardMult->SetXTitle("N_{ch tracks}");
400   fhRegForwardMult->Sumw2();
401   list->Add( fhRegForwardMult );           // At(25)
402   
403   //Sum pT in FORWARD region vs. multiplicity
404   fhRegForwardSumPtvsMult = new TH2F("hRegForwardSumPtvsMult", "Forward #Sigma p_{T} vs. Multiplicity ", bins, min, max, 21, -0.5,   20.5);
405   fhRegForwardSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
406   fhRegForwardSumPtvsMult->SetXTitle("N_{charge}");
407   fhRegForwardSumPtvsMult->Sumw2();
408   list->Add( fhRegForwardSumPtvsMult );     // At(26);
409   
410   //Multiplicity in BACKWARD region
411   fhRegBackwardMult = new TH2F("hRegBackwardMult", "N_{ch}^{backward}", bins, min, max, 21, -0.5,   20.5);
412   fhRegBackwardMult->SetXTitle("N_{ch tracks}");
413   fhRegBackwardMult->Sumw2();
414   list->Add( fhRegBackwardMult );           // At(27)
415  
416   //Sum pT in BACKWARD region vs. multiplicity
417   fhRegBackwardSumPtvsMult = new TH2F("hRegBackwardSumPtvsMult", "Backward #Sigma p_{T} vs. Multiplicity ", bins, min, max, 21, -0.5,   20.5);
418   fhRegBackwardSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
419   fhRegBackwardSumPtvsMult->SetXTitle("N_{charge}");
420   fhRegBackwardSumPtvsMult->Sumw2();
421   list->Add( fhRegBackwardSumPtvsMult );     // At(28);
422   
423   //Track pT distribution in FORWARD region vs. leading-cluster pT 
424   fhRegForwardPartPtDistVsEt = new TH2F("hRegForwardPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading cluster p_{T}", fTrackEtaCut), 100,0.,50., bins, min, max);
425   fhRegForwardPartPtDistVsEt->SetYTitle("Leading cluster p_{T}");
426   fhRegForwardPartPtDistVsEt->SetXTitle("p_{T} (GeV/c)");
427   fhRegForwardPartPtDistVsEt->Sumw2();
428   list->Add( fhRegForwardPartPtDistVsEt );  // At(29) 
429   
430   //Track pT distribution in BACKWARD region vs. leading-cluster pT 
431   fhRegBackwardPartPtDistVsEt = new TH2F("hRegBackwardPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading cluster p_{T}", fTrackEtaCut), 100,0.,50., bins, min, max);
432   fhRegBackwardPartPtDistVsEt->SetYTitle("Leading cluster p_{T}");
433   fhRegBackwardPartPtDistVsEt->SetXTitle("p_{T}");
434   fhRegBackwardPartPtDistVsEt->Sumw2();
435   list->Add( fhRegBackwardPartPtDistVsEt );  // At(30) 
436   
437   //Multiplicity in TRANSVERSE (MIN+MAX) region
438   fhRegTransMult  = new TH2F("hRegTransMult", "N_{ch}^{transv}", bins, min, max, 21, -0.5,   20.5);
439   fhRegTransMult->SetXTitle("N_{ch tracks}");
440   fhRegTransMult->Sumw2();
441   list->Add( fhRegTransMult );           // At(31)
442   
443   //Sum pT in TRANSVERSE (MIN+MAX) region vs. multiplicity
444   fhRegTransSumPtVsMult = new TH2F("hRegTransSumPtVsMult", "Transverse #Sigma p_{T} vs. Multiplicity ",bins, min, max, 21, -0.5,   20.5);
445   fhRegTransSumPtVsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
446   fhRegTransSumPtVsMult->SetXTitle("N_{charge}");
447   fhRegTransSumPtVsMult->Sumw2();
448   list->Add( fhRegTransSumPtVsMult );     // At(32);
449   
450   //Sum pT in MIN region per cluster pT bin
451   fhMinRegSumPtJetPtBin = new TH2F("hMinRegSumPtJetPtBin", "Transverse Min Reg #Sigma p_{T} per cluster pT bin",bins, min, max, 50, 0.,   20.);
452   fhMinRegSumPtJetPtBin->SetXTitle("Leading cluster p_{T}");
453   fhMinRegSumPtJetPtBin->Sumw2();
454   list->Add( fhMinRegSumPtJetPtBin );           // At(33)
455   
456   //Sum pT in MAX region per cluster pT bin
457   fhMaxRegSumPtJetPtBin = new TH2F("hMaxRegSumPtJetPtBin",      "Transverse Max Reg #Sigma p_{T} per cluster pT bin", bins, min, max, 50, 0.,   20.);
458   fhMaxRegSumPtJetPtBin->SetXTitle("Leading cluster p_{T}");
459   fhMaxRegSumPtJetPtBin->Sumw2();
460   list->Add( fhMaxRegSumPtJetPtBin );           // At(34)
461   
462   //Multiplicity in main vertex
463   fhVertexMult = new TH1F("hVertexMult",      "Multiplicity in Main Vertex", 81, -0.5 , 80.5);
464   fhVertexMult->SetXTitle("Main Vertex Multiplicity");
465   fhVertexMult->Sumw2();
466   list->Add( fhVertexMult ); //At(35)
467   
468   fh1Xsec = new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1); 
469   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
470   fh1Xsec->Sumw2();
471   list->Add( fh1Xsec );            //At(36)
472   
473   fh1Trials = new TH1F("h1Trials","trials from pyxsec.root",1,0,1);
474   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
475   fh1Trials->Sumw2();
476   list->Add( fh1Trials ); //At(37)
477   
478   
479   fListOfHistos=list; 
480 }
481
482
483 //____________________________________________________________________
484 void AliHistogramsUE::DrawUE(Int_t debug){
485
486     // To draw histograms at the end of task running 
487     // Normalize histos to region area TODO: 
488     // Normalization done at Analysis, taking into account 
489     // area variations on per-event basis (cone case)
490    
491     //HIGH WARNING!!!!!: DO NOT SCALE ANY OF THE ORIGINAL HISTOGRAMS
492     //MAKE A COPY, DRAW IT, And later sacale that copy. CAF Issue!!!!!
493     
494     Int_t binsPtInHist = fhEleadingPt->GetNbinsX();
495     Double_t minJetPtInHist = fhEleadingPt->GetXaxis()->GetBinLowEdge(1);
496     Double_t maxJetPtInHist = fhEleadingPt->GetXaxis()->GetBinUpEdge(binsPtInHist);
497      
498     //Sum pT
499     TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1100,700);
500     c1->Divide(2,2);
501     c1->cd(1);
502     TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", binsPtInHist,  minJetPtInHist, maxJetPtInHist);
503     //TH1F *h1r = new TH1F();
504     h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1);
505     //h1r->Scale( areafactor );
506     h1r->SetMarkerStyle(20);
507     h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
508     h1r->SetYTitle("P_{T}^{90, max}");
509     h1r->DrawCopy("p");
510     
511     c1->cd(2);
512     TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", binsPtInHist,  minJetPtInHist, maxJetPtInHist);
513     h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1);
514     //h2r->Scale( areafactor );
515     h2r->SetMarkerStyle(20);
516     h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
517     h2r->SetYTitle("P_{T}^{90, min}");
518     h2r->DrawCopy("p");
519     
520     c1->cd(3);
521     TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", binsPtInHist,  minJetPtInHist, maxJetPtInHist);
522     //TH1F *h41r = new TH1F("hRegForwvsDiffPt" , "", fbinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
523     //TH1F *h42r = new TH1F("hRegBackvsDiffPt" , "", fbinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
524     //h41r->Divide(fhRegForwardSumPtVsEt,fhEleadingPt,1,1);
525     //h42r->Divide(fhRegBackwardSumPtVsEt,fhEleadingPt,1,1);
526     h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1);
527     //h4r->Scale(2.); // make average
528     //h4r->Scale( areafactor );
529     h4r->SetYTitle("#DeltaP_{T}^{90}");
530     h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
531     h4r->SetMarkerStyle(20);
532     h4r->DrawCopy("p");
533
534     c1->cd(4);
535     TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading",   "",  binsPtInHist, minJetPtInHist,   maxJetPtInHist);
536     TH1F *h6r = new TH1F("hRegionMultMinVsEtleading",   "",  binsPtInHist, minJetPtInHist,   maxJetPtInHist);
537     h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1);
538     h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1);
539     //h5r->Scale( areafactor );
540     h5r->SetYTitle("N_{Tracks}^{90}");
541     h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
542     h5r->SetMarkerStyle(20);
543     h5r->DrawCopy("p");
544     h6r->SetMarkerStyle(21);
545     h6r->SetMarkerColor(2);
546     h6r->SetYTitle("N_{Tracks}^{90}");
547     h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
548     h6r->DrawCopy("p same");
549     c1->Update();
550     
551     //Get Normalization
552     //Double_t xsec = fh1Xsec->GetBinContent(1);
553     Double_t xsec = fh1Xsec->GetBinContent(1);
554     Double_t ntrials = fh1Trials->GetBinContent(1);
555     Double_t normFactor = xsec/ntrials;
556     if(debug > 1)Printf("xSec %f nTrials %f Norm %f \n",xsec,ntrials,normFactor);
557     
558     
559     //Jet pT distribution
560     TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800);
561     TH1 * copy = 0;
562     c2->Divide(2,2);
563     c2->cd(1);
564     fhEleadingPt->SetMarkerStyle(20);
565     fhEleadingPt->SetMarkerColor(2);
566     //if( normFactor > 0.) fhEleadingPt->Scale(normFactor);
567     //fhEleadingPt->Draw("p");
568     copy = fhEleadingPt->DrawCopy("p");
569     if( normFactor > 0.) copy->Scale(normFactor);
570     gPad->SetLogy();
571     
572     c2->cd(2);
573     Int_t xbin1 = fhdNdEtaPhiDist->GetYaxis()->FindFixBin(minJetPtInHist);
574     Int_t xbin2 = fhdNdEtaPhiDist->GetYaxis()->FindFixBin(maxJetPtInHist);
575     TH1D* dNdEtaPhiDistAllJets = fhdNdEtaPhiDist->ProjectionX("dNdEtaPhiDistAllJets",xbin1,xbin2);
576     dNdEtaPhiDistAllJets->SetMarkerStyle(20);
577     dNdEtaPhiDistAllJets->SetMarkerColor(2);
578     dNdEtaPhiDistAllJets->DrawCopy("p");
579     gPad->SetLogy();
580     
581     c2->cd(3);      
582     fhNJets->DrawCopy();
583     //c2->cd(4);      
584     //fhValidRegion->DrawCopy("p");
585     //fhTransRegPartPtDist  = (TH1F*)fListOfHistos->At(2);
586     //fhRegionMultMin          = (TH1F*)fListOfHistos->At(3);
587     //fhMinRegAvePt     = (TH1F*)fListOfHistos->At(4);
588     //fhMinRegSumPt     = (TH1F*)fListOfHistos->At(5);   
589     //fhMinRegMaxPtPart     = (TH1F*)fListOfHistos->At(6);
590     //fhMinRegSumPtvsMult   = (TH1F*)fListOfHistos->At(7);
591     c2->Update(); 
592     
593     //pT distributions
594     TCanvas* c3 = new TCanvas("c3"," pT dist",160,160,1200,800);
595     c3->Divide(2,2);
596     c3->cd(1);
597      //fhTransRegPartPtDist->SetMarkerStyle(20);
598      //fhTransRegPartPtDist->SetMarkerColor(2); 
599      //fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries());
600      //fhTransRegPartPtDist->DrawCopy("p");
601      //gPad->SetLogy();
602      
603
604     c3->cd(2); 
605     fhMinRegSumPt->SetMarkerStyle(20);
606     fhMinRegSumPt->SetMarkerColor(2);  
607     //fhMinRegSumPt->Scale(areafactor);
608     fhMinRegSumPt->DrawCopy("p");
609     gPad->SetLogy();
610     
611     c3->cd(3);
612     fhMinRegAvePt->SetMarkerStyle(20);
613     fhMinRegAvePt->SetMarkerColor(2);  
614     //fhMinRegAvePt->Scale(areafactor);
615     fhMinRegAvePt->DrawCopy("p");
616     gPad->SetLogy();
617     
618     c3->cd(4);
619     TH1F *h7r = new TH1F("hRegionMultMinVsMult",   "",  21, -0.5,   20.5);
620     h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1);
621     h7r->SetMarkerStyle(20);
622     h7r->SetMarkerColor(2);   
623     h7r->DrawCopy("p");
624     c3->Update();
625     
626     
627
628     //Save canvas
629     c1->SaveAs("c1.pdf");
630     AliInfo("Canvas 1 saved");
631     c2->SaveAs("c2.pdf");
632     AliInfo("Canvas 2 saved");
633     c3->SaveAs("c3.pdf");
634     AliInfo("Canvas 3 saved");
635   
636 }
637
638 //____________________________________________________________________
639 void AliHistogramsUE::FillHistogram(const char* name, Double_t fillX){
640   
641   // Fill 1D histogram with double
642   ((TH1F*)fListOfHistos->FindObject(name))->Fill(fillX);
643
644 }
645
646 //____________________________________________________________________
647 void AliHistogramsUE::FillHistogram(const char* name, Int_t fillX){
648   
649   // Fill 1D histogram with integer
650   ((TH1F*)fListOfHistos->FindObject(name))->Fill(fillX);
651
652 }
653
654 //____________________________________________________________________
655 void AliHistogramsUE::FillHistogram(const char* name, Double_t fillX, Double_t fillY){
656   
657  // Case of TH1F with weight or TH2F w/o weight
658  TObject *obj = fListOfHistos->FindObject(name);
659  if (obj->InheritsFrom("TH1F")){
660         ((TH1F*)fListOfHistos->FindObject(name))->Fill(fillX, fillY);
661  } else {
662         ((TH2F*)fListOfHistos->FindObject(name))->Fill(fillX, fillY);
663         }
664
665
666 }
667
668 //____________________________________________________________________
669 void AliHistogramsUE::FillHistogram(const char* name, Double_t fillX, Double_t fillY, Double_t weight){
670
671   // Fill 2D histogram with double and weight   
672   ((TH2F*)fListOfHistos->FindObject(name))->Fill(fillX, fillY, weight);
673
674 }
675
676 //____________________________________________________________________
677 void AliHistogramsUE::FillHistogram(const char* name, Double_t fillX, Int_t fillY, Double_t weight){
678   
679   // Fill 2D histogram with integer and weight   
680   ((TH2F*)fListOfHistos->FindObject(name))->Fill(fillX, fillY, weight);
681
682 }
683
684 //____________________________________________________________________
685 TObjArray*  AliHistogramsUE::GetHistosForPlotting(TString data, TString branches){
686
687         // Instance filled histos for plotting purpose
688         printf("Creating histograms ... \n");
689                 
690         printf("Reading file: %s\n",data.Data());
691         
692         // Read input files ----------------------------------------- 
693         TFile *fdata = new TFile(data.Data());
694         TDirectoryFile *ddata[20];
695         TList *ldata[20];
696
697         TObjArray  *arrb=branches.Tokenize(";");
698         TIter next(arrb);
699         TObject *o=0;
700         Int_t br=0;
701         while ( (o=next()) ){
702                 ddata[br] = (TDirectoryFile*)fdata->Get(Form("PWG4_UE%s",o->GetName()));
703                 if(!ddata[br]) printf("ERROR: No histo dir found! \n");
704                 ldata[br] = (TList*)ddata[br]->Get(Form("histosUE%s",o->GetName()));
705                 printf("Reading branch: %s\n",o->GetName());
706                 if(!ldata[br]) printf("ERROR: No histo list found! \n");
707                 br++;
708         }
709         
710         TObjArray *arr=new TObjArray;
711
712         TH1F *hjets[20];                // accepted leading jets
713         TH1F *hnjets[20];          // number of accepted jets
714         TH2F *hetaphi[20];         // delta-phi particle-jet correlation
715         TH2F *hptfull[20];         // particle pT all regions vs. jet pT 
716         TH2F *hpttransv[20];       // particle pT transv. regions vs. jet pT 
717         TH1F *hmax[20];         // sum pT in MAX region
718         TH1F *hmin[20];         // sum pT in MIN region
719         TH1F *hmultmax[20];     // multiplicity in MAX region
720         TH1F *hmultmin[20];       // multiplicity in MIN region
721
722         for (Int_t i =0; i<br; i++){
723
724                 //Number of jets --------------------------------------------
725                 hnjets[i] = (TH1F*) ldata[i]->FindObject("fhNJets"); 
726                 hnjets[i]->GetXaxis()->SetTitle("Number of jets");
727                 hnjets[i]->GetYaxis()->SetTitle("1/n_{ev} dN/dN_{jets}");
728                 hnjets[i]->SetMarkerStyle(20); 
729                 hnjets[i]->SetMarkerColor(1+i);
730
731
732                 //Leading jets ----------------------------------------------
733                 hjets[i] = (TH1F*) ldata[i]->FindObject("hEleadingPt"); 
734                 hjets[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
735                 hjets[i]->GetYaxis()->SetTitle("1/n_{ev} dN/dp_{T} (|#eta<0.5|)");
736                 hjets[i]->SetMarkerStyle(20); 
737                 hjets[i]->SetMarkerColor(1+i);
738                 hjets[i]->SetMinimum(0.1);
739                 hjets[i]->SetMaximum(1000.);
740
741
742                 //Transverse Region MAX -------------------------------------
743                 hmax[i] = (TH1F*) ldata[i]->FindObject("hRegionSumPtMaxVsEt");
744                 if (!hmax[i])AliInfo("Histo not found!!!");
745                 hmax[i]->GetXaxis()->SetTitle("Leading jet P_{T} (GeV/c)");
746                 hmax[i]->GetYaxis()->SetTitle("P_{T}^{90,max} (GeV/c)");
747                 hmax[i]->SetMarkerStyle(20);
748                 hmax[i]->SetMarkerColor(1+i);
749                 hmax[i]->Divide(hjets[i]); // normalize for jet spectrum
750                 hmax[i]->SetMaximum(5.);
751
752
753                 //Transverse Region MIN -------------------------------------
754                 hmin[i] = (TH1F*) ldata[i]->FindObject("hRegionSumPtMinVsEt");
755                 hmin[i]->GetXaxis()->SetTitle("Leading jet P_{T} (GeV/c)");
756                 hmin[i]->GetYaxis()->SetTitle("P_{T}^{90,min} (GeV/c)");
757                 hmin[i]->SetMarkerStyle(20);
758                 hmin[i]->SetMarkerColor(1+i);
759                 hmin[i]->SetMaximum(3.);
760                 hmin[i]->Divide(hjets[i]); // normalize for jet spectrum
761
762
763                 //Multiplicity MAX ------------------------------------------
764                 hmultmax[i] = (TH1F*) ldata[i]->FindObject("hRegionMultMaxVsEt");
765                 hmultmax[i]->GetXaxis()->SetTitle("Leading Jet P_{T} (GeV/c)");
766                 hmultmax[i]->GetYaxis()->SetTitle("N_{ch}^{90,max}");
767                 hmultmax[i]->SetMarkerStyle(20);
768                 hmultmax[i]->SetMarkerColor(1+i);
769                 hmultmax[i]->SetMaximum(10.);
770                 hmultmax[i]->Divide(hjets[i]); // normalize for jet spectrum
771
772
773                 //Multiplicity MIN ------------------------------------------
774                 hmultmin[i] = (TH1F*) ldata[i]->FindObject("hRegionMultMinVsEt");
775                 hmultmin[i]->GetXaxis()->SetTitle("Leading Jet P_{T} (GeV/c)");
776                 hmultmin[i]->GetYaxis()->SetTitle("N_{ch}^{90,min}");
777                 hmultmin[i]->SetMarkerStyle(20);
778                 hmultmin[i]->SetMarkerColor(1+i);
779                 hmultmin[i]->SetMaximum(3.);
780                 hmultmin[i]->Divide(hjets[i]); // normalize for jet spectrum
781             
782
783                 // Phi-correlation with leading jet --------------------------
784                 hetaphi[i] = (TH2F*) ldata[i]->FindObject("hdNdEtaPhiDist");
785                 hetaphi[i]->GetXaxis()->SetTitle("#Delta #phi (w.r.t. leading jet)");
786                 hetaphi[i]->SetMarkerStyle(20);
787             
788
789                 // pT distribution in full region vs. jet pT --------------------------
790                 hptfull[i] = (TH2F*) ldata[i]->FindObject("hFullRegPartPtDistVsEt");
791                 hptfull[i]->GetYaxis()->SetTitle("Leading Jet P_{T} (GeV/c)");
792                 hptfull[i]->GetXaxis()->SetTitle("Track P_{T} (GeV/c)");
793             
794
795                 // pT distribution in transv region vs. jet pT --------------------------
796                 hpttransv[i] = (TH2F*) ldata[i]->FindObject("hTransRegPartPtDistVsEt");
797                 hpttransv[i]->GetYaxis()->SetTitle("Leading Jet P_{T} (GeV/c)");
798                 hpttransv[i]->GetXaxis()->SetTitle("Track P_{T} (GeV/c)");
799             
800
801                 // Return Histos
802                 arr->Add(hnjets[i]);      //at 0 number of jets
803                 arr->Add(hjets[i]);     //at 1 leading jets
804                 arr->Add(hmax[i]);      //at 2 sum pT MAX
805                 arr->Add(hmin[i]);      //at 3 sumpT MIN
806                 arr->Add(hmultmax[i]);   //at 4 multiplicity MAX 
807                 arr->Add(hmultmin[i]);   //at 5 multiplicity MIN
808                 arr->Add(hetaphi[i]);     //at 6 phi correlation
809                 arr->Add(hptfull[i]);     //at 7 pT distr in full region
810                 arr->Add(hpttransv[i]);   //at 8 pT distr in transv region
811
812         }
813         return arr;
814 }
815
816 //_______________________________________________________________________________________
817 void AliHistogramsUE::SetStyle(){
818
819   // Set plotting style
820   gPad->SetFrameFillColor(0);
821   gPad->SetFillColor(0);
822   gPad->SetBorderSize(2);
823   gPad->SetGridy();
824   gPad->SetFrameBorderMode(0);
825   //gStyle->SetOptStat(0);
826   gStyle->SetOptTitle(0);
827
828 }
829
830
831 //____________________________________________________________________
832 TList* AliHistogramsUE::GetHistograms(){
833
834   // Return list of relevant histograms 
835   return fListOfHistos;
836
837 }
838
839 //____________________________________________________________________
840 void AliHistogramsUE::PlotBranchesUE(TString file, TString branches, Double_t minJetProjection){
841
842   // Function to be called by external macro to plot analysis from different jet branches 
843   
844   //Count the branches 
845   TObjArray  *arr=branches.Tokenize(";");
846   TIter next(arr);
847   TObject *o=0;
848   Int_t br=0;
849   while ( (o=next()) ){
850         br++;
851         }
852   
853
854   //Canvas
855   TObjArray *arrC=CreateCanvas(9);
856  
857   //Histograms 
858   TObjArray *arrH=GetHistosForPlotting(file.Data(),branches.Data());
859   TH1F *hnjets[20];
860   TH1F *hjets[20];
861   TH1F *hmax[20];
862   TH1F *hmin[20];
863   TH1F *hmultmax[20];
864   TH1F *hmultmin[20];
865   TH2F *hetaphi[20];
866   TH2F *hptfull[20];
867   TH2F *hpttransv[20];
868
869   for (Int_t i= 0; i<br; i++){
870         hnjets[i]=((TH1F*)arrH->At(9*i));
871         hjets[i]=((TH1F*)arrH->At(1+(9*i)));
872         hmax[i]=((TH1F*)arrH->At(2+(9*i)));
873         hmin[i]=((TH1F*)arrH->At(3+(9*i)));
874         hmultmax[i]=((TH1F*)arrH->At(4+(9*i)));
875         hmultmin[i]=((TH1F*)arrH->At(5+(9*i)));
876         hetaphi[i]=((TH2F*)arrH->At(6+(9*i))); 
877         hptfull[i]=((TH2F*)arrH->At(7+(9*i))); 
878         hpttransv[i]=((TH2F*)arrH->At(8+(9*i))); 
879         }
880
881   //Define jet-pT range in projections
882   Int_t binstartprojection = hetaphi[0]->GetYaxis()->FindFixBin(minJetProjection); //be careful...  
883   Int_t binstopprojection = hetaphi[0]->GetYaxis()->GetNbins();
884
885   Double_t normphi[20]; 
886
887   for (Int_t i= 0; i<br; i++){
888         normphi[i]=hjets[i]->Integral(binstartprojection,binstopprojection);    
889         hnjets[i]->Scale(1./(hnjets[i]->GetBinWidth(1)*hnjets[i]->GetEntries()));
890         hjets[i]->Scale(1./(hjets[i]->GetBinWidth(1)*hnjets[i]->GetEntries()));
891         }
892
893   //LEGENDS
894   //Legend jets
895   TLegend *leg=new TLegend(0.5,0.6,0.89,0.89);  // for jet spectrum
896   leg->SetFillColor(0);
897   leg->SetHeader("Jet Finders:");
898   //Legend density
899   TLegend *legd = new TLegend(0.1077586,0.6016949,0.4971264,0.8919492,NULL,"brNDC");
900   legd->SetFillColor(0);
901   legd->SetHeader("Jet Finders:");
902   //Legend pT distributions
903   TLegend *legp = new TLegend(0.1364943,0.1292373,0.5258621,0.4194915,NULL,"brNDC");
904   legp->SetFillColor(0);
905   legp->SetHeader("Jet Finders:");
906
907   arr=branches.Tokenize(";");
908   TIter next1(arr);
909   o=0;
910   Int_t brleg=0;
911   while ( (o=next1()) ){
912         leg->AddEntry(hjets[brleg],Form("UE%s",o->GetName()),"p");
913         legd->AddEntry(hjets[brleg],Form("UE%s",o->GetName()),"p");
914         legp->AddEntry(hjets[brleg],Form("UE%s",o->GetName()),"p");
915         brleg++;
916         }
917
918   //1) NUMBER OF CLUSTERS 
919   TCanvas *c0=((TCanvas*)arrC->At(0)); 
920   c0->cd();
921   SetStyle();
922   gPad->SetLogy();
923   hnjets[0]->SetMaximum(1.4);
924   hnjets[0]->Draw("E1");
925   for (Int_t i= 1; i<br; i++){
926         hnjets[i]->Draw("same");
927         }
928   leg->Draw("same");
929
930   //2) LEADING CLUSTERS pT 
931   TCanvas *c1=((TCanvas*)arrC->At(1)); 
932   c1->cd();
933   SetStyle();
934   gPad->SetLogy();
935   hjets[0]->Draw("E1");
936   hjets[0]->SetMaximum(1.4);
937   for (Int_t i= 1; i<br; i++){
938         hjets[i]->Draw("same");
939         }
940   leg->Draw("same");
941
942   //3) SUM-pT IN MAX REGION
943   TCanvas *c2=((TCanvas*)arrC->At(2)); 
944   c2->cd();
945   SetStyle();
946   hmax[0]->Draw("E1");
947   for (Int_t i= 1; i<br; i++){
948         hmax[i]->Draw("same");
949         }
950   leg->Draw("same");
951
952   //4) SUM-pT IN MIN REGION
953   TCanvas *c3=((TCanvas*)arrC->At(3)); 
954   c3->cd();
955   SetStyle();
956   hmin[0]->GetYaxis()->SetRangeUser(0.,2.5);
957   hmin[0]->Draw("E1");
958   for (Int_t i= 1; i<br; i++){
959         hmin[i]->Draw("same");
960         }
961   leg->Draw("same");
962
963   //5) MULTIPLICITY IN MAX REGION
964   TCanvas *c4=((TCanvas*)arrC->At(4));
965   c4->cd();
966   SetStyle();
967   hmultmax[0]->GetYaxis()->SetRangeUser(0.,5.);
968   hmultmax[0]->Draw("E1");
969   for (Int_t i= 1; i<br; i++){
970         hmultmax[i]->Draw("same");
971         }
972   leg->Draw("same");
973
974   //6) MULTIPLICITY IN MIN REGION
975   TCanvas *c5=((TCanvas*)arrC->At(5));
976   c5->cd();
977   SetStyle();
978   hmultmin[0]->GetYaxis()->SetRangeUser(0.,2.5);
979   hmultmin[0]->Draw("E1");
980   for (Int_t i= 1; i<br; i++){
981         hmultmin[i]->Draw("same");
982         }
983   leg->Draw("same");
984
985   //7) JET-TRACK CORRELATION
986   TCanvas *c6=((TCanvas*)arrC->At(6));
987   c6->cd();
988   SetStyle();
989   gPad->SetLogy();
990   TH1D *tmpetaphi[20];
991   for (Int_t i= 0; i<br; i++){
992         tmpetaphi[i]=new TH1D();
993         tmpetaphi[i]=hetaphi[i]->ProjectionX(Form("data%d",i),binstartprojection);
994         tmpetaphi[i]->GetYaxis()->SetTitle("1/n_{lj} dN/d#Delta #phi");
995         tmpetaphi[i]->Scale(1./(hetaphi[i]->GetBinWidth(1)*normphi[i]));
996         tmpetaphi[i]->SetMarkerColor(i+1);
997         tmpetaphi[i]->GetXaxis()->SetLimits(-3.*TMath::Pi()/2.,TMath::Pi()/2.);
998         tmpetaphi[i]->GetYaxis()->SetRangeUser(0.5,20.);
999         if (i==0) tmpetaphi[i]->Draw("E1");
1000         else tmpetaphi[i]->Draw("same");
1001         // evaluate mean multiplicity in transverse regions
1002         Double_t err1=0.;
1003         Double_t err2=0.;
1004         Double_t err3=0.;
1005         Double_t mean=(tmpetaphi[i]->IntegralAndError(1,5,err1)+tmpetaphi[i]->IntegralAndError(27,36,err2)+tmpetaphi[i]->IntegralAndError(58,62,err3))/(20.);
1006         err1=TMath::Sqrt(err1*err1+err2*err2+err3*err3)/20.;
1007         Printf("Branch: %d  MeanTransvMult: %f err: %f",i,mean,err1);
1008         }
1009   legd->Draw("same");
1010
1011   //8) TRACK-pT DISTRIBUTION IN FULL REGION
1012   TCanvas *c7=((TCanvas*)arrC->At(7));
1013   c7->cd();
1014   SetStyle();
1015   gPad->SetLogy();
1016   gPad->SetLogx();
1017   gPad->SetTitle("Full region (projection X)");
1018   TH1D *tmpfull[20];
1019   for (Int_t i= 0; i<br; i++){
1020         tmpfull[i]=new TH1D();
1021         tmpfull[i]=hptfull[i]->ProjectionX(Form("full%d",i),binstartprojection);
1022         tmpfull[i]->GetYaxis()->SetTitle("1/n_{lj} dN/dp_{T}");
1023         tmpfull[i]->Scale(1./(tmpfull[i]->GetBinWidth(1)*normphi[i]));
1024         tmpfull[i]->SetMarkerStyle(20);
1025         tmpfull[i]->SetMarkerColor(i+1);
1026         tmpfull[i]->SetMaximum(100.);
1027         if (i==0)tmpfull[i]->Draw("E1");
1028         else tmpfull[i]->Draw("same"); 
1029         }
1030   legp->Draw("same");
1031
1032   //9) TRACK-pT DISTRIBUTION IN TRANSVERSE (MIN+MAX) REGION
1033   TCanvas *c8=((TCanvas*)arrC->At(8));
1034   c8->cd();
1035   SetStyle();
1036   gPad->SetLogy();
1037   gPad->SetLogx();
1038   gPad->SetTitle("Transverse regions (projection X)");
1039   TH1D *tmptransv[20];
1040   for (Int_t i= 0; i<br; i++){
1041         tmptransv[i]=new TH1D();
1042         tmptransv[i]=hpttransv[i]->ProjectionX(Form("transv%d",i),binstartprojection);
1043         tmptransv[i]->GetYaxis()->SetTitle("1/n_{lj} dN/dp_{T}");
1044         tmptransv[i]->Scale(1./(tmptransv[i]->GetBinWidth(1)*normphi[i]));
1045         tmptransv[i]->SetMarkerStyle(20);
1046         tmptransv[i]->SetMarkerColor(i+1);
1047         tmptransv[i]->SetMaximum(10.);
1048         if (i==0)tmptransv[i]->Draw("E1");
1049         else tmptransv[i]->Draw("same"); 
1050         }
1051   legp->Draw("same");
1052 }