moved some TLists to members of class
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskPartonDisc.cxx
1 /************************************************************************* 
2 * Copyright(c) 1998-2008, 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 ////////////////////////////////////////////////////////
17 //                                                    //
18 // Analysis task for parton discrimination studies    //
19 //                                                    // 
20 // Author:                                            //
21 // Hermes Leon Vargas  (hleon@ikf.uni-frankfurt.de)   //
22 //                                                    // 
23 ////////////////////////////////////////////////////////
24
25 #include "TROOT.h"
26 #include "TDirectory.h"
27 #include "TKey.h"
28 #include "TList.h"
29 #include "TSystem.h"
30 #include "TChain.h"
31 #include "TTree.h"
32 #include "TFile.h"
33 #include "TString.h"
34 #include "TH1F.h"
35 #include "TH1I.h"
36 #include "TH2F.h"
37 #include "TH2I.h"
38 #include "TH3F.h"
39 #include "TProfile.h"
40 #include "TCanvas.h"
41 #include "AliAnalysisTask.h"
42 #include "AliAnalysisManager.h"
43 #include "AliAODHandler.h"
44 #include "AliVEvent.h"
45 #include "AliAODEvent.h"
46 #include "AliAODInputHandler.h"
47 #include "AliVParticle.h"
48 #include "AliAODMCParticle.h"
49 #include "AliGenEventHeader.h"
50 #include "AliGenPythiaEventHeader.h"
51 #include "AliGenDPMjetEventHeader.h"
52 #include "AliAODMCHeader.h"
53 #include "AliMCEventHandler.h"
54 #include "AliMCEvent.h"
55 #include <algorithm> 
56 #include "AliLog.h"
57 #include "AliAODVertex.h"
58 #include "AliAODHeader.h"
59 #include "AliESDtrack.h"
60 #include "AliAnalysisHelperJetTasks.h"
61 #include "AliAnalysisTaskPartonDisc.h"
62 #include "AliAODVZERO.h"
63 #include "TRandom3.h"
64 #include "TDatime.h"
65
66 // Analysis task for parton discrimination
67
68 ClassImp(AliAnalysisTaskPartonDisc)
69
70 Double_t *AliAnalysisTaskPartonDisc::fgContainer = 0x0; 
71
72 //________________________________________________________________________
73 AliAnalysisTaskPartonDisc::AliAnalysisTaskPartonDisc() 
74   : AliAnalysisTaskSE(), fAOD(0), fUseAODMC(kFALSE), fPhojetMC(kFALSE), fBranchMC("jetsMC"), fBranchRec("jetsREC"), fBranchSecRec(""), fSqrts(0),  fNtX(0), fJetRadius(0.), fFlavorRadius(0.), fFilterBit(0xFF), fOutputList(0), fJetPt(0), fJetPtSec(0), fJetPtMC(0), fJetEta(0), fJetEtaSec(0), fJetPhi(0), fJetPhiSec(0), fJetEtaMC(0), fJetPhiMC(0), fPtAODMC(0), fPtAOD(0), fEtaAODMC(0), fPhiAODMC(0), fEtaAOD(0), fPhiAOD(0), fFlavor(0), fNJetsMC(0), fNJetsRD(0), fNJetsRDSeco(0), fJetsMultPtMC(0), fJetsMultPtRD(0), fNChTrRD(0), fProfNChTrRD(0), fFracQQ(0), fFracGQ(0), fFracGG(0), fFracOutGoingQQ(0), fFracOutGoingGQ(0), fFracOutGoingGG(0), fh1Xsec(0), fh1Trials(0), fMpdg(0), fProcessJetPt(0), fFlavorLead(0), fProcessLeadJetPt(0), fPDGMothLPart(0), fFlavProc(0), fAvgTrials(1), fUseAODJetInput(kFALSE), fMinTrackPtInNTX(0), fMaxTrackPtInNTX(0), fSCMRD(0), fMinpTVal(0), fZVertex(0), fh1Events(0), fUseOnlyMC(kFALSE), fCheckMCStatus(kTRUE), fEvtCount(0), fNAccJetsMC(0), fNAccJetsRD(0), fNAccJetsRDSeco(0), fEnablePrints(kFALSE), fRecJetPtInclusive(0), fMCJetPtInclusive(0), fRecJetPtLeading(0), fMCJetPtLeading(0), fSecRecJetPtInclusive(0), fSecRecJetPtLeading(0), fHasPerpCone(kTRUE), fEtaPerpCoord(0), fPhiPerpCoord(0), fPtPerpCoord(0), fJetEvent(kFALSE), fPerpCone(0), fNChTrMCPerp(0), fNChTrRecPerp(0), fSCMMCPerp(0), fSCMRecPerp(0), fIsHIevent(kFALSE), fCurrentJetMinPtNT90(0), fBckgSbsJet(0), fCurrentJetMinPtNT90Recalc(0), fNChTrCorrMCQuark(0), fNChTrCorrMCGluon(0), fNChTrCorrMCPerp(0), fIsPossibleToSubstBckg(kTRUE), fNChTrRecECorr(0), fNChTrRecPerpECorr(0), fRefMult(0), fCurrentJetCharge(0), fRefMultWOJet(0), fVZEROMult(0), fMultWOJetVZero(0), fVZero(0), fRefMultFullV0(0), fRefMultV0Corr(0), fFullV0V0Corr(0), fNTXV0MultPt(0), fNTXCBMultPt(0), fMinpTValUE(2.0), fRefMultFullV0UJ(0), fRefMultV0CorrUJ(0), fFullV0V0CorrUJ(0), fMultWOJetVZeroUJ(0), fRefMultWOJetUJ(0), fMaxpTValUE(2.0), fRefAODTrackCount(0), fRefAODTrackCountUJ(0), fTrackCountWOJet(0), fTrackCountWOJetUJ(0), fTrackCountWOJetUJMC(0), fFullV0V0CorrUJMC(0), fMinpTValMC(2.0), fIncExcR(0.0), fForceNotTR(kFALSE), fNotExtDiJEx(kFALSE), fMinTrackPtInNTXRecalc(0), fMaxTrackPtInNTXRecalc(0), fPtDistInJetConeRaw(0), fPtDistInPerpConeRaw(0), fPtInPerpCon(0), fMinTrackPtInNTXR(0), fMaxTrackPtInNTXR(0), fEventCent(0), fJetEtaAll(0), fJetEtaOnlyTPCcut(0), fNChTrRecECorrPPMult(0), fNChTrRecPerpECorrPPMult(0), fForceSkipSJ(kFALSE), fJetPtCentPbPbRaw(0), fJetPtCentPbPbCorr(0), fJetAcceptance(0.5), fIncreasingExcl(kFALSE), fTotTracksCone(0), fTotTracksInCone(0), fTTrackRandomRejection(0), fJTrackRandomRejection(0), fMinPtInGlobMult(0)
75 {
76   // Constructor
77
78   for(Int_t a=0; a<16;a++)
79     {
80       fJetFlags[a]=kTRUE;
81       if(a<12)
82         {
83           fNChTr[a]=0;
84           fHistPtParton[a]=0;
85           fSCM[a]=0;
86           if(a<8)
87             {
88               fNChTrRDMult[a]=0;
89               fNAccJetsRDMult[a]=0;
90               fTotalJetCharge[a]=0;
91               fSCMRDMult[a]=0;
92               fNChTrRDMultMC[a]=0;
93               fSCMRDMultMC[a]=0;
94               fNChTrRDMultSE[a]=0;
95               fNAccJetsRDMultSE[a]=0;
96               fTotalJetChargeSE[a]=0;
97               fSCMRDMultSE[a]=0;
98               fNChTrRDMultOJ[a]=0;
99               fSCMRDMultOJ[a]=0;
100               fNChTrRDMultSEOJ[a]=0;
101               fSCMRDMultSEOJ[a]=0;
102               fNChTrRDMultOJMC[a]=0;
103               fSCMRDMultOJMC[a]=0;
104               fNChTrRDMultSEOJMC[a]=0;
105               fSCMRDMultSEOJMC[a]=0;
106               fNChTrRecPerpMultSEOJ[a]=0;
107               fJEtaMCMultOJ[a]=0;  
108               fJEtaMCMultSEOJ[a]=0;
109               fJEtaRDMultOJ[a]=0;           
110               fJEtaRDMultSEOJ[a]=0;         
111               fJetPtMCMultOJ[a]=0;         
112               fJetPtMCMultSEOJ[a]=0;       
113               fJetPtRDMultOJ[a]=0;        
114               fJetPtRDMultSEOJ[a]=0;       
115               fEntriesQuark[a]=0;          
116               fEntriesGluon[a]=0;          
117             }
118           if(a<6)
119             {
120               fProcessPDG[a]=0;
121               fFragPion[a]=0;
122               fFragKaon[a]=0;
123               fFragProton[a]=0;
124               fFragChargedR4[a]=0;
125               fFragChargedR3[a]=0;
126               fFragChargedR2[a]=0;
127               fHistContainerR4[a]=0;
128               fHistContainerR3[a]=0;
129               fHistContainerR2[a]=0;
130               if(a<3)
131                 {
132                   fJetEtaJetPt[a]=0;
133                   if(a<2)
134                     {
135                       fFragCandidates[a]=0;
136                       fMinTrackPtInNTXh[a]=0;
137                       fMaxTrackPtInNTXh[a]=0; 
138                     }
139                 }
140             }
141         }
142     }
143 }
144 //________________________________________________________________________
145 AliAnalysisTaskPartonDisc::AliAnalysisTaskPartonDisc(const char *name) 
146   : AliAnalysisTaskSE(name), fAOD(0), fUseAODMC(kFALSE), fPhojetMC(kFALSE), fBranchMC("jetsMC"), fBranchRec("jetsREC"), fBranchSecRec(""), fSqrts(0),  fNtX(0), fJetRadius(0.), fFlavorRadius(0.), fFilterBit(0xFF), fOutputList(0), fJetPt(0), fJetPtSec(0), fJetPtMC(0), fJetEta(0), fJetEtaSec(0), fJetPhi(0), fJetPhiSec(0), fJetEtaMC(0), fJetPhiMC(0), fPtAODMC(0), fPtAOD(0), fEtaAODMC(0), fPhiAODMC(0), fEtaAOD(0), fPhiAOD(0), fFlavor(0), fNJetsMC(0), fNJetsRD(0), fNJetsRDSeco(0), fJetsMultPtMC(0), fJetsMultPtRD(0), fNChTrRD(0), fProfNChTrRD(0), fFracQQ(0), fFracGQ(0), fFracGG(0), fFracOutGoingQQ(0), fFracOutGoingGQ(0), fFracOutGoingGG(0), fh1Xsec(0), fh1Trials(0), fMpdg(0), fProcessJetPt(0), fFlavorLead(0), fProcessLeadJetPt(0), fPDGMothLPart(0), fFlavProc(0), fAvgTrials(1), fUseAODJetInput(kFALSE), fMinTrackPtInNTX(0), fMaxTrackPtInNTX(0), fSCMRD(0), fMinpTVal(0), fZVertex(0), fh1Events(0), fUseOnlyMC(kFALSE), fCheckMCStatus(kTRUE), fEvtCount(0), fNAccJetsMC(0), fNAccJetsRD(0), fNAccJetsRDSeco(0), fEnablePrints(kFALSE), fRecJetPtInclusive(0), fMCJetPtInclusive(0), fRecJetPtLeading(0), fMCJetPtLeading(0), fSecRecJetPtInclusive(0), fSecRecJetPtLeading(0), fHasPerpCone(kTRUE), fEtaPerpCoord(0), fPhiPerpCoord(0), fPtPerpCoord(0), fJetEvent(kFALSE), fPerpCone(0), fNChTrMCPerp(0), fNChTrRecPerp(0), fSCMMCPerp(0), fSCMRecPerp(0), fIsHIevent(kFALSE), fCurrentJetMinPtNT90(0), fBckgSbsJet(0), fCurrentJetMinPtNT90Recalc(0), fNChTrCorrMCQuark(0), fNChTrCorrMCGluon(0), fNChTrCorrMCPerp(0), fIsPossibleToSubstBckg(kTRUE), fNChTrRecECorr(0), fNChTrRecPerpECorr(0), fRefMult(0), fCurrentJetCharge(0), fRefMultWOJet(0), fVZEROMult(0), fMultWOJetVZero(0), fVZero(0), fRefMultFullV0(0), fRefMultV0Corr(0), fFullV0V0Corr(0), fNTXV0MultPt(0), fNTXCBMultPt(0), fMinpTValUE(2.0), fRefMultFullV0UJ(0), fRefMultV0CorrUJ(0), fFullV0V0CorrUJ(0), fMultWOJetVZeroUJ(0), fRefMultWOJetUJ(0), fMaxpTValUE(2.0), fRefAODTrackCount(0), fRefAODTrackCountUJ(0), fTrackCountWOJet(0), fTrackCountWOJetUJ(0), fTrackCountWOJetUJMC(0), fFullV0V0CorrUJMC(0), fMinpTValMC(2.0), fIncExcR(0.0), fForceNotTR(kFALSE), fNotExtDiJEx(kFALSE), fMinTrackPtInNTXRecalc(0), fMaxTrackPtInNTXRecalc(0), fPtDistInJetConeRaw(0), fPtDistInPerpConeRaw(0), fPtInPerpCon(0), fMinTrackPtInNTXR(0), fMaxTrackPtInNTXR(0), fEventCent(0), fJetEtaAll(0), fJetEtaOnlyTPCcut(0), fNChTrRecECorrPPMult(0), fNChTrRecPerpECorrPPMult(0), fForceSkipSJ(kFALSE), fJetPtCentPbPbRaw(0), fJetPtCentPbPbCorr(0), fJetAcceptance(0.5), fIncreasingExcl(kFALSE), fTotTracksCone(0), fTotTracksInCone(0), fTTrackRandomRejection(0), fJTrackRandomRejection(0), fMinPtInGlobMult(0)
147 {
148
149   // Constructor
150
151   for(Int_t a=0; a<16;a++)
152     {
153       fJetFlags[a]=kTRUE;
154       if(a<12)
155         {
156           fNChTr[a]=0;
157           fHistPtParton[a]=0;
158           fSCM[a]=0;
159           if(a<8)
160             {
161               fNChTrRDMult[a]=0;
162               fNAccJetsRDMult[a]=0;
163               fTotalJetCharge[a]=0;
164               fSCMRDMult[a]=0;
165               fNChTrRDMultMC[a]=0;
166               fSCMRDMultMC[a]=0;
167               fNChTrRDMultSE[a]=0;
168               fNAccJetsRDMultSE[a]=0;
169               fTotalJetChargeSE[a]=0;
170               fSCMRDMultSE[a]=0;
171               fNChTrRDMultOJ[a]=0;
172               fSCMRDMultOJ[a]=0;
173               fNChTrRDMultSEOJ[a]=0;
174               fSCMRDMultSEOJ[a]=0;
175               fNChTrRDMultOJMC[a]=0;
176               fSCMRDMultOJMC[a]=0;
177               fNChTrRDMultSEOJMC[a]=0;
178               fSCMRDMultSEOJMC[a]=0;
179               fNChTrRecPerpMultSEOJ[a]=0;               
180               fJEtaMCMultOJ[a]=0;  
181               fJEtaMCMultSEOJ[a]=0;
182               fJEtaRDMultOJ[a]=0;           
183               fJEtaRDMultSEOJ[a]=0;         
184               fJetPtMCMultOJ[a]=0;         
185               fJetPtMCMultSEOJ[a]=0;       
186               fJetPtRDMultOJ[a]=0;        
187               fJetPtRDMultSEOJ[a]=0;       
188               fEntriesQuark[a]=0;          
189               fEntriesGluon[a]=0;          
190             }
191           if(a<6)
192             {
193               fProcessPDG[a]=0;
194               fFragPion[a]=0;
195               fFragKaon[a]=0;
196               fFragProton[a]=0;
197               fFragChargedR4[a]=0;
198               fFragChargedR3[a]=0;
199               fFragChargedR2[a]=0;
200               fHistContainerR4[a]=0;
201               fHistContainerR3[a]=0;
202               fHistContainerR2[a]=0;
203               if(a<3)
204                 {
205                   fJetEtaJetPt[a]=0;
206                   if(a<2)
207                     {
208                       fFragCandidates[a]=0;
209                       fMinTrackPtInNTXh[a]=0;
210                       fMaxTrackPtInNTXh[a]=0; 
211                     }
212                 }
213             }
214         }
215     }
216
217   // Define input and output slots here
218   // Input slot #0 works with a TChain
219   DefineInput(0, TChain::Class());
220   // Output slot #0 id reserved by the base class for AOD
221   // Output slot #1 writes into a TH1 container
222   DefineOutput(1, TList::Class());
223 }
224 //________________________________________________________________________
225 Bool_t AliAnalysisTaskPartonDisc::UserNotify()
226 {
227   //
228   // read the cross sections
229   // and number of trials from pyxsec.root
230   // from AliAnalysisTaskJetSpectrum2
231   //
232
233   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
234   Float_t xsection = 0;
235   Float_t ftrials  = 1;
236   Int_t nevents = 0;
237
238   fAvgTrials = 1;
239   if(tree)
240     {
241       TFile *curfile = tree->GetCurrentFile();
242       if (!curfile) 
243         {
244           Error("Notify","No current file");
245           return kFALSE;
246         }
247       if(!fh1Xsec||!fh1Trials)
248         {
249           Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
250           return kFALSE;
251         }
252       AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
253       fh1Xsec->Fill("<#sigma>",xsection);
254       // construct a poor man average trials (per event!?)
255       Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
256       if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
257       // number of events read out to create the AOD
258       NumberOfReadEventsAOD(curfile->GetName(),nevents);
259       fh1Events->Fill("#sum{nevents}",nevents); //  filled once per file
260     }  
261   return kTRUE;
262
263 }
264 //________________________________________________________________________
265 void AliAnalysisTaskPartonDisc::UserCreateOutputObjects()
266 {
267   // Create histograms
268   // Called once  
269
270   fOutputList = new TList();
271   fOutputList->SetOwner(kTRUE);
272
273   fJetPt = new TH1F("fJetPt", "p_{T} distribution of reco jets", 60, 0., 300.);
274   fJetPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
275   fJetPt->GetYaxis()->SetTitle("dN/dp_{T} (c/GeV)");
276   fJetPt->Sumw2();
277   fOutputList->Add(fJetPt);
278
279   fJetPtSec = new TH1F("fJetPtSec", "p_{T} distribution of reco jets, seconday rec branch", 60, 0., 300.);
280   fJetPtSec->GetXaxis()->SetTitle("p_{T} (GeV/c)");
281   fJetPtSec->GetYaxis()->SetTitle("dN/dp_{T} (c/GeV)");
282   fJetPtSec->Sumw2();
283   fOutputList->Add(fJetPtSec);
284
285   fJetPtMC = new TH1F("fJetPtMC", "p_{T} distribution of mc jets", 60, 0., 300.);
286   fJetPtMC->GetXaxis()->SetTitle("p_{T} (GeV/c)");
287   fJetPtMC->GetYaxis()->SetTitle("dN/dp_{T} (c/GeV)");
288   fJetPtMC->Sumw2();
289   fOutputList->Add(fJetPtMC);
290
291   fJetEta = new TH2F("fJetEta", "Eta distribution of reconstructed jets", 50, -1.5, 1.5, 50, -1.5, 1.5);
292   fJetEta->GetXaxis()->SetTitle("#eta");
293   fJetEta->GetYaxis()->SetTitle("#eta");
294   fJetEta->Sumw2();
295   fOutputList->Add(fJetEta);
296
297   fJetEtaSec = new TH2F("fJetEtaSec", "Eta distribution of reconstructed jets, secondary branch", 50, -1.5, 1.5, 50, -1.5, 1.5);
298   fJetEtaSec->GetXaxis()->SetTitle("#eta");
299   fJetEtaSec->GetYaxis()->SetTitle("#eta");
300   fJetEtaSec->Sumw2();
301   fOutputList->Add(fJetEtaSec);
302
303   fJetPhi = new TH2F("fJetPhi", "Phi distribution of reconstructed jets", 50, 0., TMath::TwoPi(), 50, 0., TMath::TwoPi());
304   fJetPhi->GetXaxis()->SetTitle("#phi");
305   fJetPhi->GetYaxis()->SetTitle("#phi");
306   fJetPhi->Sumw2();
307   fOutputList->Add(fJetPhi);
308
309   fJetPhiSec = new TH2F("fJetPhiSec", "Phi distribution of reconstructed jets, secondary branch", 50, 0., TMath::TwoPi(), 50, 0., TMath::TwoPi());
310   fJetPhiSec->GetXaxis()->SetTitle("#phi");
311   fJetPhiSec->GetYaxis()->SetTitle("#phi");
312   fJetPhiSec->Sumw2();
313   fOutputList->Add(fJetPhiSec);
314
315   fJetEtaMC = new TH2F("fJetEtaMC", "Eta distribution of MC jets", 50, -1.5, 1.5, 50, -1.5, 1.5);
316   fJetEtaMC->GetXaxis()->SetTitle("#eta");
317   fJetEtaMC->GetYaxis()->SetTitle("#eta");
318   fJetEtaMC->Sumw2();
319   fOutputList->Add(fJetEtaMC);
320
321   fJetPhiMC = new TH2F("fJetPhiMC", "Phi distribution of MC jets", 50, 0., TMath::TwoPi(), 50, 0., TMath::TwoPi());
322   fJetPhiMC->GetXaxis()->SetTitle("#phi");
323   fJetPhiMC->GetYaxis()->SetTitle("#phi");
324   fJetPhiMC->Sumw2();
325   fOutputList->Add(fJetPhiMC);
326
327   fPtAODMC = new TH2F("fPtAODMC", "p_{T} distribution of mc tracks", 200, 0., 100., 200, 0., 100.);
328   fPtAODMC->GetXaxis()->SetTitle("p_{T} (GeV/c)");
329   fPtAODMC->GetYaxis()->SetTitle("p_{T} (GeV/c)");
330   fPtAODMC->Sumw2();
331   fOutputList->Add(fPtAODMC);
332
333   fPtAOD = new TH2F("fPtAOD", "p_{T} distribution of aod tracks", 200, 0., 100., 200, 0., 100.);
334   fPtAOD->GetXaxis()->SetTitle("p_{T} (GeV/c)");
335   fPtAOD->GetYaxis()->SetTitle("p_{T} (GeV/c)");
336   fPtAOD->Sumw2();
337   fOutputList->Add(fPtAOD);
338
339   fEtaAODMC = new TH2F("fEtaAODMC", "Eta distribution of MC aod tracks", 50, -1.5, 1.5, 50, -1.5, 1.5);
340   fEtaAODMC->GetXaxis()->SetTitle("#eta");
341   fEtaAODMC->GetYaxis()->SetTitle("#eta");
342   fEtaAODMC->Sumw2();
343   fOutputList->Add(fEtaAODMC);
344
345   fPhiAODMC = new TH2F("fPhiAODMC", "Phi distribution of MC aod tracks", 50, 0., TMath::TwoPi(), 50, 0., TMath::TwoPi());
346   fPhiAODMC->GetXaxis()->SetTitle("#phi");
347   fPhiAODMC->GetYaxis()->SetTitle("#phi");
348   fPhiAODMC->Sumw2();
349   fOutputList->Add(fPhiAODMC);
350
351   fEtaAOD = new TH2F("fEtaAOD", "Eta distribution of aod tracks", 50, -1.5, 1.5, 50, -1.5, 1.5);
352   fEtaAOD->GetXaxis()->SetTitle("#eta");
353   fEtaAOD->GetYaxis()->SetTitle("#eta");
354   fEtaAOD->Sumw2();
355   fOutputList->Add(fEtaAOD);
356
357   fPhiAOD = new TH2F("fPhiAOD", "Phi distribution of aod tracks", 50, 0.,TMath::TwoPi(), 50, 0.,TMath::TwoPi());
358   fPhiAOD->GetXaxis()->SetTitle("#phi");
359   fPhiAOD->GetYaxis()->SetTitle("#phi");
360   fPhiAOD->Sumw2();
361   fOutputList->Add(fPhiAOD);
362
363   fFlavor = new TH2F("fFlavor", "Flavor distribution of jets", 27, -5.5, 21.5, 60, 0., 300.);
364   fFlavor->GetXaxis()->SetTitle("PDG code");
365   fFlavor->GetYaxis()->SetTitle("p_{T}^{JET}");
366   fFlavor->Sumw2();
367   fOutputList->Add(fFlavor);
368
369   fNJetsMC = new TH2F("fNJetsMC", "Number of simulated jets per event, as a function of p_T", 101, -0.5, 100.5, 101, -0.5, 100.5);
370   fNJetsMC->GetXaxis()->SetTitle("Number of jets");
371   fNJetsMC->GetYaxis()->SetTitle("Number of jets");
372   fNJetsMC->Sumw2();
373   fOutputList->Add(fNJetsMC);
374
375   fNJetsRD = new TH2F("fNJetsRD", "Number of jets per event in real data", 101, -0.5, 100.5, 101, -0.5, 100.5);
376   fNJetsRD->GetXaxis()->SetTitle("Number of jets");
377   fNJetsRD->GetYaxis()->SetTitle("Number of jets");
378   fNJetsRD->Sumw2();
379   fOutputList->Add(fNJetsRD);
380
381   fNJetsRDSeco = new TH2F("fNJetsRDSeco", "Number of jets per event in real data, secondary branch", 101, -0.5, 100.5, 101, -0.5, 100.5);
382   fNJetsRDSeco->GetXaxis()->SetTitle("Number of jets");
383   fNJetsRDSeco->GetYaxis()->SetTitle("Number of jets");
384   fNJetsRDSeco->Sumw2();
385   fOutputList->Add(fNJetsRDSeco);
386
387   fJetsMultPtMC = new TH2F("fJetsMultPtMC", "Jet multiplicity associated to jet pT, as a function of p_T MC", 60, 0., 300., 7, -0.5, 6.5);
388   fJetsMultPtMC->GetXaxis()->SetTitle("p_{T}^{JET}");
389   fJetsMultPtMC->GetYaxis()->SetTitle("Number of jets");
390   fJetsMultPtMC->Sumw2();
391   fOutputList->Add(fJetsMultPtMC);
392
393   fJetsMultPtRD = new TH2F("fJetsMultPtRD", "Jet multiplicity associated to jet pT, as a function of p_T Reco Data",60, 0., 300., 7, -0.5, 6.5);
394   fJetsMultPtRD->GetXaxis()->SetTitle("p_{T}^{JET}");
395   fJetsMultPtRD->GetYaxis()->SetTitle("Number of jets");
396   fJetsMultPtRD->Sumw2();
397   fOutputList->Add(fJetsMultPtRD);
398
399   fNChTrRD = new TH2F("fNChTrRD","Number of tracks to recover transverse energy, jet_{p_{T}}",101,-0.5,100.5, 60, 0., 300.);
400   fNChTrRD->GetXaxis()->SetTitle("NTracks");
401   fNChTrRD->GetYaxis()->SetTitle("p_{T}^{JET}");
402   fNChTrRD->Sumw2();
403   fOutputList->Add(fNChTrRD);
404
405   fProfNChTrRD = new TProfile("fProfNChTrRD","Number of tracks to recover transverse energy, jet_{p_{T}}", 50, 0, 250);
406   fProfNChTrRD->GetXaxis()->SetTitle("p_{T}^{JET}");
407   fProfNChTrRD->GetYaxis()->SetTitle("NTracks");
408   fProfNChTrRD->Sumw2();
409   fOutputList->Add(fProfNChTrRD);
410
411   fFracQQ = new TH1F("fFracQQ","QQ",1000,0.,0.5);
412   fFracQQ->GetXaxis()->SetTitle("x_{T}");
413   fFracQQ->GetYaxis()->SetTitle("Entries");
414   fFracQQ->Sumw2();
415   fOutputList->Add(fFracQQ);
416
417   fFracGQ = new TH1F("fFracGQ","GQ",1000,0.,0.5);
418   fFracGQ->GetXaxis()->SetTitle("x_{T}");
419   fFracGQ->GetYaxis()->SetTitle("Entries");
420   fFracGQ->Sumw2();
421   fOutputList->Add(fFracGQ);
422
423   fFracGG = new TH1F("fFracGG","GG",1000,0.,0.5);
424   fFracGG->GetXaxis()->SetTitle("x_{T}");
425   fFracGG->GetYaxis()->SetTitle("Entries");
426   fFracGG->Sumw2();
427   fOutputList->Add(fFracGG);
428
429   fFracOutGoingQQ = new TH1F("fFracOutGoingQQ","QQ",1000,0.,0.5);
430   fFracOutGoingQQ->GetXaxis()->SetTitle("x_{T}");
431   fFracOutGoingQQ->GetYaxis()->SetTitle("Entries");
432   fFracOutGoingQQ->Sumw2();
433   fOutputList->Add(fFracOutGoingQQ);
434
435   fFracOutGoingGQ = new TH1F("fFracOutGoingGQ","GQ",1000,0.,0.5);
436   fFracOutGoingGQ->GetXaxis()->SetTitle("x_{T}");
437   fFracOutGoingGQ->GetYaxis()->SetTitle("Entries");
438   fFracOutGoingGQ->Sumw2();
439   fOutputList->Add(fFracOutGoingGQ);
440
441   fFracOutGoingGG = new TH1F("fFracOutGoingGG","GG",1000,0.,0.5);
442   fFracOutGoingGG->GetXaxis()->SetTitle("x_{T}");
443   fFracOutGoingGG->GetYaxis()->SetTitle("Entries");
444   fFracOutGoingGG->Sumw2();
445   fOutputList->Add(fFracOutGoingGG);
446
447   fh1Xsec =  new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1);
448   fh1Xsec->SetXTitle("<#sigma>");
449   fh1Xsec->Sumw2();
450   fOutputList->Add(fh1Xsec);  
451
452   fh1Trials = new TH1F("h1Trials","trials from pyxsec.root",1,0,1);
453   fh1Trials->SetXTitle("#sum{ntrials}");
454   fh1Trials->Sumw2();
455   fOutputList->Add(fh1Trials);
456            
457   fProcessJetPt = new TH2F("fProcessJetPt","Pythia process number, jet_{p_{T}}",70, 0.5, 70.5, 60, 0., 300.);
458   fProcessJetPt->GetXaxis()->SetTitle("Pythia process");
459   fProcessJetPt->GetYaxis()->SetTitle("p_{T}^{JET}");
460   fProcessJetPt->Sumw2();
461   fOutputList->Add(fProcessJetPt);
462
463   fFlavorLead = new TH2F("fFlavorLead", "Flavor distribution of leading jets", 27, -5.5, 21.5, 60, 0., 300.);
464   fFlavorLead->GetXaxis()->SetTitle("PDG code");
465   fFlavorLead->GetYaxis()->SetTitle("p_{T}^{JET}");
466   fFlavorLead->Sumw2();
467   fOutputList->Add(fFlavorLead);
468
469   fProcessLeadJetPt = new TH2F("fProcessLeadJetPt","Pythia process number for leading, jet_{p_{T}}",70, 0.5, 70.5, 60, 0., 300.);
470   fProcessLeadJetPt->GetXaxis()->SetTitle("Pythia process");
471   fProcessLeadJetPt->GetYaxis()->SetTitle("p_{T}^{JET}");
472   fProcessLeadJetPt->Sumw2();
473   fOutputList->Add(fProcessLeadJetPt);
474  
475   fPDGMothLPart = new TH3F("fPDGMothLPart","Mother of leading parton, leading parton, jet p_{T}", 27, -5.5, 21.5, 27, -5.5, 21.5, 60, 0., 300.);
476   fPDGMothLPart->GetXaxis()->SetTitle("Mother of leading parton");
477   fPDGMothLPart->GetYaxis()->SetTitle("Leading parton");
478   fPDGMothLPart->GetZaxis()->SetTitle("p_{T}^{JET}");
479   fPDGMothLPart->Sumw2();
480   fOutputList->Add(fPDGMothLPart);
481   
482   fFlavProc = new TH2F("fFlavProc","Flavor, Flavor status code", 27, -5.5, 21.5, 101, -0.5, 100.5);
483   fFlavProc->GetXaxis()->SetTitle("Jet flavor");
484   fFlavProc->GetYaxis()->SetTitle("Parton status code");
485   fFlavProc->Sumw2();
486   fOutputList->Add(fFlavProc);
487
488   fSCMRD = new TH2F("fSCMRD","Second Central Moment, jet_{p_{T}}",200,0.,0.2, 60, 0., 300.);
489   fSCMRD->GetXaxis()->SetTitle("<#delta R_{c}^{2}>");
490   fSCMRD->GetYaxis()->SetTitle("p_{T}^{JET}");
491   fSCMRD->Sumw2();
492   fOutputList->Add(fSCMRD);
493  
494   fZVertex = new TH2F("fZVertex","Z vertex position, Z vertex position}",40,-20.,20., 40, -20., 20.);
495   fZVertex->GetXaxis()->SetTitle("Vertex z");
496   fZVertex->GetYaxis()->SetTitle("Vertex z");
497   fZVertex->Sumw2();
498   fOutputList->Add(fZVertex);
499
500   fh1Events = new TH1F("fh1Events","nevents from PWG4_JetTasksOutput.root",1,0,1);
501   fh1Events->SetXTitle("#sum{nevents}");
502   fh1Events->Sumw2();
503   fOutputList->Add(fh1Events);
504
505   fNAccJetsMC = new TH2F("fNAccJetsMC", "Number accepted simulated jets per event", 101, -0.5, 100.5, 101, -0.5, 100.5);
506   fNAccJetsMC->GetXaxis()->SetTitle("Number of jets");
507   fNAccJetsMC->GetYaxis()->SetTitle("Number of jets");
508   fNAccJetsMC->Sumw2();
509   fOutputList->Add(fNAccJetsMC);
510
511   fNAccJetsRD = new TH2F("fNAccJetsRD", "Number of accepted jets per event in real data", 101, -0.5, 100.5, 101, -0.5, 100.5);
512   fNAccJetsRD->GetXaxis()->SetTitle("Number of jets");
513   fNAccJetsRD->GetYaxis()->SetTitle("Number of jets");
514   fNAccJetsRD->Sumw2();
515   fOutputList->Add(fNAccJetsRD);
516
517   fNAccJetsRDSeco = new TH2F("fNAccJetsRDSeco", "Number of accepted jets per event in real data, secondary branch", 101, -0.5, 100.5, 101, -0.5, 100.5);
518   fNAccJetsRDSeco->GetXaxis()->SetTitle("Number of jets");
519   fNAccJetsRDSeco->GetYaxis()->SetTitle("Number of jets");
520   fNAccJetsRDSeco->Sumw2();
521   fOutputList->Add(fNAccJetsRDSeco);
522
523   fRecJetPtInclusive = new TH1F("fRecJetPtInclusive", "p_{T} distribution of inclusive reco jets", 60, 0., 300.);
524   fRecJetPtInclusive->GetXaxis()->SetTitle("p_{T} (GeV/c)");
525   fRecJetPtInclusive->GetYaxis()->SetTitle("d#sigma (mb)");
526   fRecJetPtInclusive->Sumw2();
527   fOutputList->Add(fRecJetPtInclusive);
528
529   fMCJetPtInclusive = new TH1F("fMCJetPtInclusive", "p_{T} distribution of inclusive MC jets", 60, 0., 300.);
530   fMCJetPtInclusive->GetXaxis()->SetTitle("p_{T} (GeV/c)");
531   fMCJetPtInclusive->GetYaxis()->SetTitle("d#sigma (mb)");
532   fMCJetPtInclusive->Sumw2();
533   fOutputList->Add(fMCJetPtInclusive);
534
535   fRecJetPtLeading = new TH1F("fRecJetPtLeading", "p_{T} distribution of leading reco jets", 60, 0., 300.);
536   fRecJetPtLeading->GetXaxis()->SetTitle("p_{T} (GeV/c)");
537   fRecJetPtLeading->GetYaxis()->SetTitle("d#sigma (mb)");
538   fRecJetPtLeading->Sumw2();
539   fOutputList->Add(fRecJetPtLeading);
540
541   fMCJetPtLeading = new TH1F("fMCJetPtLeading", "p_{T} distribution of leading MC jets", 60, 0., 300.);
542   fMCJetPtLeading->GetXaxis()->SetTitle("p_{T} (GeV/c)");
543   fMCJetPtLeading->GetYaxis()->SetTitle("d#sigma (mb)");
544   fMCJetPtLeading->Sumw2();
545   fOutputList->Add(fMCJetPtLeading);
546
547   fSecRecJetPtInclusive = new TH1F("fSecRecJetPtInclusive", "p_{T} distribution of inclusive reco jets (2nd branch)", 60, 0., 300.);
548   fSecRecJetPtInclusive->GetXaxis()->SetTitle("p_{T} (GeV/c)");
549   fSecRecJetPtInclusive->GetYaxis()->SetTitle("d#sigma (mb)");
550   fSecRecJetPtInclusive->Sumw2();
551   fOutputList->Add(fSecRecJetPtInclusive);
552
553   fSecRecJetPtLeading = new TH1F("fSecRecJetPtLeading", "p_{T} distribution of leading reco jets (2nd branch)", 60, 0., 300.);
554   fSecRecJetPtLeading->GetXaxis()->SetTitle("p_{T} (GeV/c)");
555   fSecRecJetPtLeading->GetYaxis()->SetTitle("d#sigma (mb)");
556   fSecRecJetPtLeading->Sumw2();
557   fOutputList->Add(fSecRecJetPtLeading);
558
559   fNChTrMCPerp = new TH2F("fNChTrMCPerp","Number of tracks to recover transverse energy of perp. cone, jet_{p_{T} MC}",101,-0.5,100.5, 60, 0., 300.);
560   fNChTrMCPerp->GetXaxis()->SetTitle("NTracks Perp");
561   fNChTrMCPerp->GetYaxis()->SetTitle("p_{T}^{MC JET}");
562   fNChTrMCPerp->Sumw2();
563   fOutputList->Add(fNChTrMCPerp);
564
565   fNChTrRecPerp = new TH2F("fNChTrRecPerp","Number of tracks to recover transverse energy of perp. cone, jet_{p_{T} Rec}",101,-0.5,100.5, 60, 0., 300.);
566   fNChTrRecPerp->GetXaxis()->SetTitle("NTracks Perp");
567   fNChTrRecPerp->GetYaxis()->SetTitle("p_{T}^{RECO JET}");
568   fNChTrRecPerp->Sumw2();
569   fOutputList->Add(fNChTrRecPerp);
570
571   fSCMMCPerp = new TH2F("fSCMMCPerp","Second Central Moment of perp. cone, jet_{p_{T} MC}",200,0.,0.2, 60, 0., 300.);
572   fSCMMCPerp->GetXaxis()->SetTitle("<#delta R_{c}^{2}> Perp");
573   fSCMMCPerp->GetYaxis()->SetTitle("p_{T}^{JET}");
574   fSCMMCPerp->Sumw2();
575   fOutputList->Add(fSCMMCPerp);
576
577   fSCMRecPerp = new TH2F("fSCMRecPerp","Second Central Moment of perp. cone, jet_{p_{T} Reco}",200,0.,0.2, 60, 0., 300.);
578   fSCMRecPerp->GetXaxis()->SetTitle("<#delta R_{c}^{2}> Perp");
579   fSCMRecPerp->GetYaxis()->SetTitle("p_{T}^{JET}");
580   fSCMRecPerp->Sumw2();
581   fOutputList->Add(fSCMRecPerp);
582
583   fNChTrCorrMCQuark = new TH2F("fNChTrCorrMCQuark","Number of tracks to recover corrected transverse energy, MC quarks",101,-0.5,100.5, 60, 0., 300.);
584   fNChTrCorrMCQuark->GetXaxis()->SetTitle("NTracks");
585   fNChTrCorrMCQuark->GetYaxis()->SetTitle("p_{T}^{MC Corr. JET}");
586   fNChTrCorrMCQuark->Sumw2();
587   fOutputList->Add(fNChTrCorrMCQuark);
588
589   fNChTrCorrMCGluon = new TH2F("fNChTrCorrMCGluon","Number of tracks to recover corrected transverse energy, MC gluons",101,-0.5,100.5, 60, 0., 300.);
590   fNChTrCorrMCGluon->GetXaxis()->SetTitle("NTracks");
591   fNChTrCorrMCGluon->GetYaxis()->SetTitle("p_{T}^{MC Corr. JET}");
592   fNChTrCorrMCGluon->Sumw2();
593   fOutputList->Add(fNChTrCorrMCGluon);
594
595   fNChTrCorrMCPerp = new TH2F("fNChTrCorrMCPerp","Number of tracks to recover perp. cone. after corrected jet pT",101,-0.5,100.5, 60, 0., 300.);
596   fNChTrCorrMCPerp->GetXaxis()->SetTitle("NTracks");
597   fNChTrCorrMCPerp->GetYaxis()->SetTitle("p_{T}^{MC Corr. JET}");
598   fNChTrCorrMCPerp->Sumw2();
599   fOutputList->Add(fNChTrCorrMCPerp);
600
601   // 9 selection bins: (nuevo)
602   // 1st. Proton collisions                  fill 1  Bin1 [0.5,1.5)
603   // 2nd. PbPb collisions, Bin  0-10         fill 2  Bin2 [1.5,2.5)
604   // 3rd. PbPb collisions, Bin 10-20         fill 3  Bin3 [2.5,3.5)
605   // 4rd. PbPb collisions, Bin 20-30         fill 4  Bin4 [3.5,4.5)
606   // 5th. PbPb collisions, Bin 30-40         fill 5  Bin5 [4.5,5.5)
607   // 6th. PbPb collisions, Bin 40-50         fill 6  Bin6 [5.5,6.5)
608   // 7th. PbPb collisions, Bin 50-60         fill 7  Bin7 [6.5,7.5)
609   // 8th. PbPb collisions, Bin 60-70         fill 8  Bin8 [7.5,8.5)
610   // 9th. PbPb collisions, Bin 70-80         fill 9  Bin9 [8.5,9.5)
611   // 10th. PbPb collisions, Bin 80-100.1    fill 10  Bin10 [9.5,10.5)
612
613   fNChTrRecECorr = new TH3F("fNChTrRecECorr","NTX in ener. corr. jet , corr. jet pT, centrality",101,-0.5,100.5, 60, 0., 300.,10,0.5,10.5);
614   fNChTrRecECorr->GetXaxis()->SetTitle("NTracks");
615   fNChTrRecECorr->GetYaxis()->SetTitle("p_{T}^{JET}");
616   fNChTrRecECorr->GetZaxis()->SetTitle("Selection Bin");
617   fNChTrRecECorr->Sumw2();
618   fOutputList->Add(fNChTrRecECorr);
619
620   fNChTrRecPerpECorr = new TH3F("fNChTrRecPerpECorr","Tracks above min in perp.cone , corr. jet pT, centrality",101,-0.5,100.5, 60, 0., 300.,10,0.5,10.5);
621   fNChTrRecPerpECorr->GetXaxis()->SetTitle("NTracks");
622   fNChTrRecPerpECorr->GetYaxis()->SetTitle("p_{T}^{JET}");
623   fNChTrRecPerpECorr->GetZaxis()->SetTitle("Selection Bin");
624   fNChTrRecPerpECorr->Sumw2();
625   fOutputList->Add(fNChTrRecPerpECorr);
626   
627   fRefMult = new TH1F("fRefMult", "Reference multiplicity in the AOD", 301, -0.5, 300.5);
628   fRefMult->GetXaxis()->SetTitle("Reference multiplicity");
629   fRefMult->Sumw2();
630   fOutputList->Add(fRefMult);
631
632   fRefMultWOJet = new TH2F("fRefMultWOJet", "Reference multiplicity in the AOD, multiplicity without jets", 301, -0.5, 300.5, 301, -0.5, 300.5);
633   fRefMultWOJet->GetXaxis()->SetTitle("Reference multiplicity");
634   fRefMultWOJet->GetYaxis()->SetTitle("Multiplicity without jets");
635   fRefMultWOJet->Sumw2();
636   fOutputList->Add(fRefMultWOJet);
637
638   fVZEROMult = new TH2F("fVZEROMult", "Multiplicity V0A and V0C", 501, -0.5, 500.5, 501, -0.5, 500.5);
639   fVZEROMult->GetXaxis()->SetTitle("Multiplicity V0A");
640   fVZEROMult->GetYaxis()->SetTitle("Multiplicity V0C");
641   fVZEROMult->Sumw2();
642   fOutputList->Add(fVZEROMult);
643
644   fMultWOJetVZero = new TH2F("fMultWOJetVZero", "Multiplicity without jets and VZERO mult.",301, -0.5, 300.5, 1001, -0.5, 1000.5);
645   fMultWOJetVZero->GetXaxis()->SetTitle("Multiplicity without jets TPC");
646   fMultWOJetVZero->GetYaxis()->SetTitle("Multiplicity full V0");
647   fMultWOJetVZero->Sumw2();
648   fOutputList->Add(fMultWOJetVZero);
649
650   fRefMultFullV0 = new TH2F("fRefMultFullV0", "Reference multiplicity in the AOD, multiplicity from full V0",301, -0.5, 300.5, 1001, -0.5, 1000.5);
651   fRefMultFullV0->GetXaxis()->SetTitle("Reference multiplicity in AOD");
652   fRefMultFullV0->GetYaxis()->SetTitle("Multiplicity full V0");
653   fRefMultFullV0->Sumw2();
654   fOutputList->Add(fRefMultFullV0);
655
656   fRefMultV0Corr = new TH2F("fRefMultV0Corr", "Reference multiplicity in the AOD, multiplicity from corrected V0",301, -0.5, 300.5, 1001, -0.5, 1000.5);
657   fRefMultV0Corr->GetXaxis()->SetTitle("Reference multiplicity in AOD");
658   fRefMultV0Corr->GetYaxis()->SetTitle("Multiplicity V0 no jets");
659   fRefMultV0Corr->Sumw2();
660   fOutputList->Add(fRefMultV0Corr);
661
662   fFullV0V0Corr = new TH2F("fFullV0V0Corr", "Multiplicity from full V0, multiplicity from corrected V0",1001, -0.5, 1000.5, 1001, -0.5, 1000.5);
663   fFullV0V0Corr->GetXaxis()->SetTitle("Multiplicity from full V0");
664   fFullV0V0Corr->GetYaxis()->SetTitle("Multiplicity V0 no jets");
665   fFullV0V0Corr->Sumw2();
666   fOutputList->Add(fFullV0V0Corr);
667
668   fNTXV0MultPt = new TH3F("fNTXV0MultPt", "NTX, Multiplicity from corrected V0, jet pT",101,-0.5,100.5, 1001, -0.5, 1000.5, 60, 0., 300.);
669   fNTXV0MultPt->GetXaxis()->SetTitle("NTracks");
670   fNTXV0MultPt->GetYaxis()->SetTitle("Multiplicity V0 no jets");
671   fNTXV0MultPt->GetZaxis()->SetTitle("p_{T}^{JET}");
672   fNTXV0MultPt->Sumw2();
673   fOutputList->Add(fNTXV0MultPt);
674
675   fNTXCBMultPt = new TH3F("fNTXCBMultPt", "NTX, Multiplicity from corrected Central Barrel, jet pT",101,-0.5,100.5, 301, -0.5, 300.5, 60, 0., 300.);
676   fNTXCBMultPt->GetXaxis()->SetTitle("NTracks");
677   fNTXCBMultPt->GetYaxis()->SetTitle("Multiplicity corrected Central Barrel");
678   fNTXCBMultPt->GetZaxis()->SetTitle("p_{T}^{JET}");
679   fNTXCBMultPt->Sumw2();
680   fOutputList->Add(fNTXCBMultPt);
681
682   fRefMultFullV0UJ = new TH2F("fRefMultFullV0UJ", "Reference multiplicity in the AOD, multiplicity from full V0, 1 jet event",301, -0.5, 300.5, 1001, -0.5, 1000.5);
683   fRefMultFullV0UJ->GetXaxis()->SetTitle("Reference multiplicity in AOD");
684   fRefMultFullV0UJ->GetYaxis()->SetTitle("Multiplicity full V0");
685   fRefMultFullV0UJ->Sumw2();
686   fOutputList->Add(fRefMultFullV0UJ);
687
688   fRefMultV0CorrUJ = new TH2F("fRefMultV0CorrUJ", "Reference multiplicity in the AOD, multiplicity from corrected V0, 1 jet event",301, -0.5, 300.5, 1001, -0.5, 1000.5);
689   fRefMultV0CorrUJ->GetXaxis()->SetTitle("Reference multiplicity in AOD");
690   fRefMultV0CorrUJ->GetYaxis()->SetTitle("Multiplicity V0 no jets");
691   fRefMultV0CorrUJ->Sumw2();
692   fOutputList->Add(fRefMultV0CorrUJ);
693
694   fFullV0V0CorrUJ = new TH2F("fFullV0V0CorrUJ", "Multiplicity from full V0, multiplicity from corrected V0, 1 jet event",1001, -0.5, 1000.5, 1001, -0.5, 1000.5);
695   fFullV0V0CorrUJ->GetXaxis()->SetTitle("Multiplicity from full V0");
696   fFullV0V0CorrUJ->GetYaxis()->SetTitle("Multiplicity V0 no jets");
697   fFullV0V0CorrUJ->Sumw2();
698   fOutputList->Add(fFullV0V0CorrUJ);
699
700   fMultWOJetVZeroUJ = new TH2F("fMultWOJetVZeroUJ", "Multiplicity without jets and VZERO mult., 1 jet event",301, -0.5, 300.5, 1001, -0.5, 1000.5);
701   fMultWOJetVZeroUJ->GetXaxis()->SetTitle("Multiplicity without jets TPC");
702   fMultWOJetVZeroUJ->GetYaxis()->SetTitle("Multiplicity full V0");
703   fMultWOJetVZeroUJ->Sumw2();
704   fOutputList->Add(fMultWOJetVZeroUJ);
705
706   fRefMultWOJetUJ = new TH2F("fRefMultWOJetUJ", "Reference multiplicity in the AOD, multiplicity without jets, 1 jet event", 301, -0.5, 300.5, 301, -0.5, 300.5);
707   fRefMultWOJetUJ->GetXaxis()->SetTitle("Reference multiplicity");
708   fRefMultWOJetUJ->GetYaxis()->SetTitle("Multiplicity without jets");
709   fRefMultWOJetUJ->Sumw2();
710   fOutputList->Add(fRefMultWOJetUJ);
711
712   fRefAODTrackCount = new TH2F("fRefAODTrackCount", "Reference multiplicity in the AOD, my own referece mult.", 301, -0.5, 300.5, 301, -0.5, 300.5);
713   fRefAODTrackCount->GetXaxis()->SetTitle("AOD Reference multiplicity");
714   fRefAODTrackCount->GetYaxis()->SetTitle("My Reference multiplicity");
715   fRefAODTrackCount->Sumw2();
716   fOutputList->Add(fRefAODTrackCount);
717
718   fRefAODTrackCountUJ = new TH2F("fRefAODTrackCountUJ", "Reference multiplicity in the AOD, my own referece mult., 1 jet event", 301, -0.5, 300.5, 301, -0.5, 300.5);
719   fRefAODTrackCountUJ->GetXaxis()->SetTitle("AOD Reference multiplicity");
720   fRefAODTrackCountUJ->GetYaxis()->SetTitle("My Reference multiplicity");
721   fRefAODTrackCountUJ->Sumw2();
722   fOutputList->Add(fRefAODTrackCountUJ);
723
724   fTrackCountWOJet = new TH2F("fTrackCountWOJet", "My own total referece mult., soft mult", 151, -0.5, 150.5, 151, -0.5, 150.5);
725   fTrackCountWOJet->GetXaxis()->SetTitle("Total TPC multiplicity");
726   fTrackCountWOJet->GetYaxis()->SetTitle("Soft TPC multiplicity");
727   fTrackCountWOJet->Sumw2();
728   fOutputList->Add(fTrackCountWOJet);
729
730   fTrackCountWOJetUJ = new TH2F("fTrackCountWOJetUJ", "My own total referece mult., soft mult, 1 jet", 151, -0.5, 150.5, 151, -0.5, 150.5);
731   fTrackCountWOJetUJ->GetXaxis()->SetTitle("Total TPC multiplicity");
732   fTrackCountWOJetUJ->GetYaxis()->SetTitle("Soft TPC multiplicity");
733   fTrackCountWOJetUJ->Sumw2();
734   fOutputList->Add(fTrackCountWOJetUJ);
735
736   fTrackCountWOJetUJMC = new TH2F("fTrackCountWOJetUJMC", "My own total referece mult., soft mult, 1 jet, MC!", 151, -0.5, 150.5, 151, -0.5, 150.5);
737   fTrackCountWOJetUJMC->GetXaxis()->SetTitle("Total TPC (eta) multiplicity");
738   fTrackCountWOJetUJMC->GetYaxis()->SetTitle("Soft TPC (eta) multiplicity");
739   fTrackCountWOJetUJMC->Sumw2();
740   fOutputList->Add(fTrackCountWOJetUJMC);
741
742   fFullV0V0CorrUJMC = new TH2F("fFullV0V0CorrUJMC", "Multiplicity from full V0, multiplicity from corrected V0, 1 jet event, MC!",1001, -0.5, 1000.5, 1001, -0.5, 1000.5);
743   fFullV0V0CorrUJMC->GetXaxis()->SetTitle("Multiplicity from full V0 (acceptance)");
744   fFullV0V0CorrUJMC->GetYaxis()->SetTitle("Multiplicity V0(acceptance) no jets");
745   fFullV0V0CorrUJMC->Sumw2();
746   fOutputList->Add(fFullV0V0CorrUJMC);
747
748   fMinTrackPtInNTXRecalc = new TH3F("fMinTrackPtInNTXRecalc", "Minimum track pT for the jets after pT correction, raw jet pT", 200, 0., 100., 60, 0., 300.,10,0.5,10.5);
749   fMinTrackPtInNTXRecalc->GetXaxis()->SetTitle("p_{T}^{TRACK} (GeV/c)");
750   fMinTrackPtInNTXRecalc->GetYaxis()->SetTitle("p_{T}^{JET} (GeV/c)");
751   fMinTrackPtInNTXRecalc->GetZaxis()->SetTitle("Selection Bin");  // 9 selections bins as fNChTrRecECorr
752   fMinTrackPtInNTXRecalc->Sumw2();
753   fOutputList->Add(fMinTrackPtInNTXRecalc);
754
755   fMaxTrackPtInNTXRecalc = new TH2F("fMaxTrackPtInNTXRecalc", "Maximum track pT for the jets after pT correction, raw jet pT", 200, 0., 100., 60, 0., 300.);
756   fMaxTrackPtInNTXRecalc->GetXaxis()->SetTitle("p_{T}^{TRACK} (GeV/c)");
757   fMaxTrackPtInNTXRecalc->GetYaxis()->SetTitle("p_{T}^{JET} (GeV/c)");
758   fMaxTrackPtInNTXRecalc->Sumw2();
759   fOutputList->Add(fMaxTrackPtInNTXRecalc);
760
761   fPtDistInJetConeRaw = new TH3F("fPtDistInJetConeRaw","pT of tracks in cone, raw jet pT bin, centrality", 200, 0., 100., 8, 0.5, 8.5, 10, 0.5, 10.5);
762   fPtDistInJetConeRaw->GetXaxis()->SetTitle("p_{T}^{TRACK} (GeV/c)");
763   fPtDistInJetConeRaw->GetYaxis()->SetTitle("p_{T}^{JET} Bin");
764   fPtDistInJetConeRaw->GetZaxis()->SetTitle("Centrality Bin");
765   fPtDistInJetConeRaw->Sumw2();
766   fOutputList->Add(fPtDistInJetConeRaw);
767
768   fPtDistInPerpConeRaw = new TH3F("fPtDistInPerpConeRaw","pT of tracks in cone, raw jet pT bin, centrality", 200, 0., 100., 8, 0.5, 8.5, 10, 0.5, 10.5);
769   fPtDistInPerpConeRaw->GetXaxis()->SetTitle("p_{T}^{TRACK} (GeV/c)");
770   fPtDistInPerpConeRaw->GetYaxis()->SetTitle("p_{T}^{JET} Bin");
771   fPtDistInPerpConeRaw->GetZaxis()->SetTitle("Centrality Bin");
772   fPtDistInPerpConeRaw->Sumw2();
773   fOutputList->Add(fPtDistInPerpConeRaw);
774
775   fPtInPerpCon = new TH3F("fPtInPerpCon","Summed pT of perpendicular cone, raw jet pT bin, centrality", 200, 0., 100., 8, 0.5, 8.5, 10, 0.5, 10.5);
776   fPtInPerpCon->GetXaxis()->SetTitle("p_{T}^{PERP.CONE} (GeV/c)");
777   fPtInPerpCon->GetYaxis()->SetTitle("p_{T}^{JET} Bin");
778   fPtInPerpCon->GetZaxis()->SetTitle("Centrality Bin");
779   fPtInPerpCon->Sumw2();
780   fOutputList->Add(fPtInPerpCon);
781
782   fJetEtaAll = new TH1F("fJetEtaAll", "Eta distribution of reconstructed jets, no cuts", 50, -1.5, 1.5);
783   fJetEtaAll->GetXaxis()->SetTitle("#eta");
784   fJetEtaAll->GetYaxis()->SetTitle("entries");
785   fJetEtaAll->Sumw2();
786   fOutputList->Add(fJetEtaAll);
787
788   fJetEtaOnlyTPCcut = new TH1F("fJetEtaOnlyTPCcut", "Eta distribution of reconstructed jets, only tpc acceptance cut", 50, -1.5, 1.5);
789   fJetEtaOnlyTPCcut->GetXaxis()->SetTitle("#eta");
790   fJetEtaOnlyTPCcut->GetYaxis()->SetTitle("entries");
791   fJetEtaOnlyTPCcut->Sumw2();
792   fOutputList->Add(fJetEtaOnlyTPCcut);
793
794   // 9 multiplicity bins
795   // 1st.     <5    TPC tracks       fill 1  Bin1 [0.5,1.5)
796   // 2nd. >= 5  <10 TPC tracks       fill 2  Bin2 [1.5,2.5)
797   // 3rd. >= 10 <15 TPC tracks       fill 3  Bin3 [2.5,3.5)
798   // 4rd. >= 15 <20 TPC tracks       fill 4  Bin4 [3.5,4.5)
799   // 5th. >= 20 <30 TPC tracks       fill 5  Bin5 [4.5,5.5)
800   // 6th. >= 30 <40 TPC tracks       fill 6  Bin6 [5.5,6.5)
801   // 7th. >= 40 <50 TPC tracks       fill 7  Bin7 [6.5,7.5)
802   // 8th.    >50    TPC tracks       fill 8  Bin8 [7.5,8.5)
803
804   fNChTrRecECorrPPMult = new TH3F("fNChTrRecECorrPPMult","NTX in ener. corr. jet , corr. jet pT, pp mult.",101,-0.5,100.5, 60, 0., 300.,8,0.5,8.5);
805   fNChTrRecECorrPPMult->GetXaxis()->SetTitle("NTracks_Corrected");
806   fNChTrRecECorrPPMult->GetYaxis()->SetTitle("p_{T}^{JET}");
807   fNChTrRecECorrPPMult->GetZaxis()->SetTitle("Multiplicity Bin");
808   fNChTrRecECorrPPMult->Sumw2();
809   fOutputList->Add(fNChTrRecECorrPPMult);
810
811   fNChTrRecPerpECorrPPMult = new TH3F("fNChTrRecPerpECorrPPMult","Tracks above min in perp.cone , corr. jet pT, centrality",101,-0.5,100.5, 60, 0., 300.,8,0.5,8.5);
812   fNChTrRecPerpECorrPPMult->GetXaxis()->SetTitle("NTracks_Corrected");
813   fNChTrRecPerpECorrPPMult->GetYaxis()->SetTitle("p_{T}^{JET}");
814   fNChTrRecPerpECorrPPMult->GetZaxis()->SetTitle("Multiplicity Bin");
815   fNChTrRecPerpECorrPPMult->Sumw2();
816   fOutputList->Add(fNChTrRecPerpECorrPPMult);
817
818   fJetPtCentPbPbRaw = new TH2F("fJetPtCentPbPbRaw", "raw p_{T} distribution of reco jets", 60, 0., 300.,10,0.5,10.5);
819   fJetPtCentPbPbRaw->GetXaxis()->SetTitle("p_{T} (GeV/c)");
820   fJetPtCentPbPbRaw->GetYaxis()->SetTitle("Selection Bin");
821   fJetPtCentPbPbRaw->Sumw2();
822   fOutputList->Add(fJetPtCentPbPbRaw);
823
824   fJetPtCentPbPbCorr = new TH2F("fJetPtCentPbPbCorr", "Corrected p_{T} distribution of reco jets", 60, 0., 300.,10,0.5,10.5);
825   fJetPtCentPbPbCorr->GetXaxis()->SetTitle("p_{T} (GeV/c)");
826   fJetPtCentPbPbCorr->GetYaxis()->SetTitle("Selection Bin");
827   fJetPtCentPbPbCorr->Sumw2();
828   fOutputList->Add(fJetPtCentPbPbCorr);
829
830   fTotTracksCone = new TH3F("fTotTracksCone","Total number of tracks in the cone, raw jet pT bin, centrality", 200, 0.5, 200.5, 8, 0.5, 8.5, 10, 0.5, 10.5);
831   fTotTracksCone->GetXaxis()->SetTitle("NTracks");
832   fTotTracksCone->GetYaxis()->SetTitle("p_{T}^{JET} Bin");
833   fTotTracksCone->GetZaxis()->SetTitle("Centrality Bin");
834   fTotTracksCone->Sumw2();
835   fOutputList->Add(fTotTracksCone);
836
837   for(Int_t ipt=0;ipt<12;ipt++)
838     {
839       fNChTr[ipt] = new TH2F(Form("fNChTr[%i]",ipt),"Number of tracks to recover transverse energy, jet_{p_{T}}",101,-0.5,100.5, 60, 0., 300.);
840       fNChTr[ipt]->GetXaxis()->SetTitle("NTracks");
841       fNChTr[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
842       fNChTr[ipt]->Sumw2();
843       fOutputList->Add(fNChTr[ipt]);
844
845       fHistPtParton[ipt] = new TH1F(Form("fHistPtParton[%i]",ipt),"pT distribution of jets",50,0.,250.);
846       fHistPtParton[ipt]->GetXaxis()->SetTitle("p_{T}^{JET}");
847       fHistPtParton[ipt]->GetYaxis()->SetTitle("Entries");
848       fHistPtParton[ipt]->Sumw2();
849       fOutputList->Add(fHistPtParton[ipt]);
850
851       fSCM[ipt] = new TH2F(Form("fSCM[%i]",ipt),"Second Central Moment, jet_{p_{T}}",200,0.,0.2, 60, 0., 300.);
852       fSCM[ipt]->GetXaxis()->SetTitle("<#delta R_{c}^{2}>");
853       fSCM[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
854       fSCM[ipt]->Sumw2();
855       fOutputList->Add(fSCM[ipt]);
856
857       if(ipt<8) 
858         { 
859           fNChTrRDMult[ipt] = new TH2F(Form("fNChTrRDMult[%i]",ipt),"Number of tracks to recover transverse energy, jet_{p_{T}}",101,-0.5,100.5, 60, 0., 300.);
860           fNChTrRDMult[ipt]->GetXaxis()->SetTitle("NTracks");
861           fNChTrRDMult[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
862           fNChTrRDMult[ipt]->Sumw2();
863           fOutputList->Add(fNChTrRDMult[ipt]);
864           
865           fNAccJetsRDMult[ipt] = new TH1F(Form("fNAccJetsRDMult[%i]",ipt),"Number of accepted jets per event in real data", 101, -0.5, 100.5);
866           fNAccJetsRDMult[ipt]->GetXaxis()->SetTitle("Number of jets");
867           fNAccJetsRDMult[ipt]->Sumw2();
868           fOutputList->Add(fNAccJetsRDMult[ipt]);
869           
870           fTotalJetCharge[ipt] = new TH1F(Form("fTotalJetCharge[%i]",ipt),"Charge in the jet", 41, -20.5, 20.5);
871           fTotalJetCharge[ipt]->GetXaxis()->SetTitle("Charge in jet");
872           fTotalJetCharge[ipt]->Sumw2();
873           fOutputList->Add(fTotalJetCharge[ipt]);
874           
875           fSCMRDMult[ipt] = new TH2F(Form("fSCMRDMult[%i]",ipt),"Second Central Moment, jet_{p_{T}}",200,0.,0.2, 60, 0., 300.);
876           fSCMRDMult[ipt]->GetXaxis()->SetTitle("<#delta R_{c}^{2}>");
877           fSCMRDMult[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
878           fSCMRDMult[ipt]->Sumw2();
879           fOutputList->Add(fSCMRDMult[ipt]);
880           
881           fNChTrRDMultMC[ipt] = new TH2F(Form("fNChTrRDMultMC[%i]",ipt),"Number of tracks to recover transverse energy, jet_{p_{T}}",101,-0.5,100.5, 60, 0., 300.);
882           fNChTrRDMultMC[ipt]->GetXaxis()->SetTitle("NTracks");
883           fNChTrRDMultMC[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
884           fNChTrRDMultMC[ipt]->Sumw2();
885           fOutputList->Add(fNChTrRDMultMC[ipt]);
886           
887           fSCMRDMultMC[ipt] = new TH2F(Form("fSCMRDMultMC[%i]",ipt),"Second Central Moment, jet_{p_{T}}",200,0.,0.2, 60, 0., 300.);
888           fSCMRDMultMC[ipt]->GetXaxis()->SetTitle("<#delta R_{c}^{2}>");
889           fSCMRDMultMC[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
890           fSCMRDMultMC[ipt]->Sumw2();
891           fOutputList->Add(fSCMRDMultMC[ipt]);
892
893           //Second multiplicity estimator, removing jets and an area
894           fNChTrRDMultSE[ipt] = new TH2F(Form("fNChTrRDMultSE[%i]",ipt),"Number of tracks to recover transverse energy, jet_{p_{T}}",101,-0.5,100.5, 60, 0., 300.);
895           fNChTrRDMultSE[ipt]->GetXaxis()->SetTitle("NTracks");
896           fNChTrRDMultSE[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
897           fNChTrRDMultSE[ipt]->Sumw2();
898           fOutputList->Add(fNChTrRDMultSE[ipt]);
899           
900           fNAccJetsRDMultSE[ipt] = new TH1F(Form("fNAccJetsRDMultSE[%i]",ipt),"Number of accepted jets per event in real data", 101, -0.5, 100.5);
901           fNAccJetsRDMultSE[ipt]->GetXaxis()->SetTitle("Number of jets");
902           fNAccJetsRDMultSE[ipt]->Sumw2();
903           fOutputList->Add(fNAccJetsRDMultSE[ipt]);
904           
905           fTotalJetChargeSE[ipt] = new TH1F(Form("fTotalJetChargeSE[%i]",ipt),"Charge in the jet", 41, -20.5, 20.5);
906           fTotalJetChargeSE[ipt]->GetXaxis()->SetTitle("Charge in jet");
907           fTotalJetChargeSE[ipt]->Sumw2();
908           fOutputList->Add(fTotalJetChargeSE[ipt]);
909           
910           fSCMRDMultSE[ipt] = new TH2F(Form("fSCMRDMultSE[%i]",ipt),"Second Central Moment, jet_{p_{T}}",200,0.,0.2, 60, 0., 300.);
911           fSCMRDMultSE[ipt]->GetXaxis()->SetTitle("<#delta R_{c}^{2}>");
912           fSCMRDMultSE[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
913           fSCMRDMultSE[ipt]->Sumw2();
914           fOutputList->Add(fSCMRDMultSE[ipt]);
915
916           fNChTrRDMultOJ[ipt] = new TH2F(Form("fNChTrRDMultOJ[%i]",ipt),"Number of tracks to recover transverse energy, jet_{p_{T}, 1 jet}",101,-0.5,100.5, 60, 0., 300.);
917           fNChTrRDMultOJ[ipt]->GetXaxis()->SetTitle("NTracks");
918           fNChTrRDMultOJ[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
919           fNChTrRDMultOJ[ipt]->Sumw2();
920           fOutputList->Add(fNChTrRDMultOJ[ipt]);
921
922           fSCMRDMultOJ[ipt] = new TH2F(Form("fSCMRDMultOJ[%i]",ipt),"Second Central Moment, jet_{p_{T}}, 1 jet",200,0.,0.2, 60, 0., 300.);
923           fSCMRDMultOJ[ipt]->GetXaxis()->SetTitle("<#delta R_{c}^{2}>");
924           fSCMRDMultOJ[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
925           fSCMRDMultOJ[ipt]->Sumw2();
926           fOutputList->Add(fSCMRDMultOJ[ipt]);
927
928           fNChTrRDMultSEOJ[ipt] = new TH2F(Form("fNChTrRDMultSEOJ[%i]",ipt),"Number of tracks to recover transverse energy, jet_{p_{T}}, 1 jet",101,-0.5,100.5, 60, 0., 300.);
929           fNChTrRDMultSEOJ[ipt]->GetXaxis()->SetTitle("NTracks");
930           fNChTrRDMultSEOJ[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
931           fNChTrRDMultSEOJ[ipt]->Sumw2();
932           fOutputList->Add(fNChTrRDMultSEOJ[ipt]);
933
934           fSCMRDMultSEOJ[ipt] = new TH2F(Form("fSCMRDMultSEOJ[%i]",ipt),"Second Central Moment, jet_{p_{T}}, 1 jet",200,0.,0.2, 60, 0., 300.);
935           fSCMRDMultSEOJ[ipt]->GetXaxis()->SetTitle("<#delta R_{c}^{2}>");
936           fSCMRDMultSEOJ[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
937           fSCMRDMultSEOJ[ipt]->Sumw2();
938           fOutputList->Add(fSCMRDMultSEOJ[ipt]);
939
940           fNChTrRDMultOJMC[ipt] = new TH2F(Form("fNChTrRDMultOJMC[%i]",ipt),"Number of tracks to recover transverse energy, jet_{p_{T}, 1 jet, MC}",101,-0.5,100.5, 60, 0., 300.);
941           fNChTrRDMultOJMC[ipt]->GetXaxis()->SetTitle("NTracks");
942           fNChTrRDMultOJMC[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
943           fNChTrRDMultOJMC[ipt]->Sumw2();
944           fOutputList->Add(fNChTrRDMultOJMC[ipt]);
945
946           fSCMRDMultOJMC[ipt] = new TH2F(Form("fSCMRDMultOJMC[%i]",ipt),"Second Central Moment, jet_{p_{T}}, 1 jet, MC",200,0.,0.2, 60, 0., 300.);
947           fSCMRDMultOJMC[ipt]->GetXaxis()->SetTitle("<#delta R_{c}^{2}>");
948           fSCMRDMultOJMC[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
949           fSCMRDMultOJMC[ipt]->Sumw2();
950           fOutputList->Add(fSCMRDMultOJMC[ipt]);
951
952           fNChTrRDMultSEOJMC[ipt] = new TH2F(Form("fNChTrRDMultSEOJMC[%i]",ipt),"Number of tracks to recover transverse energy, jet_{p_{T}}, 1 jet, MC",101,-0.5,100.5, 60, 0., 300.);
953           fNChTrRDMultSEOJMC[ipt]->GetXaxis()->SetTitle("NTracks");
954           fNChTrRDMultSEOJMC[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
955           fNChTrRDMultSEOJMC[ipt]->Sumw2();
956           fOutputList->Add(fNChTrRDMultSEOJMC[ipt]);
957
958           fSCMRDMultSEOJMC[ipt] = new TH2F(Form("fSCMRDMultSEOJMC[%i]",ipt),"Second Central Moment, jet_{p_{T}}, 1 jet, MC",200,0.,0.2, 60, 0., 300.);
959           fSCMRDMultSEOJMC[ipt]->GetXaxis()->SetTitle("<#delta R_{c}^{2}>");
960           fSCMRDMultSEOJMC[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
961           fSCMRDMultSEOJMC[ipt]->Sumw2();
962           fOutputList->Add(fSCMRDMultSEOJMC[ipt]);
963
964           fNChTrRecPerpMultSEOJ[ipt] = new TH2F(Form("fNChTrRecPerpMultSEOJ[%i]",ipt),"Number of tracks above the min pT used in NTX_Raw, jet_{p_{T}}, 1 jet",101,-0.5,100.5, 60, 0., 300.);
965           fNChTrRecPerpMultSEOJ[ipt]->GetXaxis()->SetTitle("NTracks_{Exc.}");
966           fNChTrRecPerpMultSEOJ[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
967           fNChTrRecPerpMultSEOJ[ipt]->Sumw2();
968           fOutputList->Add(fNChTrRecPerpMultSEOJ[ipt]);
969
970           fJEtaMCMultOJ[ipt] = new TH1F(Form("fJEtaMCMultOJ[%i]",ipt), "Eta distribution of MC jets, V0-like mult. bins, 1 jet", 50, -1.5, 1.5);
971           fJEtaMCMultOJ[ipt]->GetXaxis()->SetTitle("#eta");
972           fJEtaMCMultOJ[ipt]->GetYaxis()->SetTitle("entries");
973           fJEtaMCMultOJ[ipt]->Sumw2();
974           fOutputList->Add(fJEtaMCMultOJ[ipt]);
975
976           fJEtaMCMultSEOJ[ipt] = new TH1F(Form("fJEtaMCMultSEOJ[%i]",ipt), "Eta distribution of MC jets, TPC-like mult. bins, 1 jet", 50, -1.5, 1.5);
977           fJEtaMCMultSEOJ[ipt]->GetXaxis()->SetTitle("#eta");
978           fJEtaMCMultSEOJ[ipt]->GetYaxis()->SetTitle("entries");
979           fJEtaMCMultSEOJ[ipt]->Sumw2();
980           fOutputList->Add(fJEtaMCMultSEOJ[ipt]);
981
982           fJEtaRDMultOJ[ipt] = new TH1F(Form("fJEtaRDMultOJ[%i]",ipt), "Eta distribution of reco jets, V0 mult. bins, 1 jet", 50, -1.5, 1.5);
983           fJEtaRDMultOJ[ipt]->GetXaxis()->SetTitle("#eta");
984           fJEtaRDMultOJ[ipt]->GetYaxis()->SetTitle("entries");
985           fJEtaRDMultOJ[ipt]->Sumw2();
986           fOutputList->Add(fJEtaRDMultOJ[ipt]);
987
988           fJEtaRDMultSEOJ[ipt] = new TH1F(Form("fJEtaRDMultSEOJ[%i]",ipt), "Eta distribution of reco jets, TPC mult. bins, 1 jet", 50, -1.5, 1.5);
989           fJEtaRDMultSEOJ[ipt]->GetXaxis()->SetTitle("#eta");
990           fJEtaRDMultSEOJ[ipt]->GetYaxis()->SetTitle("entries");
991           fJEtaRDMultSEOJ[ipt]->Sumw2();
992           fOutputList->Add(fJEtaRDMultSEOJ[ipt]);
993
994           fJetPtMCMultOJ[ipt] = new TH1F(Form("fJetPtMCMultOJ[%i]",ipt), "pT distribution of MC jets, V0-like mult. bins, 1 jet", 60, 0., 300.);
995           fJetPtMCMultOJ[ipt]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
996           fJetPtMCMultOJ[ipt]->GetYaxis()->SetTitle("entries");
997           fJetPtMCMultOJ[ipt]->Sumw2();
998           fOutputList->Add(fJetPtMCMultOJ[ipt]);
999
1000           fJetPtMCMultSEOJ[ipt] = new TH1F(Form("fJetPtMCMultSEOJ[%i]",ipt), "pT distribution of MC jets, TPC-like mult. bins, 1 jet", 60, 0., 300.);
1001           fJetPtMCMultSEOJ[ipt]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1002           fJetPtMCMultSEOJ[ipt]->GetYaxis()->SetTitle("entries");
1003           fJetPtMCMultSEOJ[ipt]->Sumw2();
1004           fOutputList->Add(fJetPtMCMultSEOJ[ipt]);
1005
1006           fJetPtRDMultOJ[ipt] = new TH1F(Form("fJetPtRDMultOJ[%i]",ipt), "pT distribution of reco jets, V0 mult. bins, 1 jet", 60, 0., 300.);
1007           fJetPtRDMultOJ[ipt]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1008           fJetPtRDMultOJ[ipt]->GetYaxis()->SetTitle("entries");
1009           fJetPtRDMultOJ[ipt]->Sumw2();
1010           fOutputList->Add(fJetPtRDMultOJ[ipt]);
1011
1012           fJetPtRDMultSEOJ[ipt] = new TH1F(Form("fJetPtRDMultSEOJ[%i]",ipt), "pT distribution of reco jets, TPC mult. bins, 1 jet", 60, 0., 300.);
1013           fJetPtRDMultSEOJ[ipt]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1014           fJetPtRDMultSEOJ[ipt]->GetYaxis()->SetTitle("entries");
1015           fJetPtRDMultSEOJ[ipt]->Sumw2();
1016           fOutputList->Add(fJetPtRDMultSEOJ[ipt]);
1017
1018           fEntriesQuark[ipt] = new TH2F(Form("fEntriesQuark[%i]",ipt),"NTX of quarks in multiplicity bins, 1 jet, MC",101,-0.5,100.5, 60, 0., 300.);
1019           fEntriesQuark[ipt]->GetXaxis()->SetTitle("NTracks");
1020           fEntriesQuark[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
1021           fEntriesQuark[ipt]->Sumw2();
1022           fOutputList->Add(fEntriesQuark[ipt]);
1023
1024           fEntriesGluon[ipt] = new TH2F(Form("fEntriesGluon[%i]",ipt),"NTX of gluons in multiplicity bins, 1 jet, MC",101,-0.5,100.5, 60, 0., 300.);
1025           fEntriesGluon[ipt]->GetXaxis()->SetTitle("NTracks");
1026           fEntriesGluon[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
1027           fEntriesGluon[ipt]->Sumw2();
1028           fOutputList->Add(fEntriesGluon[ipt]);
1029         } // end if <8
1030
1031       if(ipt<6)  // only entries for reconstructed || simulated jets
1032         { 
1033           fProcessPDG[ipt] = new TH2I(Form("fProcessPDG[%i]",ipt),"Pythia process and jet p_T", 60, 0., 300., 70, 0.5, 70.5);
1034           fProcessPDG[ipt]->GetXaxis()->SetTitle("p_{T}^{JET}");
1035           fProcessPDG[ipt]->GetYaxis()->SetTitle("Pythia process code");
1036           fProcessPDG[ipt]->Sumw2();
1037           fOutputList->Add(fProcessPDG[ipt]);
1038
1039           fFragPion[ipt] = new TH2F(Form("fFragPion[%i]",ipt),"Fragmentation in pions",35,0.,8.,50,0.,250.);
1040           fFragPion[ipt]->GetXaxis()->SetTitle("#xi=ln[Jet_{E_{T}}/H_{p_{T}}]");
1041           fFragPion[ipt]->GetYaxis()->SetTitle("Jet_{E_{T}}");
1042           fFragPion[ipt]->Sumw2();
1043           fOutputList->Add(fFragPion[ipt]);
1044           
1045           fFragKaon[ipt] = new TH2F(Form("fFragKaon[%i]",ipt),"Fragmentation in kaons",35,0.,8.,50,0.,250.);
1046           fFragKaon[ipt]->GetXaxis()->SetTitle("#xi=ln[Jet_{E_{T}}/H_{p_{T}}]");
1047           fFragKaon[ipt]->GetYaxis()->SetTitle("Jet_{E_{T}}");
1048           fFragKaon[ipt]->Sumw2();
1049           fOutputList->Add(fFragKaon[ipt]);
1050
1051           fFragProton[ipt] = new TH2F(Form("fFragProton[%i]",ipt),"Fragmentation in protons",35,0.,8.,50,0.,250.);
1052           fFragProton[ipt]->GetXaxis()->SetTitle("#xi=ln[Jet_{E_{T}}/H_{p_{T}}]");
1053           fFragProton[ipt]->GetYaxis()->SetTitle("Jet_{E_{T}}");
1054           fFragProton[ipt]->Sumw2();
1055           fOutputList->Add(fFragProton[ipt]);
1056
1057           fFragChargedR4[ipt] = new TH2F(Form("fFragChargedR4[%i]",ipt),"Fragmentation in charged particles",35,0.,8.,50,0.,250.);
1058           fFragChargedR4[ipt]->GetXaxis()->SetTitle("#xi=ln[Jet_{E_{T}}/H_{p_{T}}]");
1059           fFragChargedR4[ipt]->GetYaxis()->SetTitle("Jet_{E_{T}}");
1060           fFragChargedR4[ipt]->Sumw2();
1061           fOutputList->Add(fFragChargedR4[ipt]);
1062           
1063           fFragChargedR3[ipt] = new TH2F(Form("fFragChargedR3[%i]",ipt),"Fragmentation in charged particles",35,0.,8.,50,0.,250.);
1064           fFragChargedR3[ipt]->GetXaxis()->SetTitle("#xi=ln[Jet_{E_{T}}/H_{p_{T}}]");
1065           fFragChargedR3[ipt]->GetYaxis()->SetTitle("Jet_{E_{T}}");
1066           fFragChargedR3[ipt]->Sumw2();
1067           fOutputList->Add(fFragChargedR3[ipt]);
1068
1069           fFragChargedR2[ipt] = new TH2F(Form("fFragChargedR2[%i]",ipt),"Fragmentation in charged particles",35,0.,8.,50,0.,250.);
1070           fFragChargedR2[ipt]->GetXaxis()->SetTitle("#xi=ln[Jet_{E_{T}}/H_{p_{T}}]");
1071           fFragChargedR2[ipt]->GetYaxis()->SetTitle("Jet_{E_{T}}");
1072           fFragChargedR2[ipt]->Sumw2();
1073           fOutputList->Add(fFragChargedR2[ipt]);
1074
1075           // do not add the temporary containers
1076           fHistContainerR4[ipt] = new TH2F(Form("fHistContainerR4[%i]",ipt),"Temporary fragmentation in charged particles",35,0.,8.,50,0.,250.);
1077           fHistContainerR4[ipt]->GetXaxis()->SetTitle("#xi=ln[Jet_{E_{T}}/H_{p_{T}}]");
1078           fHistContainerR4[ipt]->GetYaxis()->SetTitle("Jet_{E_{T}}");
1079           fHistContainerR4[ipt]->Sumw2();
1080           
1081           fHistContainerR3[ipt] = new TH2F(Form("fHistContainerR3[%i]",ipt),"Temporary fragmentation in charged particles",35,0.,8.,50,0.,250.);
1082           fHistContainerR3[ipt]->GetXaxis()->SetTitle("#xi=ln[Jet_{E_{T}}/H_{p_{T}}]");
1083           fHistContainerR3[ipt]->GetYaxis()->SetTitle("Jet_{E_{T}}");
1084           fHistContainerR3[ipt]->Sumw2();
1085           
1086           fHistContainerR2[ipt] = new TH2F(Form("fHistContainerR2[%i]",ipt),"Temporary fragmentation in charged particles",35,0.,8.,50,0.,250.);
1087           fHistContainerR2[ipt]->GetXaxis()->SetTitle("#xi=ln[Jet_{E_{T}}/H_{p_{T}}]");
1088           fHistContainerR2[ipt]->GetYaxis()->SetTitle("Jet_{E_{T}}");
1089           fHistContainerR2[ipt]->Sumw2();
1090
1091           if(ipt<3)
1092             {
1093               fJetEtaJetPt[ipt] = new TH1F(Form("fJetEtaJetPt[%i]",ipt), "Eta distribution of reconstructed jets, all cut, with pT upper boundary", 50, -1.5, 1.5);
1094               fJetEtaJetPt[ipt]->GetXaxis()->SetTitle("#eta");
1095               fJetEtaJetPt[ipt]->GetYaxis()->SetTitle("entries");
1096               fJetEtaJetPt[ipt]->Sumw2();
1097               fOutputList->Add(fJetEtaJetPt[ipt]);
1098
1099               if(ipt<2)
1100                 {
1101                   fFragCandidates[ipt] = new TH2F(Form("fFragCandidates[%i]",ipt),"Parton identified candidates",35,0.,8.,50,0.,250.);
1102                   fFragCandidates[ipt]->GetXaxis()->SetTitle("#xi=ln[Jet_{E_{T}}/H_{p_{T}}]");
1103                   fFragCandidates[ipt]->GetYaxis()->SetTitle("Jet_{E_{T}}");
1104                   fFragCandidates[ipt]->Sumw2();
1105                   fOutputList->Add(fFragCandidates[ipt]);
1106                   
1107                   fMinTrackPtInNTXh[ipt] = new TH3F(Form("fMinTrackPtInNTXh[%i]",ipt), "Minimum track pT for the jets", 200, 0., 100., 60, 0., 300.,10,0.5,10.5);
1108                   fMinTrackPtInNTXh[ipt]->GetXaxis()->SetTitle("p_{T}^{TRACK}");
1109                   fMinTrackPtInNTXh[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
1110                   fMinTrackPtInNTXh[ipt]->GetZaxis()->SetTitle("Selection Bin"); //9 selection bins as fNChTrRecECorr
1111                   fMinTrackPtInNTXh[ipt]->Sumw2();
1112                   fOutputList->Add(fMinTrackPtInNTXh[ipt]);
1113                   
1114                   fMaxTrackPtInNTXh[ipt] = new TH2F(Form("fMaxTrackPtInNTXh[%i]",ipt), "Maximum track pT for the jets", 200, 0., 100., 60, 0., 300.);
1115                   fMaxTrackPtInNTXh[ipt]->GetXaxis()->SetTitle("p_{T}^{TRACK}");
1116                   fMaxTrackPtInNTXh[ipt]->GetYaxis()->SetTitle("p_{T}^{JET}");
1117                   fMaxTrackPtInNTXh[ipt]->Sumw2();
1118                   fOutputList->Add(fMaxTrackPtInNTXh[ipt]);  
1119                 } // index < 2 
1120             } // index < 3 
1121         } // index < 6 
1122     } // index < 12 
1123
1124   fPerpCone = new AliAODJet();
1125   fBckgSbsJet = new Double_t[3];
1126
1127   PostData(1, fOutputList); // Post data for ALL output slots >0 here, to get at least an empty histogram
1128 }
1129
1130 //________________________________________________________________________
1131 void AliAnalysisTaskPartonDisc::UserExec(Option_t *) 
1132 {
1133   // Main loop
1134   // Called for each event
1135
1136   //  fAOD
1137   if(fUseAODJetInput)
1138     {    
1139       fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
1140       if(!fAOD)
1141         {
1142           Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
1143           return;
1144         }
1145       // fetch the header
1146     }
1147   else
1148     {
1149       //  assume that the AOD is in the general output...
1150       fAOD  = AODEvent();
1151       if(!fAOD)
1152         {
1153           Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
1154           return;
1155         }
1156     }
1157
1158   // fin de test para fAOD
1159
1160   if(!fInputEvent)
1161     {
1162       Error("UserExec","No event found!");
1163       return;
1164     }
1165    
1166   AliAODHandler *aodHandler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1167   if(!aodHandler)
1168     {
1169       AliError("No AOD Handler!");
1170       return;
1171     }
1172
1173   fEventCent=900; //dummy val for debugging
1174
1175   if(fIsHIevent)
1176     {
1177       AliAODHeader *aodHeader = fAOD->GetHeader();
1178       fEventCent = aodHeader->GetCentrality();
1179     }
1180
1181   // Random number generator seeded per event
1182   TDatime dt;
1183   UInt_t curtime = dt.Get();
1184   UInt_t procid = gSystem->GetPid();
1185   UInt_t seed = curtime-procid;
1186   gRandom->SetSeed(seed);
1187
1188   // Jet eta exclusion
1189   if(fIncreasingExcl)
1190     fJetAcceptance = 0.5 - fIncExcR; // if the increase is 0.1 -> only jets within |eta|<0.4 
1191
1192   // First test of reference multiplicity
1193   Int_t refMultiplicity = fAOD->GetHeader()->GetRefMultiplicity();
1194   fRefMult->Fill(refMultiplicity);
1195
1196   // Multiplicity from V0 (V0A+V0C)
1197   fVZero = fAOD->GetVZEROData();
1198   Float_t multV0A = 0.0;
1199   Float_t multV0C = 0.0;
1200   Float_t multFullV0 = 0.0;
1201   if(fVZero)
1202     {
1203       multV0A = fVZero->GetMTotV0A();
1204       multV0C = fVZero->GetMTotV0C();
1205       multFullV0 = multV0A+multV0C; 
1206     }
1207   fVZEROMult->Fill(multV0A,multV0C);
1208
1209   fEvtCount++;
1210   Double_t jfr = fJetRadius;   // radius used during jet finding
1211   Int_t ntx = fNtX;   // NTX value
1212   const Int_t maxJetNum=6; // maximum number of generated jets to process
1213   AliAODJet genJets[6];  // containers for the
1214   AliAODJet recJets[6];  // correlation gen-reco
1215   Int_t nGenJets=0; 
1216   Int_t nRecJets=0;
1217   Int_t genJetsFlavor[6]={0};    // flavor of the generated jets
1218   Int_t evtype = 0; //pythia event type
1219   // Variables para la variable de estructura
1220   Double_t deltaPhiPt = 0.0;
1221   Double_t deltaEtaPt = 0.0;
1222   Double_t deltaPhiSqPt = 0.0;
1223   Double_t deltaEtaSqPt = 0.0;
1224   Double_t totalTrackPt = 0.0; 
1225   Double_t firstMomDeltPhi = 0.0;
1226   Double_t firstMomDeltEta = 0.0;
1227   Double_t secondMomDeltPhi = 0.0;
1228   Double_t secondMomDeltEta = 0.0;
1229   Double_t secondCentralPhi = 0.0;
1230   Double_t secondCentralEta = 0.0;
1231   Double_t secondCentralR = 0.0; 
1232
1233   // Variables para la variable de estructura
1234   // del cono perpendicular
1235   Double_t deltaPhiPtPerp = 0.0;
1236   Double_t deltaEtaPtPerp = 0.0;
1237   Double_t deltaPhiSqPtPerp = 0.0;
1238   Double_t deltaEtaSqPtPerp = 0.0;
1239   Double_t totalTrackPtPerp = 0.0;
1240   Double_t firstMomDeltPhiPerp = 0.0;
1241   Double_t firstMomDeltEtaPerp = 0.0;
1242   Double_t secondMomDeltPhiPerp = 0.0;
1243   Double_t secondMomDeltEtaPerp = 0.0;
1244   Double_t secondCentralPhiPerp = 0.0;
1245   Double_t secondCentralEtaPerp = 0.0;
1246   Double_t secondCentralRPerp = 0.0;  
1247
1248   Double_t perpendicularPt;
1249   Float_t px,py,pz,en; // jet 4-vector à la UA1
1250   Float_t pTbs, etabs, phibs; // energy corrected jet properties
1251
1252   // Process the MC info from the AOD
1253   if(fUseAODMC)
1254     {
1255       // Get the MC array
1256       TClonesArray *mcarray = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1257       if(!mcarray)
1258         {
1259           AliError("ERROR:No MC info in the AOD input");
1260           return;
1261         } 
1262
1263       AliMCEvent* mcEvent = MCEvent();
1264       if(mcEvent)
1265         {
1266           if(!fPhojetMC) // if it is pythia
1267             evtype = GetMCEventType(mcEvent);
1268           // From pythia 6.2 manual pp 414
1269           // QCD Hard Processes
1270           // 11 f_{i}+f_{j} -> f_{i}+f_{j} com77, ben84, eic84, chi90 
1271           // 12 f_{i}+barf_{i} -> f_{k}+barf_{k}
1272           // 13 f_{i}+barf_{i} -> g+g
1273           // 28 f_{i}+g -> f_{i}+g
1274           // 53 g+g -> f_{k}+barf_{k}
1275           // 68 g+g -> g+g
1276           if(fPhojetMC) // if it is phojet
1277             evtype = GetPhojetEventType(mcEvent); 
1278         }
1279       if(!mcEvent) // if a pure AOD event
1280         {
1281           AliDebug(2,Form("%s:%d No MCEvent \n",(char*)__FILE__,__LINE__));  
1282           AliDebug(2,Form("Trying to get the MC header \n"));  
1283           AliAODMCHeader *genEvH = static_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader"));
1284           if(!genEvH)
1285             {
1286               AliDebug(2,Form(" %s:%d No Pythia header!",(char*)__FILE__,__LINE__));  
1287               evtype = 0;
1288             }
1289           if(genEvH)
1290             evtype = genEvH->GetEventType();
1291         }   
1292       // Get the branch with the MC jets
1293       TClonesArray *aodMCJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchMC.Data()));
1294       if(!aodMCJets)
1295         {
1296           AliDebug(2,Form("%s:%d no MC Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchMC.Data()));  
1297           return;
1298         }
1299       AliDebug(2,Form("There are %d MC jets in this event\n", aodMCJets->GetEntries())); 
1300       Int_t mcjets =  aodMCJets->GetEntries();
1301       fNJetsMC->Fill(mcjets,mcjets); // number of jets FROM the branch, filled per event, this gives the event normalization...     
1302       HasOverlapedCones(aodMCJets); // Procedure for tagging usable jets
1303                                     // Up 16 jets are flagged
1304
1305       // Loop over AODMC physical primary charged particles
1306       // for the complete event
1307       Int_t tracksMC = mcarray->GetEntriesFast();         
1308       AliDebug(2,Form("There are %i tracks in the mcarray",tracksMC));  
1309       Double_t aodMCTrackEta = 0.0; 
1310       perpendicularPt = 0.0;
1311       px=0.0;
1312       py=0.0;
1313       pz=0.0;
1314       en=0.0;
1315       pTbs=0.0;
1316       etabs=0.0;
1317       phibs=0.0;
1318       fBckgSbsJet[0]=0.0;
1319       fBckgSbsJet[1]=0.0;
1320       fBckgSbsJet[2]=0.0;
1321       Int_t softRefMcNoJets = 0;
1322       Int_t myTotalMultiplicityMc = 0;
1323       Int_t v0LikeTotalMcMult = 0;
1324       Double_t randomNum = 0.;
1325       for(Int_t aodMCTrack = 0; aodMCTrack < tracksMC; aodMCTrack++ )
1326         {
1327           AliAODMCParticle *mctrackf = (AliAODMCParticle*) mcarray->At(aodMCTrack);
1328           if(!mctrackf) continue;
1329           if(!mctrackf->IsPhysicalPrimary()) continue;
1330           if(mctrackf->Charge()==0||mctrackf->Charge()==-99) continue;
1331           randomNum = gRandom->Rndm();
1332           if(randomNum<fTTrackRandomRejection) continue; //rechaza fTTrackRandomRejection
1333
1334           //Lo del V0, voy a contar particulas primarias cargadas
1335           if(mctrackf->Pt()>fMinpTValMC) // cut off en MC para MIPs
1336             {
1337               //V0A
1338               if(((mctrackf->Eta())>(2.8))&&((mctrackf->Eta())<(5.1)))
1339                 v0LikeTotalMcMult++;
1340               //V0C
1341               if(((mctrackf->Eta())>(-3.7))&&((mctrackf->Eta())<(-1.7)))
1342                 v0LikeTotalMcMult++;
1343             }
1344           //Fin de lo del V0
1345           aodMCTrackEta = TMath::Abs(mctrackf->Eta());
1346           if(aodMCTrackEta>0.9) continue;
1347           fPtAODMC->Fill(mctrackf->Pt(),mctrackf->Pt());
1348           fEtaAODMC->Fill(mctrackf->Eta(),mctrackf->Eta());
1349           fPhiAODMC->Fill(mctrackf->Phi(),mctrackf->Phi());
1350           if(fJetEvent) // if has an accepted jet, calculate the perpendicular cone
1351             {
1352               if(HasPerpendicularCone()) // If there is a perpendicular cone available
1353                 {
1354                   if(mctrackf->Pt()>fMinpTVal)
1355                     {
1356                     if(GetDeltaR(fEtaPerpCoord,fPhiPerpCoord,mctrackf->Eta(),mctrackf->Phi())<fJetRadius)
1357                       perpendicularPt = perpendicularPt + mctrackf->Pt();
1358                     }
1359                 }
1360             } // end IF jet event
1361           if(mctrackf->Pt()>fMinPtInGlobMult) // Min pT used in multiplicity estimation
1362             {
1363               myTotalMultiplicityMc++; // total multiplicity TPC like
1364               if(mctrackf->Pt()<fMinpTValUE) continue; // pT cut  fMinpTValUE
1365               if(mctrackf->Pt()>fMaxpTValUE) continue; // pT cut  fMaxpTValUE
1366               if(!IsTrackInsideExcludedArea(mctrackf->Eta(), mctrackf->Phi(), aodMCJets)) 
1367                 softRefMcNoJets++;
1368             }
1369         } // end loop over particles
1370
1371       Int_t correctedV0LikeMult= v0LikeTotalMcMult-GetV0LikeExcludedMultMC(aodMCJets,mcarray);
1372
1373       //estimadores
1374       if(mcjets==1) // correlation for only monojet events
1375         {
1376           fFullV0V0CorrUJMC->Fill(v0LikeTotalMcMult,correctedV0LikeMult);
1377           fTrackCountWOJetUJMC->Fill(myTotalMultiplicityMc,softRefMcNoJets); 
1378         }
1379
1380       if(fJetEvent) // if has an accepted jet, calculate the perpendicular cone
1381         {
1382           if(HasPerpendicularCone()) // If there is a perpendicular cone available
1383             {
1384               px = perpendicularPt*TMath::Cos(fPhiPerpCoord);
1385               py = perpendicularPt*TMath::Sin(fPhiPerpCoord);
1386               pz = perpendicularPt/TMath::Tan(2.0*TMath::ATan(TMath::Exp(-fEtaPerpCoord)));
1387               en = TMath::Sqrt(px*px + py*py + pz*pz);
1388               fPerpCone->SetPxPyPzE(px, py, pz, en);
1389             }
1390           if(!HasPerpendicularCone())
1391             AliDebug(2,"No perpendicular cone!!!");  
1392         }
1393
1394
1395       fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
1396  
1397       Int_t flavor = 0;     // flavor of the jet      
1398       Int_t nTracksPerc;    // ntx for the original jet
1399       Int_t nTracksPercBckgSubst;  // ntx for the energy corrected jet
1400       Double_t jetPt=0;
1401       Int_t pdgOfMCt;
1402       Float_t trackxi;
1403       Double_t jetXt;
1404       Double_t jetPts[7]={0};  // to store the pt of the jets
1405       Int_t mcJetCounter=0;    // counter of MC jets
1406       Int_t nTracksAboveThresholdPerp=0;  // n tracks of the perpendicular cone
1407       Int_t nTrUpThrPerpBckSubs=0;  // n tracks of the perpendicular cone, after the minimum pT recalculation
1408       fIsPossibleToSubstBckg = kTRUE; // Initialize before the loop
1409       if(fJetEvent) // si tiene jets validos
1410         {
1411           if(!HasPerpendicularCone()) // pero no encontro un cono perpendicular libre
1412             fIsPossibleToSubstBckg = kFALSE; // if not perpendicular cone, set to kFALSE, so no perpendicular calculations available
1413         }
1414       // Loop to fill a pT spectrum of the mc jets
1415       Int_t imcj=0; // index for montecarlo jets to correlate
1416       for (Int_t indxmc = 0; indxmc < mcjets; indxmc++) 
1417         {
1418           AliAODJet *mcjet = dynamic_cast<AliAODJet*>(aodMCJets->At(indxmc));
1419           if (!mcjet) 
1420             {
1421               AliDebug(2,Form("ERROR: Could not receive jet %d\n", indxmc));  
1422               continue;
1423             }
1424           
1425           ///////////////////////////////////////////////////////////////////////////////
1426           ///// Part for Chritians plot of inclusive and leading jets comp at 2.76 TeV //
1427           if(!IsInsideAcceptance(mcjet))  // old condition
1428             continue;
1429           if(indxmc==0) // leading jet
1430             fMCJetPtLeading->Fill(mcjet->Pt());
1431           fMCJetPtInclusive->Fill(mcjet->Pt()); // all
1432           ///// End of Christians Plot MC
1433           ///////////////////////////////////////////////////////////////////////////////
1434           
1435           if(indxmc>15)
1436             continue;
1437
1438           if(!fJetFlags[indxmc]) // If the jet is flaged kFALSE, not usable
1439             continue;
1440
1441           //Initialize variables for this jet
1442           //adiciones para la variable de estructura
1443           nTracksPerc = 0;
1444           nTracksPercBckgSubst = 0;
1445           fMinTrackPtInNTX=200.0;  //Initialize for each jet, overflown
1446           fMaxTrackPtInNTX=200.0;  //Initialize for each jet, overflown
1447           fMinTrackPtInNTXR=200.0;  //Initialize for each jet, overflown
1448           fMaxTrackPtInNTXR=200.0;  //Initialize for each jet, overflown
1449           deltaPhiPt = 0.0;
1450           deltaEtaPt = 0.0;
1451           deltaPhiSqPt = 0.0;
1452           deltaEtaSqPt = 0.0;
1453           totalTrackPt = 0.0;
1454           firstMomDeltPhi = 0.0;
1455           firstMomDeltEta = 0.0;
1456           secondMomDeltPhi = 0.0;
1457           secondMomDeltEta = 0.0;
1458           secondCentralPhi = 0.0;
1459           secondCentralEta = 0.0;
1460           secondCentralR = 0.0; 
1461           
1462           if(imcj<maxJetNum)
1463             genJets[imcj]= *mcjet;
1464           if(mcJetCounter<maxJetNum)
1465             jetPts[mcJetCounter]=mcjet->Pt();
1466           mcJetCounter++;  // number of jets in the acceptance
1467           jetPt = mcjet->Pt();
1468           flavor =  GetJetFlavour(mcjet,tracksMC,mcarray);
1469           if(imcj<maxJetNum)
1470             genJetsFlavor[imcj] = flavor;
1471           fJetPtMC->Fill(mcjet->Pt());
1472           fJetEtaMC->Fill(mcjet->Eta(),mcjet->Eta());
1473           fJetPhiMC->Fill(mcjet->Phi(),mcjet->Phi());
1474           fFlavor->Fill(flavor,jetPt);
1475           AliDebug(4,Form("Sabor del jet con pt=%f es :%d \n",jetPt,flavor)); 
1476           nTracksPerc = GetNumberOfMcChargedTracks(ntx,mcjet,tracksMC,mcarray,jfr); // este fija el min track pT, si es posible substraer el bckg
1477           if(fIsPossibleToSubstBckg&&!IsEqualRel(fCurrentJetMinPtNT90, 7000.)) //calculating only if there is a perpendicular cone available //IsEqualRel(jetpT, 0.0) //fCurrentJetMinPtNT90!=7000.
1478             {                                                     //and only if the method worked
1479               AliDebug(4,Form("For this jet and I have a perpendicular cone available")); 
1480               // Aqui lo que debo contar es el numero de tracks arriba del min pT del jet correspondiente
1481               // que es fCurrentJetMinPtNT90
1482               nTracksAboveThresholdPerp = GetNMcChargedTracksAboveThreshold(fPerpCone,tracksMC,mcarray,jfr);
1483             }
1484           // Corrected jet (pT)
1485           if(fIsPossibleToSubstBckg)   // for the current jet 
1486             {                                
1487               pTbs= mcjet->Pt()-fPerpCone->Pt();
1488               etabs= mcjet->Eta();
1489               phibs= mcjet->Phi();
1490               fBckgSbsJet[0]=pTbs; //pT
1491               fBckgSbsJet[1]=etabs; //eta
1492               fBckgSbsJet[2]=phibs; //phi
1493               // Now re-calculate nt90 for the energy corrected jet
1494               nTracksPercBckgSubst = GetRecalcNTXMc(ntx,mcjet,tracksMC,mcarray,jfr);
1495               // Now re-calculate the perpendicular cone NT90 background
1496               if(!IsEqualRel(fCurrentJetMinPtNT90Recalc, 7000.)) //calculating for the current jet, if the previos method worked //!IsEqualRel(fCurrentJetMinPtNT90, 7000.) //fCurrentJetMinPtNT90Recalc!=7000.
1497                 {
1498                   // Aqui lo que debo contar es el numero de tracks arriba del min pT del jet correspondiente
1499                   // despues de la correccion de energia del jet
1500                   nTrUpThrPerpBckSubs = GetRecalcNMcChTrUpThr(fPerpCone,tracksMC,mcarray,jfr);
1501                 }
1502             }
1503
1504           //check cross sections incoming partons
1505           jetXt= 2*jetPt/fSqrts;
1506           if(evtype==11||evtype==12||evtype==13) //QQ
1507             fFracQQ->Fill(jetXt);
1508           if(evtype==28) //GQ
1509             fFracGQ->Fill(jetXt);
1510           if(evtype==53||evtype==68) //GG
1511             fFracGG->Fill(jetXt);
1512
1513           //check cross sections outgoing partons
1514           if(evtype==11||evtype==12||evtype==53) //QQ
1515             fFracOutGoingQQ->Fill(jetXt);
1516           if(evtype==28) //GQ
1517             fFracOutGoingGQ->Fill(jetXt);
1518           if(evtype==13||evtype==68) //GG
1519             fFracOutGoingGG->Fill(jetXt);
1520
1521           fProcessJetPt->Fill(evtype,jetPt);  // pythia process, filled for each jet in acceptance
1522
1523           //Fill jet flavor as a function of pT and the pythia process but only leading jet
1524           if(imcj==0)  //leading jet
1525             {
1526               fFlavorLead->Fill(flavor,jetPt);
1527               fProcessLeadJetPt->Fill(evtype,jetPt); 
1528             }
1529           AliDebug(4,Form("Before the check of comparison")); 
1530           // To check tracks related to this MC jet
1531           // RefTracks check
1532           Bool_t rTrkFlagMC = kFALSE;
1533           Int_t trkinmcjet = mcjet->GetRefTracks()->GetEntriesFast();
1534           if(trkinmcjet!=0&&!fForceNotTR)
1535             rTrkFlagMC = kTRUE;
1536           AliDebug(4,Form("Number of tracks in RefTracks MC jet:%i \n",trkinmcjet));  
1537           if(rTrkFlagMC)  // If there are tracks refs available
1538             {
1539               AliDebug(4,Form("Checking composition in MC with track refs")); 
1540               for(Int_t aodMCT = 0; aodMCT < trkinmcjet; aodMCT++ )
1541                 {
1542                   pdgOfMCt=0;
1543                   trackxi=0;
1544                   AliAODMCParticle *mctrack = (AliAODMCParticle*) mcjet->GetRefTracks()->At(aodMCT);
1545                   if(!mctrack) continue;
1546                   if(!mctrack->IsPhysicalPrimary()) continue;
1547                   if(mctrack->Charge()==0||mctrack->Charge()==-99) continue;
1548                   if(mctrack->Pt()<fMinpTVal) continue; // MC no cut in the case of track reference, should be in, NO, cut anyhow to be safe
1549                   randomNum = gRandom->Rndm();
1550                   if(randomNum<fJTrackRandomRejection) continue; //rechaza fJTrackRandomRejection
1551                   deltaPhiPt += DeltaPhiMC(mcjet, mctrack)*mctrack->Pt();
1552                   deltaEtaPt += DeltaEtaMC(mcjet, mctrack)*mctrack->Pt();
1553                   deltaPhiSqPt += DeltaPhiSqMC(mcjet, mctrack)*mctrack->Pt();
1554                   deltaEtaSqPt += DeltaEtaSqMC(mcjet, mctrack)*mctrack->Pt();
1555                   totalTrackPt += mctrack->Pt();
1556
1557                   pdgOfMCt=abs(mctrack->GetPdgCode());
1558                   if(!IsEqualRel(mctrack->Pt(), 0.0)) //!IsEqualRel(fCurrentJetMinPtNT90Recalc, 7000.) // mctrack->Pt()!=0
1559                     trackxi= log(jetPt/mctrack->Pt());
1560                   switch(abs(flavor))
1561                     {
1562                     case 1:
1563                       if(pdgOfMCt==321)
1564                         fFragKaon[0]->Fill(trackxi,jetPt);
1565                       if(pdgOfMCt==211)
1566                         fFragPion[0]->Fill(trackxi,jetPt);
1567                       if(pdgOfMCt==2212)
1568                         fFragProton[0]->Fill(trackxi,jetPt);
1569                       break;
1570                     case 2:
1571                       if(pdgOfMCt==321)
1572                         fFragKaon[1]->Fill(trackxi,jetPt);
1573                       if(pdgOfMCt==211)
1574                         fFragPion[1]->Fill(trackxi,jetPt);
1575                       if(pdgOfMCt==2212)
1576                         fFragProton[1]->Fill(trackxi,jetPt);
1577                       break;
1578                     case 3:
1579                       if(pdgOfMCt==321)
1580                         fFragKaon[2]->Fill(trackxi,jetPt);
1581                       if(pdgOfMCt==211)
1582                         fFragPion[2]->Fill(trackxi,jetPt);
1583                       if(pdgOfMCt==2212)
1584                         fFragProton[2]->Fill(trackxi,jetPt);
1585                       break;
1586                     case 4:
1587                       if(pdgOfMCt==321)
1588                         fFragKaon[3]->Fill(trackxi,jetPt);
1589                       if(pdgOfMCt==211)
1590                         fFragPion[3]->Fill(trackxi,jetPt);
1591                       if(pdgOfMCt==2212)
1592                         fFragProton[3]->Fill(trackxi,jetPt);
1593                       break;
1594                     case 5:
1595                       if(pdgOfMCt==321)
1596                         fFragKaon[4]->Fill(trackxi,jetPt);
1597                       if(pdgOfMCt==211)
1598                         fFragPion[4]->Fill(trackxi,jetPt);
1599                       if(pdgOfMCt==2212)
1600                         fFragProton[4]->Fill(trackxi,jetPt);
1601                       break;
1602                     case 21:
1603                       if(pdgOfMCt==321)
1604                         fFragKaon[5]->Fill(trackxi,jetPt);
1605                       if(pdgOfMCt==211)
1606                         fFragPion[5]->Fill(trackxi,jetPt);
1607                       if(pdgOfMCt==2212)
1608                         fFragProton[5]->Fill(trackxi,jetPt);
1609                       break;          
1610                     default:
1611                       break;
1612                     } // end switch flavor
1613                 }// end loop over AODMC particles of trackrefs
1614               if(!IsEqualRel(totalTrackPt, 0.0)) //!IsEqualRel(totalTrackPt, 0.0) //totalTrackPt!=0.0
1615                 {
1616                   firstMomDeltPhi = deltaPhiPt/totalTrackPt;
1617                   firstMomDeltEta = deltaEtaPt/totalTrackPt;
1618                   secondMomDeltPhi = deltaPhiSqPt/totalTrackPt;
1619                   secondMomDeltEta = deltaEtaSqPt/totalTrackPt;
1620                   secondCentralPhi = secondMomDeltPhi - firstMomDeltPhi*firstMomDeltPhi;
1621                   secondCentralEta = secondMomDeltEta - firstMomDeltEta*firstMomDeltEta;
1622                   secondCentralR = secondCentralPhi + secondCentralEta;
1623                 } // end if totalTrackPt!=0.0
1624               if(IsEqualRel(totalTrackPt, 0.0))  //totalTrackPt==0.0
1625                 secondCentralR = 10.0; //overflow
1626             }// end version with ref tracks (flag check)
1627             
1628           if(!rTrkFlagMC)  // No ref tracks available
1629             {
1630               AliDebug(4,Form("Checking composition in MC without track refs")); 
1631               for(Int_t aodMCT = 0; aodMCT < tracksMC; aodMCT++ )
1632                 {
1633                   pdgOfMCt=0;
1634                   trackxi=0;
1635                   AliAODMCParticle *mctrack = (AliAODMCParticle*) mcarray->At(aodMCT);
1636                   if(!mctrack) continue;
1637                   if(!mctrack->IsPhysicalPrimary()) continue;
1638                   if(mctrack->Charge()==0||mctrack->Charge()==-99) continue;
1639                   if(!IsMCTrackInsideThisJet(mctrack, mcjet, jfr)) continue;
1640                   if(mctrack->Pt()<fMinpTVal) continue; // MC: HERE PT CUT, NO TRACK REF
1641                   randomNum = gRandom->Rndm();
1642                   if(randomNum<fJTrackRandomRejection) continue; //rechaza fJTrackRandomRejection
1643                   deltaPhiPt += DeltaPhiMC(mcjet, mctrack)*mctrack->Pt();
1644                   deltaEtaPt += DeltaEtaMC(mcjet, mctrack)*mctrack->Pt();
1645                   deltaPhiSqPt += DeltaPhiSqMC(mcjet, mctrack)*mctrack->Pt();
1646                   deltaEtaSqPt += DeltaEtaSqMC(mcjet, mctrack)*mctrack->Pt();
1647                   totalTrackPt += mctrack->Pt();
1648
1649                   pdgOfMCt=abs(mctrack->GetPdgCode());
1650                   if(!IsEqualRel(mctrack->Pt(), 0.0)) //!IsEqualRel(mctrack->Pt(), 0.0) // mctrack->Pt()!=0
1651                     trackxi= log(jetPt/mctrack->Pt());
1652                   switch(flavor)
1653                     {
1654                     case 1:
1655                       if(pdgOfMCt==321)
1656                         fFragKaon[0]->Fill(trackxi,jetPt);
1657                       if(pdgOfMCt==211)
1658                         fFragPion[0]->Fill(trackxi,jetPt);
1659                       if(pdgOfMCt==2212)
1660                         fFragProton[0]->Fill(trackxi,jetPt);
1661                       break;
1662                     case 2:
1663                       if(pdgOfMCt==321)
1664                         fFragKaon[1]->Fill(trackxi,jetPt);
1665                       if(pdgOfMCt==211)
1666                         fFragPion[1]->Fill(trackxi,jetPt);
1667                       if(pdgOfMCt==2212)
1668                         fFragProton[1]->Fill(trackxi,jetPt);
1669                       break;
1670                     case 3:
1671                       if(pdgOfMCt==321)
1672                         fFragKaon[2]->Fill(trackxi,jetPt);
1673                       if(pdgOfMCt==211)
1674                         fFragPion[2]->Fill(trackxi,jetPt);
1675                       if(pdgOfMCt==2212)
1676                         fFragProton[2]->Fill(trackxi,jetPt);
1677                       break;
1678                     case 4:
1679                       if(pdgOfMCt==321)
1680                         fFragKaon[3]->Fill(trackxi,jetPt);
1681                       if(pdgOfMCt==211)
1682                         fFragPion[3]->Fill(trackxi,jetPt);
1683                       if(pdgOfMCt==2212)
1684                         fFragProton[3]->Fill(trackxi,jetPt);
1685                       break;
1686                     case 5:
1687                       if(pdgOfMCt==321)
1688                         fFragKaon[4]->Fill(trackxi,jetPt);
1689                       if(pdgOfMCt==211)
1690                         fFragPion[4]->Fill(trackxi,jetPt);
1691                       if(pdgOfMCt==2212)
1692                         fFragProton[4]->Fill(trackxi,jetPt);
1693                       break;
1694                     case 21:
1695                       if(pdgOfMCt==321)
1696                         fFragKaon[5]->Fill(trackxi,jetPt);
1697                       if(pdgOfMCt==211)
1698                         fFragPion[5]->Fill(trackxi,jetPt);
1699                       if(pdgOfMCt==2212)
1700                     fFragProton[5]->Fill(trackxi,jetPt);
1701                       break;          
1702                     default:
1703                       break;
1704                     } // end switch flavor
1705                 }// end loop over AODMC particles
1706               if(!IsEqualRel(totalTrackPt, 0.0)) //!IsEqualRel(totalTrackPt, 0.0) // totalTrackPt!=0.0
1707                 {
1708                   firstMomDeltPhi = deltaPhiPt/totalTrackPt;
1709                   firstMomDeltEta = deltaEtaPt/totalTrackPt;
1710                   secondMomDeltPhi = deltaPhiSqPt/totalTrackPt;
1711                   secondMomDeltEta = deltaEtaSqPt/totalTrackPt;
1712                   secondCentralPhi = secondMomDeltPhi - firstMomDeltPhi*firstMomDeltPhi;
1713                   secondCentralEta = secondMomDeltEta - firstMomDeltEta*firstMomDeltEta;
1714                   secondCentralR = secondCentralPhi + secondCentralEta;
1715                 } // end if totalTrackPt!=0.0
1716               if(IsEqualRel(totalTrackPt, 0.0)) //!IsEqualRel(totalTrackPt, 0.0) //totalTrackPt==0.0
1717                 secondCentralR = 10.0; //overflow
1718             } //End old version (no ref tracks)
1719
1720           if(fIsPossibleToSubstBckg)
1721             {
1722               // To make sure, re-initialize
1723               deltaPhiPtPerp = 0.0;
1724               deltaEtaPtPerp = 0.0;
1725               deltaPhiSqPtPerp = 0.0;
1726               deltaEtaSqPtPerp = 0.0;
1727               totalTrackPtPerp = 0.0;
1728               firstMomDeltPhiPerp = 0.0;
1729               firstMomDeltEtaPerp = 0.0;
1730               secondMomDeltPhiPerp = 0.0;
1731               secondMomDeltEtaPerp = 0.0;
1732               secondCentralPhiPerp = 0.0;
1733               secondCentralEtaPerp = 0.0;
1734               secondCentralRPerp = 0.0;
1735
1736               AliDebug(4,Form("Checking SCM in MC for the perpendicular cone")); 
1737               for(Int_t aodMCperp = 0; aodMCperp < tracksMC; aodMCperp++ )
1738                 {
1739                   AliAODMCParticle *mctrackperp = (AliAODMCParticle*) mcarray->At(aodMCperp);
1740                   if(!mctrackperp) continue;
1741                   if(!mctrackperp->IsPhysicalPrimary()) continue;
1742                   if(mctrackperp->Charge()==0||mctrackperp->Charge()==-99) continue;
1743                   if(!IsMCTrackInsideThisJet(mctrackperp, fPerpCone, jfr)) continue;
1744                   if(mctrackperp->Pt()<fMinpTVal) continue; // MC: HERE PT CUT   
1745                   randomNum = gRandom->Rndm();
1746                   if(randomNum<fJTrackRandomRejection) continue; //rechaza fJTrackRandomRejection
1747                   deltaPhiPtPerp += DeltaPhiMC(fPerpCone, mctrackperp)*mctrackperp->Pt();
1748                   deltaEtaPtPerp += DeltaEtaMC(fPerpCone, mctrackperp)*mctrackperp->Pt();
1749                   deltaPhiSqPtPerp += DeltaPhiSqMC(fPerpCone, mctrackperp)*mctrackperp->Pt();
1750                   deltaEtaSqPtPerp += DeltaEtaSqMC(fPerpCone, mctrackperp)*mctrackperp->Pt();
1751                   totalTrackPtPerp += mctrackperp->Pt();
1752                 }// end loop over AODMC particles
1753               if(!IsEqualRel(totalTrackPtPerp, 0.0)) //!IsEqualRel(totalTrackPt, 0.0) // totalTrackPtPerp!=0.0
1754                 {
1755                   firstMomDeltPhiPerp = deltaPhiPtPerp/totalTrackPtPerp;
1756                   firstMomDeltEtaPerp = deltaEtaPtPerp/totalTrackPtPerp;
1757                   secondMomDeltPhiPerp = deltaPhiSqPtPerp/totalTrackPtPerp;
1758                   secondMomDeltEtaPerp = deltaEtaSqPtPerp/totalTrackPtPerp;
1759                   secondCentralPhiPerp = secondMomDeltPhiPerp - firstMomDeltPhiPerp*firstMomDeltPhiPerp;
1760                   secondCentralEtaPerp = secondMomDeltEtaPerp - firstMomDeltEtaPerp*firstMomDeltEtaPerp;
1761                   secondCentralRPerp = secondCentralPhiPerp + secondCentralEtaPerp;
1762                 } // end if totalTrackPt!=0.0
1763               if(IsEqualRel(totalTrackPtPerp, 0.0)) //!IsEqualRel(totalTrackPtPerp, 0.0) //totalTrackPtPerp==0.0
1764                 secondCentralRPerp = 10.0; //overflow
1765             }
1766           ///// end of adding the SCM for the perpendicular cone
1767
1768           if(mcjets==1) // if only one jet in the whole event, and inside acceptance
1769             {
1770               // reference multiplicity stuff in pp, also filled in PbPb, but does not matter.
1771               // set to: V0 like corrected multiplicity: correctedV0LikeMult
1772               if(correctedV0LikeMult<25)
1773                 {
1774                   fNChTrRDMultOJMC[0]->Fill(nTracksPerc,jetPt);
1775                   fSCMRDMultOJMC[0]->Fill(secondCentralR,jetPt);
1776                   fJEtaMCMultOJ[0]->Fill(mcjet->Eta());
1777                   fJetPtMCMultOJ[0]->Fill(mcjet->Pt());
1778                 }
1779               if(correctedV0LikeMult>=25&&correctedV0LikeMult<50)
1780                 {
1781                   fNChTrRDMultOJMC[1]->Fill(nTracksPerc,jetPt);
1782                   fSCMRDMultOJMC[1]->Fill(secondCentralR,jetPt);
1783                   fJEtaMCMultOJ[1]->Fill(mcjet->Eta());
1784                   fJetPtMCMultOJ[1]->Fill(mcjet->Pt());
1785                 }
1786               if(correctedV0LikeMult>=50&&correctedV0LikeMult<90)
1787                 {
1788                   fNChTrRDMultOJMC[2]->Fill(nTracksPerc,jetPt);
1789                   fSCMRDMultOJMC[2]->Fill(secondCentralR,jetPt);
1790                   fJEtaMCMultOJ[2]->Fill(mcjet->Eta());
1791                   fJetPtMCMultOJ[2]->Fill(mcjet->Pt());
1792                 }
1793               if(correctedV0LikeMult>=90&&correctedV0LikeMult<120)
1794                 {
1795                   fNChTrRDMultOJMC[3]->Fill(nTracksPerc,jetPt);
1796                   fSCMRDMultOJMC[3]->Fill(secondCentralR,jetPt);
1797                   fJEtaMCMultOJ[3]->Fill(mcjet->Eta());
1798                   fJetPtMCMultOJ[3]->Fill(mcjet->Pt());
1799                 }
1800               if(correctedV0LikeMult>=120&&correctedV0LikeMult<150)
1801                 {
1802                   fNChTrRDMultOJMC[4]->Fill(nTracksPerc,jetPt);
1803                   fSCMRDMultOJMC[4]->Fill(secondCentralR,jetPt);
1804                   fJEtaMCMultOJ[4]->Fill(mcjet->Eta());
1805                   fJetPtMCMultOJ[4]->Fill(mcjet->Pt());
1806                 }
1807               if(correctedV0LikeMult>=150&&correctedV0LikeMult<200)
1808                 {
1809                   fNChTrRDMultOJMC[5]->Fill(nTracksPerc,jetPt);
1810                   fSCMRDMultOJMC[5]->Fill(secondCentralR,jetPt);
1811                   fJEtaMCMultOJ[5]->Fill(mcjet->Eta());
1812                   fJetPtMCMultOJ[5]->Fill(mcjet->Pt());
1813                 }
1814               if(correctedV0LikeMult>=200&&correctedV0LikeMult<300)
1815                 {
1816                   fNChTrRDMultOJMC[6]->Fill(nTracksPerc,jetPt);
1817                   fSCMRDMultOJMC[6]->Fill(secondCentralR,jetPt);
1818                   fJEtaMCMultOJ[6]->Fill(mcjet->Eta());
1819                   fJetPtMCMultOJ[6]->Fill(mcjet->Pt());
1820                 }
1821               if(correctedV0LikeMult>=300)
1822                 {
1823                   fNChTrRDMultOJMC[7]->Fill(nTracksPerc,jetPt);
1824                   fSCMRDMultOJMC[7]->Fill(secondCentralR,jetPt);
1825                   fJEtaMCMultOJ[7]->Fill(mcjet->Eta());
1826                   fJetPtMCMultOJ[7]->Fill(mcjet->Pt());
1827                 }
1828               //Results for inclusive jets
1829               // 2nd. Reference: set to: TPC tracks minus jet, minus dijet area
1830               if(softRefMcNoJets<5)
1831                 {
1832                   fNChTrRDMultSEOJMC[0]->Fill(nTracksPerc,jetPt);
1833                   fSCMRDMultSEOJMC[0]->Fill(secondCentralR,jetPt);
1834                   fJEtaMCMultSEOJ[0]->Fill(mcjet->Eta());
1835                   fJetPtMCMultSEOJ[0]->Fill(mcjet->Pt());
1836                   if(abs(flavor)==1||abs(flavor)==2||abs(flavor)==3||abs(flavor)==4||abs(flavor)==5) //if quark jet
1837                     fEntriesQuark[0]->Fill(nTracksPerc,jetPt);
1838                   if(abs(flavor)==21) //if gluon jet
1839                     fEntriesGluon[0]->Fill(nTracksPerc,jetPt);
1840                 }
1841               if(softRefMcNoJets>=5&&softRefMcNoJets<10)
1842                 {
1843                   fNChTrRDMultSEOJMC[1]->Fill(nTracksPerc,jetPt);
1844                   fSCMRDMultSEOJMC[1]->Fill(secondCentralR,jetPt);
1845                   fJEtaMCMultSEOJ[1]->Fill(mcjet->Eta());
1846                   fJetPtMCMultSEOJ[1]->Fill(mcjet->Pt()); 
1847                   if(abs(flavor)==1||abs(flavor)==2||abs(flavor)==3||abs(flavor)==4||abs(flavor)==5) //if quark jet
1848                     fEntriesQuark[1]->Fill(nTracksPerc,jetPt);
1849                   if(abs(flavor)==21) //if gluon jet
1850                     fEntriesGluon[1]->Fill(nTracksPerc,jetPt);
1851                 }
1852               if(softRefMcNoJets>=10&&softRefMcNoJets<15)
1853                 {
1854                   fNChTrRDMultSEOJMC[2]->Fill(nTracksPerc,jetPt);
1855                   fSCMRDMultSEOJMC[2]->Fill(secondCentralR,jetPt);
1856                   fJEtaMCMultSEOJ[2]->Fill(mcjet->Eta());
1857                   fJetPtMCMultSEOJ[2]->Fill(mcjet->Pt());
1858                   if(abs(flavor)==1||abs(flavor)==2||abs(flavor)==3||abs(flavor)==4||abs(flavor)==5) //if quark jet
1859                     fEntriesQuark[2]->Fill(nTracksPerc,jetPt);
1860                   if(abs(flavor)==21) //if gluon jet
1861                     fEntriesGluon[2]->Fill(nTracksPerc,jetPt);
1862                 }
1863               if(softRefMcNoJets>=15&&softRefMcNoJets<20)
1864                 {
1865                   fNChTrRDMultSEOJMC[3]->Fill(nTracksPerc,jetPt);
1866                   fSCMRDMultSEOJMC[3]->Fill(secondCentralR,jetPt);
1867                   fJEtaMCMultSEOJ[3]->Fill(mcjet->Eta());
1868                   fJetPtMCMultSEOJ[3]->Fill(mcjet->Pt());
1869                   if(abs(flavor)==1||abs(flavor)==2||abs(flavor)==3||abs(flavor)==4||abs(flavor)==5) //if quark jet
1870                     fEntriesQuark[3]->Fill(nTracksPerc,jetPt);
1871                   if(abs(flavor)==21) //if gluon jet
1872                     fEntriesGluon[3]->Fill(nTracksPerc,jetPt);
1873                 }
1874               if(softRefMcNoJets>=20&&softRefMcNoJets<30)
1875                 {
1876                   fNChTrRDMultSEOJMC[4]->Fill(nTracksPerc,jetPt);
1877                   fSCMRDMultSEOJMC[4]->Fill(secondCentralR,jetPt);
1878                   fJEtaMCMultSEOJ[4]->Fill(mcjet->Eta());
1879                   fJetPtMCMultSEOJ[4]->Fill(mcjet->Pt());
1880                   if(abs(flavor)==1||abs(flavor)==2||abs(flavor)==3||abs(flavor)==4||abs(flavor)==5) //if quark jet
1881                     fEntriesQuark[4]->Fill(nTracksPerc,jetPt);
1882                   if(abs(flavor)==21) //if gluon jet
1883                     fEntriesGluon[4]->Fill(nTracksPerc,jetPt);
1884                 }
1885               if(softRefMcNoJets>=30&&softRefMcNoJets<40)
1886                 {
1887                   fNChTrRDMultSEOJMC[5]->Fill(nTracksPerc,jetPt);
1888                   fSCMRDMultSEOJMC[5]->Fill(secondCentralR,jetPt);
1889                   fJEtaMCMultSEOJ[5]->Fill(mcjet->Eta());
1890                   fJetPtMCMultSEOJ[5]->Fill(mcjet->Pt());
1891                   if(abs(flavor)==1||abs(flavor)==2||abs(flavor)==3||abs(flavor)==4||abs(flavor)==5) //if quark jet
1892                     fEntriesQuark[5]->Fill(nTracksPerc,jetPt);
1893                   if(abs(flavor)==21) //if gluon jet
1894                     fEntriesGluon[5]->Fill(nTracksPerc,jetPt);
1895                 }
1896               if(softRefMcNoJets>=40&&softRefMcNoJets<50)
1897                 {
1898                   fNChTrRDMultSEOJMC[6]->Fill(nTracksPerc,jetPt);
1899                   fSCMRDMultSEOJMC[6]->Fill(secondCentralR,jetPt);
1900                   fJEtaMCMultSEOJ[6]->Fill(mcjet->Eta());
1901                   fJetPtMCMultSEOJ[6]->Fill(mcjet->Pt());
1902                   if(abs(flavor)==1||abs(flavor)==2||abs(flavor)==3||abs(flavor)==4||abs(flavor)==5) //if quark jet
1903                     fEntriesQuark[6]->Fill(nTracksPerc,jetPt);
1904                   if(abs(flavor)==21) //if gluon jet
1905                     fEntriesGluon[6]->Fill(nTracksPerc,jetPt);
1906                 }
1907               if(softRefMcNoJets>=50)
1908                 {           
1909                   fNChTrRDMultSEOJMC[7]->Fill(nTracksPerc,jetPt);
1910                   fSCMRDMultSEOJMC[7]->Fill(secondCentralR,jetPt);
1911                   fJEtaMCMultSEOJ[7]->Fill(mcjet->Eta());
1912                   fJetPtMCMultSEOJ[7]->Fill(mcjet->Pt());
1913                   if(abs(flavor)==1||abs(flavor)==2||abs(flavor)==3||abs(flavor)==4||abs(flavor)==5) //if quark jet
1914                     fEntriesQuark[7]->Fill(nTracksPerc,jetPt);
1915                   if(abs(flavor)==21) //if gluon jet
1916                     fEntriesGluon[7]->Fill(nTracksPerc,jetPt);
1917                 }
1918             }
1919           //End results for inclusive jets,starts parton by parton
1920           
1921           switch(abs(flavor))
1922             {
1923             case 1:
1924               fNChTr[0]->Fill(nTracksPerc,jetPt);
1925               fProcessPDG[0]->Fill(jetPt,evtype);
1926               fHistPtParton[0]->Fill(jetPt);
1927               fSCM[0]->Fill(secondCentralR,jetPt);
1928               fMinTrackPtInNTXh[0]->Fill(fMinTrackPtInNTX,jetPt,1); // 0 for pp MC
1929               fMaxTrackPtInNTXh[0]->Fill(fMaxTrackPtInNTX,jetPt); // 0 for MC
1930               if(fIsPossibleToSubstBckg)
1931                 {
1932                   fNChTrCorrMCQuark->Fill(nTracksPercBckgSubst,pTbs); 
1933                   fSCMMCPerp->Fill(secondCentralR,jetPt);
1934                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.)) //!IsEqualRel(mctrack->Pt(), 0.0) // fCurrentJetMinPtNT90!=7000.
1935                     fNChTrMCPerp->Fill(nTracksAboveThresholdPerp,jetPt);
1936                   if(!IsEqualRel(fCurrentJetMinPtNT90Recalc, 7000.)) // !IsEqualRel(fCurrentJetMinPtNT90, 7000.) // fCurrentJetMinPtNT90Recalc!=7000.
1937                     fNChTrCorrMCPerp->Fill(nTrUpThrPerpBckSubs,pTbs);
1938                 }
1939               break;
1940             case 2:
1941               fNChTr[1]->Fill(nTracksPerc,jetPt);
1942               fProcessPDG[1]->Fill(jetPt,evtype);
1943               fHistPtParton[1]->Fill(jetPt);
1944               fSCM[1]->Fill(secondCentralR,jetPt);
1945               fMinTrackPtInNTXh[0]->Fill(fMinTrackPtInNTX,jetPt,1); // 0 for pp MC
1946               fMaxTrackPtInNTXh[0]->Fill(fMaxTrackPtInNTX,jetPt); // 0 for MC
1947               if(fIsPossibleToSubstBckg)
1948                 {
1949                   fNChTrCorrMCQuark->Fill(nTracksPercBckgSubst,pTbs); 
1950                   fSCMMCPerp->Fill(secondCentralR,jetPt);
1951                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
1952                     fNChTrMCPerp->Fill(nTracksAboveThresholdPerp,jetPt);
1953                   if(!IsEqualRel(fCurrentJetMinPtNT90Recalc, 7000.))
1954                     fNChTrCorrMCPerp->Fill(nTrUpThrPerpBckSubs,pTbs);
1955                 }
1956               break;
1957             case 3:
1958               fNChTr[2]->Fill(nTracksPerc,jetPt);
1959               fProcessPDG[2]->Fill(jetPt,evtype);
1960               fHistPtParton[2]->Fill(jetPt);
1961               fSCM[2]->Fill(secondCentralR,jetPt);
1962               fMinTrackPtInNTXh[0]->Fill(fMinTrackPtInNTX,jetPt,1); // 0 for pp MC
1963               fMaxTrackPtInNTXh[0]->Fill(fMaxTrackPtInNTX,jetPt); // 0 for MC
1964               if(fIsPossibleToSubstBckg)
1965                 {
1966                   fNChTrCorrMCQuark->Fill(nTracksPercBckgSubst,pTbs); 
1967                   fSCMMCPerp->Fill(secondCentralR,jetPt);
1968                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
1969                     fNChTrMCPerp->Fill(nTracksAboveThresholdPerp,jetPt);
1970                   if(!IsEqualRel(fCurrentJetMinPtNT90Recalc, 7000.))
1971                     fNChTrCorrMCPerp->Fill(nTrUpThrPerpBckSubs,pTbs);
1972                 }
1973               break;
1974             case 4:
1975               fNChTr[3]->Fill(nTracksPerc,jetPt);
1976               fProcessPDG[3]->Fill(jetPt,evtype);
1977               fHistPtParton[3]->Fill(jetPt);
1978               fSCM[3]->Fill(secondCentralR,jetPt);
1979               fMinTrackPtInNTXh[0]->Fill(fMinTrackPtInNTX,jetPt,1); // 0 for pp MC
1980               fMaxTrackPtInNTXh[0]->Fill(fMaxTrackPtInNTX,jetPt); // 0 for MC
1981               if(fIsPossibleToSubstBckg)
1982                 {
1983                   fNChTrCorrMCQuark->Fill(nTracksPercBckgSubst,pTbs); 
1984                   fSCMMCPerp->Fill(secondCentralR,jetPt);
1985                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
1986                     fNChTrMCPerp->Fill(nTracksAboveThresholdPerp,jetPt);
1987                   if(!IsEqualRel(fCurrentJetMinPtNT90Recalc, 7000.))
1988                     fNChTrCorrMCPerp->Fill(nTrUpThrPerpBckSubs,pTbs);
1989                 }
1990               break;
1991             case 5:
1992               fNChTr[4]->Fill(nTracksPerc,jetPt);
1993               fProcessPDG[4]->Fill(jetPt,evtype);
1994               fHistPtParton[4]->Fill(jetPt);
1995               fSCM[4]->Fill(secondCentralR,jetPt);
1996               fMinTrackPtInNTXh[0]->Fill(fMinTrackPtInNTX,jetPt,1); // 0 for pp MC
1997               fMaxTrackPtInNTXh[0]->Fill(fMaxTrackPtInNTX,jetPt); // 0 for MC
1998               if(fIsPossibleToSubstBckg)
1999                 {
2000                   fNChTrCorrMCQuark->Fill(nTracksPercBckgSubst,pTbs); 
2001                   fSCMMCPerp->Fill(secondCentralR,jetPt);
2002                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
2003                     fNChTrMCPerp->Fill(nTracksAboveThresholdPerp,jetPt);
2004                   if(!IsEqualRel(fCurrentJetMinPtNT90Recalc, 7000.))
2005                     fNChTrCorrMCPerp->Fill(nTrUpThrPerpBckSubs,pTbs);
2006                 }
2007               break;
2008             case 21:
2009               fNChTr[5]->Fill(nTracksPerc,jetPt);
2010               fProcessPDG[5]->Fill(jetPt,evtype);
2011               fHistPtParton[5]->Fill(jetPt);
2012               fSCM[5]->Fill(secondCentralR,jetPt);
2013               fMinTrackPtInNTXh[0]->Fill(fMinTrackPtInNTX,jetPt,1); // 0 for pp MC
2014               fMaxTrackPtInNTXh[0]->Fill(fMaxTrackPtInNTX,jetPt); // 0 for MC
2015               if(fIsPossibleToSubstBckg)
2016                 {
2017                   fNChTrCorrMCGluon->Fill(nTracksPercBckgSubst,pTbs); 
2018                   fSCMMCPerp->Fill(secondCentralR,jetPt);
2019                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
2020                     fNChTrMCPerp->Fill(nTracksAboveThresholdPerp,jetPt);
2021                   if(!IsEqualRel(fCurrentJetMinPtNT90Recalc, 7000.))
2022                     fNChTrCorrMCPerp->Fill(nTrUpThrPerpBckSubs,pTbs);
2023                 }
2024               break;          
2025             default:
2026               break;
2027             }
2028           AliDebug(4,Form("Sabor del jet numero:%d es: %d y se necesitaron %d tracks \n",indxmc,flavor,nTracksPerc)); 
2029           imcj++;
2030         } // MC jets for cycle
2031       nGenJets=imcj;
2032       for(Int_t u=0 ; u<mcJetCounter  ;u++)
2033         {
2034           if(u<7)
2035             fJetsMultPtMC->Fill(jetPts[u],mcJetCounter);
2036         }
2037       // if(fEnablePrints)
2038       //  {
2039       //   if(mcJetCounter>=3)
2040       //    printf("%i Jets inside acceptance at event number:%i \n",mcJetCounter,fEvtCount-1);
2041       //  }
2042       fNAccJetsMC->Fill(mcJetCounter,mcJetCounter);
2043     } // end if MC info in AOD
2044   
2045   if(!fUseOnlyMC) 
2046     {  
2047       // Primero que todo, debe de ir la seleccion de eventos reconstruidos:
2048       // 1. Que tenga un vertice reconstruido dentro de 10 cm.
2049       // Vertex info for reconstructed events
2050       AliAODVertex *pvx = fAOD->GetPrimaryVertex();
2051       if(!pvx)
2052         {
2053           AliError("No primary vertex!");
2054           return;
2055         }
2056       if(TMath::Abs(pvx->GetZ())>10.) // if the event vertex is larger than 10 cm, reject
2057         return;
2058       fZVertex->Fill(pvx->GetZ(),pvx->GetZ()); // vertex, provide number of accepted events as entries for reco jets
2059
2060       ///////////////////////////////////////
2061       // SECONDARY RECO BRANCH STUFF       // 
2062       // Get the secondary branch with the reconstructed jets 
2063       if(fBranchSecRec!="")
2064         {
2065           AliDebug(4,Form("fBranchSecRec was not default \n")); 
2066           TClonesArray *aodSecRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchSecRec.Data()));
2067           if(!aodSecRecJets)
2068             {
2069               AliError(Form("%s:%d no reconstructed Secondary Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchSecRec.Data())); 
2070               return;  //stop the analysis
2071             }
2072           AliDebug(4,Form("There are %d reconstructed jets from the secondary branch in this event \n", aodSecRecJets->GetEntries())); 
2073           Int_t recojetsSEC =  aodSecRecJets->GetEntries();
2074           fNJetsRDSeco->Fill(recojetsSEC,recojetsSEC);  // number of jets in the secondary branch
2075
2076           HasOverlapedCones(aodSecRecJets); // Procedure for tagging usable jets
2077                                             // Up 16 jets are flagged
2078           
2079           AliDebug(4,"Antes de realizar el loop jets reconstruidos del segundo branch \n"); 
2080           Int_t secondjetacccounter = 0;
2081           for (Int_t IDXS = 0; IDXS < recojetsSEC; IDXS++) 
2082             {
2083               AliDebug(4,Form("Number of current jet:%i \n",IDXS));
2084               AliAODJet *rjetsec = dynamic_cast<AliAODJet*>(aodSecRecJets->At(IDXS));
2085               if (!rjetsec) 
2086                 {
2087                   AliDebug(2,Form("ERROR: Could not receive jet %d from the second branch\n", IDXS)); 
2088                   continue;
2089                 }
2090               
2091               ///////////////////////////////////////////////////////////////////////////////
2092               ///// Part for Chritians plot of inclusive and leading jets comp at 2.76 TeV //
2093               if(!IsInsideAcceptance(rjetsec))  // old condition
2094                 continue;
2095               if(IDXS==0) // leading jet
2096                 fSecRecJetPtLeading->Fill(rjetsec->Pt());
2097               fSecRecJetPtInclusive->Fill(rjetsec->Pt()); // all
2098               ///// End of Christians Plot reco 2nd branch
2099               ///////////////////////////////////////////////////////////////////////////////
2100
2101               if(IDXS>15)
2102                 continue;
2103
2104               if(!fJetFlags[IDXS]) // If the jet is flaged kFALSE, not usable
2105                 continue;
2106
2107               fJetPtSec->Fill(rjetsec->Pt());
2108               fJetEtaSec->Fill(rjetsec->Eta(),rjetsec->Eta());
2109               fJetPhiSec->Fill(rjetsec->Phi(),rjetsec->Phi());
2110               secondjetacccounter++;
2111             }
2112           fNAccJetsRDSeco->Fill(secondjetacccounter,secondjetacccounter);
2113         }
2114       // END OF SECONDARY BRANCH STUFF     //
2115       ///////////////////////////////////////
2116
2117       // Get the branch with the reconstructed jets
2118       TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
2119       if(!aodRecJets)
2120         {
2121           AliError(Form("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data())); 
2122           return;
2123         }
2124       
2125       AliDebug(4,Form("There are %d reconstructed jets in this event\n", aodRecJets->GetEntries())); 
2126       Int_t recojets =  aodRecJets->GetEntries();
2127       fNJetsRD->Fill(recojets,recojets); // numero de jets directamente del branch
2128
2129       HasOverlapedCones(aodRecJets); // Procedure for tagging usable jets
2130                                     // Up 16 jets are flagged
2131
2132       AliDebug(4,"Antes de realizar el loop sobre AOD tracks \n"); 
2133       // Loop over AOD tracks
2134       Int_t tracksAOD = fAOD->GetNumberOfTracks(); 
2135       AliDebug(4,Form("Numero de tracks en el AOD:%d \n",tracksAOD));
2136       Double_t aodtracketa = 0.0;
2137       perpendicularPt = 0.0;
2138       px=0.0;
2139       py=0.0;
2140       pz=0.0;
2141       en=0.0;
2142       pTbs=0.0;
2143       etabs=0.0;
2144       phibs=0.0;
2145       fBckgSbsJet[0]=0.0;
2146       fBckgSbsJet[1]=0.0;
2147       fBckgSbsJet[2]=0.0;
2148       Int_t refNJMult = 0;
2149       Int_t myTotalMultRef = 0; 
2150       Int_t myTotalSoftMultRef = 0; 
2151       for(Int_t aodT = 0; aodT < tracksAOD; aodT++ )
2152         {
2153           AliAODTrack *aodtrack = fAOD->GetTrack(aodT);
2154           if(!aodtrack) continue;
2155           aodtracketa = TMath::Abs(aodtrack->Eta());
2156           if(aodtracketa>0.9) continue;
2157           if(!aodtrack->TestFilterBit(fFilterBit)) continue; //track filter selection
2158           fEtaAOD->Fill(aodtrack->Eta(),aodtrack->Eta());
2159           fPhiAOD->Fill(aodtrack->Phi(),aodtrack->Phi());
2160           fPtAOD->Fill(aodtrack->Pt(),aodtrack->Pt());
2161           if(fJetEvent) // if has an accepted jet, calculate the perpendicular cone
2162             {
2163               if(HasPerpendicularCone()) // If there is a perpendicular cone available
2164                 {
2165                   if(aodtrack->Pt()>fMinpTVal)
2166                     {
2167                       if(GetDeltaR(fEtaPerpCoord,fPhiPerpCoord,aodtrack->Eta(),aodtrack->Phi())<fJetRadius)
2168                         perpendicularPt = perpendicularPt + aodtrack->Pt();
2169                     }
2170                 }
2171             } // end if jet event
2172           //Total TPC multiplicity of primaries
2173           if(aodtrack->Pt()>fMinPtInGlobMult) // Min pT used in multiplicity estimation
2174             {
2175               myTotalMultRef++;
2176               if(aodtrack->Pt()<fMinpTValUE) continue; // pT cut  fMinpTValUE
2177               if(aodtrack->Pt()>fMaxpTValUE) continue; // pT cut  fMaxpTValUE
2178               myTotalSoftMultRef++;  
2179               if(!IsTrackInsideExcludedArea(aodtrack->Eta(), aodtrack->Phi(), aodRecJets)) 
2180                 refNJMult++;
2181             }
2182         } // end track loop over the event...
2183
2184       fRefMultWOJet->Fill(refMultiplicity,refNJMult); 
2185       fMultWOJetVZero->Fill(refNJMult,multFullV0);
2186       Double_t v0CorrMult = multFullV0 - GetV0ExcludedMultiplicity(aodRecJets);
2187       fRefMultFullV0->Fill(refMultiplicity,multFullV0);
2188       fRefMultV0Corr->Fill(refMultiplicity,v0CorrMult);
2189       fFullV0V0Corr->Fill(multFullV0,v0CorrMult);
2190       fRefAODTrackCount->Fill(refMultiplicity,myTotalMultRef);
2191       fTrackCountWOJet->Fill(myTotalMultRef,refNJMult);
2192
2193       if(recojets==1) // correlation for only monojet events
2194         {
2195           fRefMultFullV0UJ->Fill(refMultiplicity,multFullV0);
2196           fRefMultV0CorrUJ->Fill(refMultiplicity,v0CorrMult);
2197           fFullV0V0CorrUJ->Fill(multFullV0,v0CorrMult);
2198           fMultWOJetVZeroUJ->Fill(refNJMult,multFullV0);
2199           fRefMultWOJetUJ->Fill(refMultiplicity,refNJMult);
2200           fRefAODTrackCountUJ->Fill(refMultiplicity,myTotalMultRef);
2201           fTrackCountWOJetUJ->Fill(myTotalMultRef,refNJMult); 
2202         }
2203
2204       if(fJetEvent) // if has an accepted jet, calculate the perpendicular cone
2205         {
2206           if(HasPerpendicularCone()) // If there is a perpendicular cone available
2207             {
2208               px = perpendicularPt*TMath::Cos(fPhiPerpCoord);
2209               py = perpendicularPt*TMath::Sin(fPhiPerpCoord);
2210               pz = perpendicularPt/TMath::Tan(2.0*TMath::ATan(TMath::Exp(-fEtaPerpCoord)));
2211               en = TMath::Sqrt(px*px + py*py + pz*pz);
2212               fPerpCone->SetPxPyPzE(px, py, pz, en);
2213             }
2214         }
2215             
2216       // Loop to fill a pT spectrum of the reco jets
2217       Int_t irecj=0; // index for reconstructed jets to correlate
2218       Int_t nrectracks[6]={0};
2219       Double_t ptrecjet[6]={0};
2220       Double_t scmr[6]={0};
2221       Double_t aodtrackxi=0;
2222       Int_t ntxreco;
2223       Int_t nTRecAboveThresholdPerp=0; 
2224       Int_t ntxrecoRecalc;
2225       Int_t nTRecAboveThresholdPerpRecalc=0; 
2226             
2227       for(Int_t i=0; i<6; i++) // Reset per event
2228         {
2229           fHistContainerR4[i]->Reset();
2230           fHistContainerR3[i]->Reset();
2231           fHistContainerR2[i]->Reset();
2232         }
2233       
2234       Double_t jetPtsR[7]={0};  // to store the pt of the jets
2235       Int_t rJetCounter=0;    // counter of accepted reco jets 
2236       fIsPossibleToSubstBckg = kTRUE; // Initialize before the loop
2237       if(fJetEvent) // si tiene jets validos
2238         {
2239           if(!HasPerpendicularCone()) // pero no encontro un cono perpendicular libre
2240             fIsPossibleToSubstBckg = kFALSE; // if not perpendicular cone, set to kFALSE, so no perpendicular calculations available
2241         }
2242       
2243       AliDebug(4,"Antes de realizar el loop jets reconstruidos \n"); 
2244       for (Int_t indxrec = 0; indxrec < recojets; indxrec++) 
2245         {
2246           AliDebug(4,Form("Number of current jet:%i \n",indxrec));
2247           ntxreco = 0;
2248           ntxrecoRecalc = 0;
2249           fMinTrackPtInNTX=200.0;  //Initialize for each jet, overflown
2250           fMaxTrackPtInNTX=200.0;  //Initialize for each jet, overflown
2251           fMinTrackPtInNTXR=200.0;  //Initialize for each jet, overflown
2252           fMaxTrackPtInNTXR=200.0;  //Initialize for each jet, overflown
2253           deltaPhiPt = 0.0;
2254           deltaEtaPt = 0.0;
2255           deltaPhiSqPt = 0.0;
2256           deltaEtaSqPt = 0.0;
2257           totalTrackPt = 0.0;
2258           firstMomDeltPhi = 0.0;
2259           firstMomDeltEta = 0.0;
2260           secondMomDeltPhi = 0.0;
2261           secondMomDeltEta = 0.0;
2262           secondCentralPhi = 0.0;
2263           secondCentralEta = 0.0;
2264           secondCentralR = 0.0;
2265           fTotTracksInCone=0; // Underflown at initialization per jet
2266
2267           AliAODJet *rjet = dynamic_cast<AliAODJet*>(aodRecJets->At(indxrec));
2268           if (!rjet) 
2269             {
2270               AliDebug(2,Form("ERROR: Could not receive jet %d\n", indxrec)); 
2271               continue;
2272             }
2273           fJetEtaAll->Fill(rjet->Eta());// all jets
2274
2275           ///////////////////////////////////////////////////////////////////////////////
2276           ///// Part for Chritians plot of inclusive and leading jets comp at 2.76 TeV //
2277           if(!IsInsideAcceptance(rjet))  // old condition
2278             continue;
2279           if(indxrec==0) // leading jet
2280             fRecJetPtLeading->Fill(rjet->Pt());
2281           fRecJetPtInclusive->Fill(rjet->Pt()); // all
2282           fJetEtaOnlyTPCcut->Fill(rjet->Eta());// only eta acceptance cut for TPC
2283           ///// End of Christians Plot reco
2284           ///////////////////////////////////////////////////////////////////////////////
2285
2286           if(indxrec>15)
2287             continue;
2288
2289           if(!fJetFlags[indxrec]) // If the jet is flaged kFALSE, not usable
2290             continue;
2291           
2292           AliDebug(4,Form("Jet #%i is in the acceptance \n",indxrec));
2293           if(rJetCounter<7)
2294             jetPtsR[rJetCounter]=rjet->Pt();
2295           rJetCounter++;
2296           fJetPt->Fill(rjet->Pt());
2297           fJetEta->Fill(rjet->Eta(),rjet->Eta());
2298           fJetPhi->Fill(rjet->Phi(),rjet->Phi());
2299
2300           if(rjet->Pt()>10.)
2301             fJetEtaJetPt[0]->Fill(rjet->Eta());
2302           if(rjet->Pt()>30.)
2303             fJetEtaJetPt[1]->Fill(rjet->Eta());
2304           if(rjet->Pt()>50.)
2305             fJetEtaJetPt[2]->Fill(rjet->Eta());
2306           
2307           // Reco RefTracks check
2308           Bool_t rTrkFlagRec = kFALSE;
2309           Int_t trkinrecjet = rjet->GetRefTracks()->GetEntriesFast();
2310           if(trkinrecjet!=0&&!fForceNotTR)
2311             rTrkFlagRec = kTRUE;
2312           AliDebug(4,Form("Number of tracks in RefTracks reco jet:%i \n",trkinrecjet));
2313           if(rTrkFlagRec)
2314             {
2315               // Check the properties of the tracks in this jet with track refs
2316               AliDebug(4,Form("Checking composition in Reco jets with track refs")); 
2317               for(Int_t aodT = 0; aodT <trkinrecjet; aodT++ )
2318                 {
2319                   aodtrackxi=0;
2320                   AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(rjet->GetRefTracks()->At(aodT));
2321                   if(!aodtrack)
2322                     {
2323                       AliError("Error, no AOD Track!");
2324                       continue;
2325                     }
2326                   if(!aodtrack->TestFilterBit(fFilterBit))
2327                     {
2328                       //                      printf("Rejecting track from track refs due to wrong filterbit! \n");
2329                       continue; //track filter selection
2330                     }
2331                   deltaPhiPt += DeltaPhiTrack(rjet, aodtrack)*aodtrack->Pt();
2332                   deltaEtaPt += DeltaEtaTrack(rjet, aodtrack)*aodtrack->Pt();
2333                   deltaPhiSqPt += DeltaPhiSqTrack(rjet, aodtrack)*aodtrack->Pt();
2334                   deltaEtaSqPt += DeltaEtaSqTrack(rjet, aodtrack)*aodtrack->Pt();
2335                   totalTrackPt += aodtrack->Pt();
2336                   fTotTracksInCone++; // Counting tracks                  
2337                   if(!IsEqualRel(aodtrack->Pt(), 0.0)) //!IsEqualRel(totalTrackPtPerp, 0.0) //aodtrack->Pt()!=0
2338                     aodtrackxi= log(rjet->Pt()/aodtrack->Pt());
2339                   if(irecj<maxJetNum)
2340                     {
2341                       fHistContainerR4[irecj]->Fill(aodtrackxi,rjet->Pt());
2342                       if(!IsTrackInsideThisJet(aodtrack, rjet, 0.3)) continue;
2343                       fHistContainerR3[irecj]->Fill(aodtrackxi,rjet->Pt());
2344                       if(!IsTrackInsideThisJet(aodtrack, rjet, 0.2)) continue;
2345                       fHistContainerR2[irecj]->Fill(aodtrackxi,rjet->Pt());
2346                     }       
2347                 } //end loop over track references
2348               if(!IsEqualRel(totalTrackPt, 0.0)) //!IsEqualRel(totalTrackPt, 0.0)  // totalTrackPt!=0.0
2349                 {
2350                   firstMomDeltPhi = deltaPhiPt/totalTrackPt;
2351                   firstMomDeltEta = deltaEtaPt/totalTrackPt;
2352                   secondMomDeltPhi = deltaPhiSqPt/totalTrackPt;
2353                   secondMomDeltEta = deltaEtaSqPt/totalTrackPt;
2354                   secondCentralPhi = secondMomDeltPhi - firstMomDeltPhi*firstMomDeltPhi;
2355                   secondCentralEta = secondMomDeltEta - firstMomDeltEta*firstMomDeltEta;
2356                   secondCentralR = secondCentralPhi + secondCentralEta;
2357                 } // end if totalTrackPt!=0.0
2358               if(IsEqualRel(totalTrackPt, 0.0)) //!IsEqualRel(totalTrackPt, 0.0) // totalTrackPt==0.0
2359                 secondCentralR = 10.0; //overflow value
2360             } // end if there are track references
2361           
2362           if(!rTrkFlagRec)
2363             {
2364               // Check properties of the tracks in this jet without track refs
2365               AliDebug(4,Form("Checking composition in Reco jets without track refs")); 
2366               for(Int_t aodT = 0; aodT < tracksAOD; aodT++ )
2367                 {
2368                   AliAODTrack *aodtrack = fAOD->GetTrack(aodT);
2369                   if(!aodtrack) continue;
2370                   if(!IsTrackInsideThisJet(aodtrack, rjet, jfr)) continue;
2371                   if(!aodtrack->TestFilterBit(fFilterBit)) continue; //track filter selection
2372                   if(aodtrack->Pt()<fMinpTVal) continue; //DATA: PT CUT
2373                   deltaPhiPt += DeltaPhiTrack(rjet, aodtrack)*aodtrack->Pt();
2374                   deltaEtaPt += DeltaEtaTrack(rjet, aodtrack)*aodtrack->Pt();
2375                   deltaPhiSqPt += DeltaPhiSqTrack(rjet, aodtrack)*aodtrack->Pt();
2376                   deltaEtaSqPt += DeltaEtaSqTrack(rjet, aodtrack)*aodtrack->Pt();
2377                   totalTrackPt += aodtrack->Pt();
2378                   fTotTracksInCone++; // Counting tracks          
2379                   if(!IsEqualRel(aodtrack->Pt(), 0.0)) //!IsEqualRel(totalTrackPt, 0.0) //aodtrack->Pt()!=0
2380                     aodtrackxi= log(rjet->Pt()/aodtrack->Pt());
2381                   if(irecj<maxJetNum)
2382                     {
2383                       fHistContainerR4[irecj]->Fill(aodtrackxi,rjet->Pt());
2384                       if(!IsTrackInsideThisJet(aodtrack, rjet, 0.3)) continue;
2385                       fHistContainerR3[irecj]->Fill(aodtrackxi,rjet->Pt());
2386                       if(!IsTrackInsideThisJet(aodtrack, rjet, 0.2)) continue;
2387                       fHistContainerR2[irecj]->Fill(aodtrackxi,rjet->Pt());
2388                     }
2389                 } // end loop over tracks
2390               if(!IsEqualRel(totalTrackPt, 0.0)) //!IsEqualRel(totalTrackPt, 0.0) //totalTrackPt!=0.0
2391                 {
2392                   firstMomDeltPhi = deltaPhiPt/totalTrackPt;
2393                   firstMomDeltEta = deltaEtaPt/totalTrackPt;
2394                   secondMomDeltPhi = deltaPhiSqPt/totalTrackPt;
2395                   secondMomDeltEta = deltaEtaSqPt/totalTrackPt;
2396                   secondCentralPhi = secondMomDeltPhi - firstMomDeltPhi*firstMomDeltPhi;
2397                   secondCentralEta = secondMomDeltEta - firstMomDeltEta*firstMomDeltEta;
2398                   secondCentralR = secondCentralPhi + secondCentralEta;
2399                 } // end if totalTrackPt!=0.0 
2400               if(IsEqualRel(totalTrackPt, 0.0)) //!IsEqualRel(totalTrackPt, 0.0) // totalTrackPt==0.0
2401                 secondCentralR = 10.0; //overflow value
2402             } // end of no track references
2403           //Esto es lo anterior, toma al jet como es, y calcula NT90      
2404           ntxreco=GetNumberOfChargedTracks(ntx,rjet, tracksAOD, fAOD, jfr); // this call fixes the minimum pT track
2405           //Y aqui calcula cuantos tracks se necesitan arriba del threshold establecido en la linea anterior
2406           //esto debe ser para cada jet. Lo unico que se calcula una sola vez es el cono perpendicular  
2407           if(fIsPossibleToSubstBckg&&!IsEqualRel(fCurrentJetMinPtNT90, 7000.)) //and only if the method worked
2408             nTRecAboveThresholdPerp = GetNRecChargedTracksAboveThreshold(fPerpCone,tracksAOD, fAOD,jfr); //here one changes NTX
2409
2410           // correct the jet pT
2411           if(fIsPossibleToSubstBckg) // If there is a perpendicular cone available, substract backg and fill the new jet pT
2412             {
2413               pTbs= rjet->Pt()-fPerpCone->Pt();
2414               etabs= rjet->Eta();
2415               phibs= rjet->Phi();
2416               fBckgSbsJet[0]=pTbs; //pT
2417               fBckgSbsJet[1]=etabs; //eta
2418               fBckgSbsJet[2]=phibs; //phi
2419               // Now re-calculate nt90 for the energy corrected jet
2420               ntxrecoRecalc = GetRecalcNTXRec(ntx,rjet, tracksAOD, fAOD, jfr); //This call saves the new min pT
2421               // Now re-calculate the perpendicular cone NT90 background
2422               if(!IsEqualRel(fCurrentJetMinPtNT90Recalc, 7000.)) //calculating for the current jet, if the previous method worked
2423                 nTRecAboveThresholdPerpRecalc = GetRecalcNRecChTrUpThr(fPerpCone,tracksAOD, fAOD,jfr);
2424             }
2425
2426           // SCM perpendicular cone
2427           if(fIsPossibleToSubstBckg)
2428             {
2429               // To make sure, re-initialize
2430               deltaPhiPtPerp = 0.0;
2431               deltaEtaPtPerp = 0.0;
2432               deltaPhiSqPtPerp = 0.0;
2433               deltaEtaSqPtPerp = 0.0;
2434               totalTrackPtPerp = 0.0;
2435               firstMomDeltPhiPerp = 0.0;
2436               firstMomDeltEtaPerp = 0.0;
2437               secondMomDeltPhiPerp = 0.0;
2438               secondMomDeltEtaPerp = 0.0;
2439               secondCentralPhiPerp = 0.0;
2440               secondCentralEtaPerp = 0.0;
2441               secondCentralRPerp = 0.0;
2442               AliDebug(4,Form("Checking SCM of perpendicular cone in Reco jets")); 
2443               for(Int_t aodTperp = 0; aodTperp < tracksAOD; aodTperp++ )
2444                 { //fPerpCone
2445                   AliAODTrack *aodtrackperprec = fAOD->GetTrack(aodTperp);
2446                   if(!aodtrackperprec) continue;
2447                   if(!IsTrackInsideThisJet(aodtrackperprec, fPerpCone, jfr)) continue;
2448                   if(!aodtrackperprec->TestFilterBit(fFilterBit)) continue; //track filter selection
2449                   if(aodtrackperprec->Pt()<fMinpTVal) continue; //DATA: PT CUT
2450                   deltaPhiPtPerp += DeltaPhiTrack(fPerpCone, aodtrackperprec)*aodtrackperprec->Pt();
2451                   deltaEtaPtPerp += DeltaEtaTrack(fPerpCone, aodtrackperprec)*aodtrackperprec->Pt();
2452                   deltaPhiSqPtPerp += DeltaPhiSqTrack(fPerpCone, aodtrackperprec)*aodtrackperprec->Pt();
2453                   deltaEtaSqPtPerp += DeltaEtaSqTrack(fPerpCone, aodtrackperprec)*aodtrackperprec->Pt();
2454                   totalTrackPtPerp += aodtrackperprec->Pt();
2455                 } // end loop over tracks
2456               if(!IsEqualRel(totalTrackPtPerp, 0.0)) //!IsEqualRel(totalTrackPt, 0.0) // totalTrackPtPerp!=0.0
2457                 {
2458                   firstMomDeltPhiPerp = deltaPhiPtPerp/totalTrackPtPerp;
2459                   firstMomDeltEtaPerp = deltaEtaPtPerp/totalTrackPtPerp;
2460                   secondMomDeltPhiPerp = deltaPhiSqPtPerp/totalTrackPtPerp;
2461                   secondMomDeltEtaPerp = deltaEtaSqPtPerp/totalTrackPtPerp;
2462                   secondCentralPhiPerp = secondMomDeltPhiPerp - firstMomDeltPhiPerp*firstMomDeltPhiPerp;
2463                   secondCentralEtaPerp = secondMomDeltEtaPerp - firstMomDeltEtaPerp*firstMomDeltEtaPerp;
2464                   secondCentralRPerp = secondCentralPhiPerp + secondCentralEtaPerp;
2465                 } // end if totalTrackPt!=0.0
2466               if(IsEqualRel(totalTrackPtPerp, 0.0)) //!IsEqualRel(totalTrackPt, 0.0) // totalTrackPtPerp==0.0
2467                 secondCentralRPerp = 10.0; //overflow
2468             } 
2469
2470           ///// end of adding the SCM for the perpendicular cone
2471
2472           if(irecj<maxJetNum)
2473             {
2474               recJets[irecj]= *rjet;
2475               nrectracks[irecj] = ntxreco;
2476               ptrecjet[irecj] = rjet->Pt();
2477               scmr[irecj] = secondCentralR;
2478               AliDebug(4,Form("Para el jet reco num: %d se necesitaron %d tracks \n",irecj,nrectracks[irecj])); 
2479             }
2480           AliDebug(4,"Before filling the histograms for this jet \n");
2481           fNChTrRD->Fill(ntxreco,rjet->Pt());
2482           fSCMRD->Fill(secondCentralR,rjet->Pt());
2483           fProfNChTrRD->Fill(rjet->Pt(),ntxreco);
2484
2485           fNTXV0MultPt->Fill(ntxreco,v0CorrMult,rjet->Pt());
2486           fNTXCBMultPt->Fill(ntxreco,refNJMult,rjet->Pt());
2487
2488           //refNJMult
2489           // reference multiplicity stuff in pp, also filled in PbPb, but does not matter.
2490           // set to: V0 corrected multiplicity: v0CorrMult
2491           if(v0CorrMult<25)
2492             {
2493               fNChTrRDMult[0]->Fill(ntxreco,rjet->Pt());
2494               fSCMRDMult[0]->Fill(secondCentralR,rjet->Pt());
2495               fTotalJetCharge[0]->Fill(fCurrentJetCharge);
2496               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2497                 {
2498                   fNChTrRDMultOJ[0]->Fill(ntxreco,rjet->Pt());
2499                   fSCMRDMultOJ[0]->Fill(secondCentralR,rjet->Pt());
2500                   fJEtaRDMultOJ[0]->Fill(rjet->Eta());
2501                   fJetPtRDMultOJ[0]->Fill(rjet->Pt());
2502                 }
2503             }
2504           if(v0CorrMult>=25&&v0CorrMult<50)
2505             {
2506               fNChTrRDMult[1]->Fill(ntxreco,rjet->Pt());
2507               fSCMRDMult[1]->Fill(secondCentralR,rjet->Pt());
2508               fTotalJetCharge[1]->Fill(fCurrentJetCharge);
2509               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2510                 {
2511                   fNChTrRDMultOJ[1]->Fill(ntxreco,rjet->Pt());
2512                   fSCMRDMultOJ[1]->Fill(secondCentralR,rjet->Pt());
2513                   fJEtaRDMultOJ[1]->Fill(rjet->Eta());
2514                   fJetPtRDMultOJ[1]->Fill(rjet->Pt());
2515                 }
2516             }
2517           if(v0CorrMult>=50&&v0CorrMult<90)
2518             {
2519               fNChTrRDMult[2]->Fill(ntxreco,rjet->Pt());
2520               fSCMRDMult[2]->Fill(secondCentralR,rjet->Pt());
2521               fTotalJetCharge[2]->Fill(fCurrentJetCharge);
2522               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2523                 {
2524                   fNChTrRDMultOJ[2]->Fill(ntxreco,rjet->Pt());
2525                   fSCMRDMultOJ[2]->Fill(secondCentralR,rjet->Pt());
2526                   fJEtaRDMultOJ[2]->Fill(rjet->Eta());
2527                   fJetPtRDMultOJ[2]->Fill(rjet->Pt());
2528                 }
2529             }
2530           if(v0CorrMult>=90&&v0CorrMult<120)
2531             {
2532               fNChTrRDMult[3]->Fill(ntxreco,rjet->Pt());
2533               fSCMRDMult[3]->Fill(secondCentralR,rjet->Pt());
2534               fTotalJetCharge[3]->Fill(fCurrentJetCharge);
2535               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2536                 {
2537                   fNChTrRDMultOJ[3]->Fill(ntxreco,rjet->Pt());
2538                   fSCMRDMultOJ[3]->Fill(secondCentralR,rjet->Pt());
2539                   fJEtaRDMultOJ[3]->Fill(rjet->Eta());
2540                   fJetPtRDMultOJ[3]->Fill(rjet->Pt());
2541                 }
2542             }
2543           if(v0CorrMult>=120&&v0CorrMult<150)
2544             {
2545               fNChTrRDMult[4]->Fill(ntxreco,rjet->Pt());
2546               fSCMRDMult[4]->Fill(secondCentralR,rjet->Pt());
2547               fTotalJetCharge[4]->Fill(fCurrentJetCharge);
2548               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2549                 {
2550                   fNChTrRDMultOJ[4]->Fill(ntxreco,rjet->Pt());
2551                   fSCMRDMultOJ[4]->Fill(secondCentralR,rjet->Pt());
2552                   fJEtaRDMultOJ[4]->Fill(rjet->Eta());
2553                   fJetPtRDMultOJ[4]->Fill(rjet->Pt());
2554                 }
2555             }
2556           if(v0CorrMult>=150&&v0CorrMult<200)
2557             {
2558               fNChTrRDMult[5]->Fill(ntxreco,rjet->Pt());
2559               fSCMRDMult[5]->Fill(secondCentralR,rjet->Pt());
2560               fTotalJetCharge[5]->Fill(fCurrentJetCharge);
2561               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2562                 {
2563                   fNChTrRDMultOJ[5]->Fill(ntxreco,rjet->Pt());
2564                   fSCMRDMultOJ[5]->Fill(secondCentralR,rjet->Pt());
2565                   fJEtaRDMultOJ[5]->Fill(rjet->Eta());
2566                   fJetPtRDMultOJ[5]->Fill(rjet->Pt());
2567                 }
2568             }
2569           if(v0CorrMult>=200&&v0CorrMult<300)
2570             {
2571               fNChTrRDMult[6]->Fill(ntxreco,rjet->Pt());
2572               fSCMRDMult[6]->Fill(secondCentralR,rjet->Pt());
2573               fTotalJetCharge[6]->Fill(fCurrentJetCharge);
2574               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2575                 {
2576                   fNChTrRDMultOJ[6]->Fill(ntxreco,rjet->Pt());
2577                   fSCMRDMultOJ[6]->Fill(secondCentralR,rjet->Pt());
2578                   fJEtaRDMultOJ[6]->Fill(rjet->Eta());
2579                   fJetPtRDMultOJ[6]->Fill(rjet->Pt());
2580                 }
2581             }
2582           if(v0CorrMult>=300)
2583             {
2584               fNChTrRDMult[7]->Fill(ntxreco,rjet->Pt());
2585               fSCMRDMult[7]->Fill(secondCentralR,rjet->Pt());
2586               fTotalJetCharge[7]->Fill(fCurrentJetCharge);
2587               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2588                 {
2589                   fNChTrRDMultOJ[7]->Fill(ntxreco,rjet->Pt());
2590                   fSCMRDMultOJ[7]->Fill(secondCentralR,rjet->Pt());
2591                   fJEtaRDMultOJ[7]->Fill(rjet->Eta());
2592                   fJetPtRDMultOJ[7]->Fill(rjet->Pt());
2593                 }
2594             }
2595
2596           // 2nd. Reference: set to: TPC tracks minus jet, minus dijet area
2597           if(refNJMult<5) //&&refNJMult>1
2598             {
2599               fNChTrRDMultSE[0]->Fill(ntxreco,rjet->Pt());
2600               fSCMRDMultSE[0]->Fill(secondCentralR,rjet->Pt());
2601               fTotalJetChargeSE[0]->Fill(fCurrentJetCharge);
2602               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2603                 {
2604                   fNChTrRDMultSEOJ[0]->Fill(ntxreco,rjet->Pt());
2605                   fSCMRDMultSEOJ[0]->Fill(secondCentralR,rjet->Pt());
2606                   fJEtaRDMultSEOJ[0]->Fill(rjet->Eta());
2607                   fJetPtRDMultSEOJ[0]->Fill(rjet->Pt());
2608                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
2609                     fNChTrRecPerpMultSEOJ[0]->Fill(nTRecAboveThresholdPerp,rjet->Pt());  
2610                   if(fIsPossibleToSubstBckg) // if it was possible to calculate a perpendicular cone
2611                     {
2612                       fNChTrRecECorrPPMult->Fill(ntxrecoRecalc,pTbs,1); //filling mult bin
2613                       fNChTrRecPerpECorrPPMult->Fill(nTRecAboveThresholdPerpRecalc,pTbs,1); //filling mult bin
2614                     }
2615                 }
2616             }
2617           if(refNJMult>=5&&refNJMult<10)
2618             {
2619               fNChTrRDMultSE[1]->Fill(ntxreco,rjet->Pt());
2620               fSCMRDMultSE[1]->Fill(secondCentralR,rjet->Pt());
2621               fTotalJetChargeSE[1]->Fill(fCurrentJetCharge);
2622               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2623                 {
2624                   fNChTrRDMultSEOJ[1]->Fill(ntxreco,rjet->Pt());
2625                   fSCMRDMultSEOJ[1]->Fill(secondCentralR,rjet->Pt());
2626                   fJEtaRDMultSEOJ[1]->Fill(rjet->Eta());
2627                   fJetPtRDMultSEOJ[1]->Fill(rjet->Pt());
2628                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
2629                     fNChTrRecPerpMultSEOJ[1]->Fill(nTRecAboveThresholdPerp,rjet->Pt());  
2630                   if(fIsPossibleToSubstBckg) // if it was possible to calculate a perpendicular cone
2631                     {
2632                       fNChTrRecECorrPPMult->Fill(ntxrecoRecalc,pTbs,2); //filling mult bin
2633                       fNChTrRecPerpECorrPPMult->Fill(nTRecAboveThresholdPerpRecalc,pTbs,2); //filling mult bin
2634                     }
2635                 }
2636             }
2637           if(refNJMult>=10&&refNJMult<15)
2638             {
2639               fNChTrRDMultSE[2]->Fill(ntxreco,rjet->Pt());
2640               fSCMRDMultSE[2]->Fill(secondCentralR,rjet->Pt());
2641               fTotalJetChargeSE[2]->Fill(fCurrentJetCharge);
2642               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2643                 {
2644                   fNChTrRDMultSEOJ[2]->Fill(ntxreco,rjet->Pt());
2645                   fSCMRDMultSEOJ[2]->Fill(secondCentralR,rjet->Pt());
2646                   fJEtaRDMultSEOJ[2]->Fill(rjet->Eta());
2647                   fJetPtRDMultSEOJ[2]->Fill(rjet->Pt());
2648                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
2649                     fNChTrRecPerpMultSEOJ[2]->Fill(nTRecAboveThresholdPerp,rjet->Pt()); 
2650                   if(fIsPossibleToSubstBckg) // if it was possible to calculate a perpendicular cone
2651                     {
2652                       fNChTrRecECorrPPMult->Fill(ntxrecoRecalc,pTbs,3); //filling mult bin
2653                       fNChTrRecPerpECorrPPMult->Fill(nTRecAboveThresholdPerpRecalc,pTbs,3); //filling mult bin
2654                     } 
2655                 }
2656             }
2657           if(refNJMult>=15&&refNJMult<20)
2658             {
2659               fNChTrRDMultSE[3]->Fill(ntxreco,rjet->Pt());
2660               fSCMRDMultSE[3]->Fill(secondCentralR,rjet->Pt());
2661               fTotalJetChargeSE[3]->Fill(fCurrentJetCharge);
2662               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2663                 {
2664                   fNChTrRDMultSEOJ[3]->Fill(ntxreco,rjet->Pt());
2665                   fSCMRDMultSEOJ[3]->Fill(secondCentralR,rjet->Pt());
2666                   fJEtaRDMultSEOJ[3]->Fill(rjet->Eta());
2667                   fJetPtRDMultSEOJ[3]->Fill(rjet->Pt());
2668                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
2669                     fNChTrRecPerpMultSEOJ[3]->Fill(nTRecAboveThresholdPerp,rjet->Pt()); 
2670                   if(fIsPossibleToSubstBckg) // if it was possible to calculate a perpendicular cone
2671                     {
2672                       fNChTrRecECorrPPMult->Fill(ntxrecoRecalc,pTbs,4); //filling mult bin
2673                       fNChTrRecPerpECorrPPMult->Fill(nTRecAboveThresholdPerpRecalc,pTbs,4); //filling mult bin
2674                     } 
2675                 }
2676             }
2677           if(refNJMult>=20&&refNJMult<30)
2678             {
2679               fNChTrRDMultSE[4]->Fill(ntxreco,rjet->Pt());
2680               fSCMRDMultSE[4]->Fill(secondCentralR,rjet->Pt());
2681               fTotalJetChargeSE[4]->Fill(fCurrentJetCharge);
2682               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2683                 {
2684                   fNChTrRDMultSEOJ[4]->Fill(ntxreco,rjet->Pt());
2685                   fSCMRDMultSEOJ[4]->Fill(secondCentralR,rjet->Pt());
2686                   fJEtaRDMultSEOJ[4]->Fill(rjet->Eta());
2687                   fJetPtRDMultSEOJ[4]->Fill(rjet->Pt());
2688                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
2689                     fNChTrRecPerpMultSEOJ[4]->Fill(nTRecAboveThresholdPerp,rjet->Pt());  
2690                   if(fIsPossibleToSubstBckg) // if it was possible to calculate a perpendicular cone
2691                     {
2692                       fNChTrRecECorrPPMult->Fill(ntxrecoRecalc,pTbs,5); //filling mult bin
2693                       fNChTrRecPerpECorrPPMult->Fill(nTRecAboveThresholdPerpRecalc,pTbs,5); //filling mult bin
2694                     }
2695                 }
2696             }
2697           if(refNJMult>=30&&refNJMult<40)
2698             {
2699               fNChTrRDMultSE[5]->Fill(ntxreco,rjet->Pt());
2700               fSCMRDMultSE[5]->Fill(secondCentralR,rjet->Pt());
2701               fTotalJetChargeSE[5]->Fill(fCurrentJetCharge);
2702               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2703                 {
2704                   fNChTrRDMultSEOJ[5]->Fill(ntxreco,rjet->Pt());
2705                   fSCMRDMultSEOJ[5]->Fill(secondCentralR,rjet->Pt());
2706                   fJEtaRDMultSEOJ[5]->Fill(rjet->Eta());
2707                   fJetPtRDMultSEOJ[5]->Fill(rjet->Pt());
2708                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
2709                     fNChTrRecPerpMultSEOJ[5]->Fill(nTRecAboveThresholdPerp,rjet->Pt());  
2710                   if(fIsPossibleToSubstBckg) // if it was possible to calculate a perpendicular cone
2711                     {
2712                       fNChTrRecECorrPPMult->Fill(ntxrecoRecalc,pTbs,6); //filling mult bin
2713                       fNChTrRecPerpECorrPPMult->Fill(nTRecAboveThresholdPerpRecalc,pTbs,6); //filling mult bin
2714                     }
2715                 }
2716             }
2717           if(refNJMult>=40&&refNJMult<50)
2718             {
2719               fNChTrRDMultSE[6]->Fill(ntxreco,rjet->Pt());
2720               fSCMRDMultSE[6]->Fill(secondCentralR,rjet->Pt());
2721               fTotalJetChargeSE[6]->Fill(fCurrentJetCharge);
2722               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2723                 {
2724                   fNChTrRDMultSEOJ[6]->Fill(ntxreco,rjet->Pt());
2725                   fSCMRDMultSEOJ[6]->Fill(secondCentralR,rjet->Pt());
2726                   fJEtaRDMultSEOJ[6]->Fill(rjet->Eta());
2727                   fJetPtRDMultSEOJ[6]->Fill(rjet->Pt());
2728                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
2729                     fNChTrRecPerpMultSEOJ[6]->Fill(nTRecAboveThresholdPerp,rjet->Pt());  
2730                   if(fIsPossibleToSubstBckg) // if it was possible to calculate a perpendicular cone
2731                     {
2732                       fNChTrRecECorrPPMult->Fill(ntxrecoRecalc,pTbs,7); //filling mult bin
2733                       fNChTrRecPerpECorrPPMult->Fill(nTRecAboveThresholdPerpRecalc,pTbs,7); //filling mult bin
2734                     }
2735                 }
2736             }
2737           if(refNJMult>=50)
2738             {
2739               fNChTrRDMultSE[7]->Fill(ntxreco,rjet->Pt());
2740               fSCMRDMultSE[7]->Fill(secondCentralR,rjet->Pt());
2741               fTotalJetChargeSE[7]->Fill(fCurrentJetCharge);
2742               if(recojets==1) // if only one jet in the whole event, and inside acceptance
2743                 {
2744                   fNChTrRDMultSEOJ[7]->Fill(ntxreco,rjet->Pt());
2745                   fSCMRDMultSEOJ[7]->Fill(secondCentralR,rjet->Pt());
2746                   fJEtaRDMultSEOJ[7]->Fill(rjet->Eta());
2747                   fJetPtRDMultSEOJ[7]->Fill(rjet->Pt());
2748                   if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
2749                     fNChTrRecPerpMultSEOJ[7]->Fill(nTRecAboveThresholdPerp,rjet->Pt());  
2750                   if(fIsPossibleToSubstBckg) // if it was possible to calculate a perpendicular cone
2751                     {
2752                       fNChTrRecECorrPPMult->Fill(ntxrecoRecalc,pTbs,8); //filling mult bin
2753                       fNChTrRecPerpECorrPPMult->Fill(nTRecAboveThresholdPerpRecalc,pTbs,8); //filling mult bin
2754                     }
2755                 }
2756             }
2757
2758           if(fIsPossibleToSubstBckg) // if it was possible to calculate a perpendicular cone
2759             {
2760               if(!IsEqualRel(fCurrentJetMinPtNT90, 7000.))
2761                 fNChTrRecPerp->Fill(nTRecAboveThresholdPerp,rjet->Pt());  // these are my previous histos
2762               fSCMRecPerp->Fill(secondCentralRPerp,rjet->Pt());  // this are my previous histos
2763               if(!fIsHIevent) //if is a proton proton event
2764                 {
2765                   fNChTrRecECorr->Fill(ntxrecoRecalc,pTbs,1); //filling proton bin
2766                   if(!IsEqualRel(fCurrentJetMinPtNT90Recalc, 7000.))
2767                     {
2768                       fNChTrRecPerpECorr->Fill(nTRecAboveThresholdPerpRecalc,pTbs,1); //filling proton bin
2769                       if((rjet->Pt()>10.)&&(rjet->Pt()<20.))
2770                         {
2771                           fPtInPerpCon->Fill(fPerpCone->Pt(),1,1);
2772                           fTotTracksCone->Fill(fTotTracksInCone,1,1);
2773                           FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 1, 1);
2774                         }
2775                       if((rjet->Pt()>20.)&&(rjet->Pt()<30.))
2776                         {
2777                           fPtInPerpCon->Fill(fPerpCone->Pt(),2,1);
2778                           fTotTracksCone->Fill(fTotTracksInCone,2,1);
2779                           FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 1, 2);
2780                         }
2781                       if((rjet->Pt()>30.)&&(rjet->Pt()<40.))
2782                         {
2783                           fPtInPerpCon->Fill(fPerpCone->Pt(),3,1);
2784                           fTotTracksCone->Fill(fTotTracksInCone,3,1);
2785                           FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 1, 3);
2786                         }
2787                       if((rjet->Pt()>40.)&&(rjet->Pt()<50.))
2788                         {
2789                           fPtInPerpCon->Fill(fPerpCone->Pt(),4,1);
2790                           fTotTracksCone->Fill(fTotTracksInCone,4,1);
2791                           FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 1, 4);
2792                         }
2793                       if((rjet->Pt()>50.)&&(rjet->Pt()<60.))
2794                         {
2795                           fPtInPerpCon->Fill(fPerpCone->Pt(),5,1);
2796                           fTotTracksCone->Fill(fTotTracksInCone,5,1);
2797                           FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 1, 5);
2798                         }
2799                       if((rjet->Pt()>60.)&&(rjet->Pt()<70.))
2800                         {
2801                           fPtInPerpCon->Fill(fPerpCone->Pt(),6,1);
2802                           fTotTracksCone->Fill(fTotTracksInCone,6,1);
2803                           FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 1, 6);
2804                         }
2805                       if((rjet->Pt()>70.)&&(rjet->Pt()<80.))
2806                         {
2807                           fPtInPerpCon->Fill(fPerpCone->Pt(),7,1);
2808                           fTotTracksCone->Fill(fTotTracksInCone,7,1);
2809                           FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 1, 7);
2810                         }
2811                       if((rjet->Pt()>80.)&&(rjet->Pt()<90.))
2812                         {
2813                           fPtInPerpCon->Fill(fPerpCone->Pt(),8,1); 
2814                           fTotTracksCone->Fill(fTotTracksInCone,8,1);
2815                           FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 1, 8);
2816                         }
2817                     }
2818                 }
2819               if(fIsHIevent) //if is a PbPb event     
2820                 {
2821                   if(fEventCent>=0&&fEventCent<10.)
2822                     {
2823                       fNChTrRecECorr->Fill(ntxrecoRecalc,pTbs,2); //filling first centrality bin
2824                       fJetPtCentPbPbRaw->Fill(rjet->Pt(),2); 
2825                       fJetPtCentPbPbCorr->Fill(pTbs,2); 
2826                       if(!IsEqualRel(fCurrentJetMinPtNT90Recalc, 7000.))
2827                         {
2828                           fNChTrRecPerpECorr->Fill(nTRecAboveThresholdPerpRecalc,pTbs,2); //filling first centrality bin
2829                           fMinTrackPtInNTXRecalc->Fill(fMinTrackPtInNTXR,rjet->Pt(),2); 
2830                           fMinTrackPtInNTXh[1]->Fill(fMinTrackPtInNTX,rjet->Pt(),2); 
2831                           if((rjet->Pt()>10.)&&(rjet->Pt()<20.))
2832                             {
2833                               fPtInPerpCon->Fill(fPerpCone->Pt(),1,2);
2834                               fTotTracksCone->Fill(fTotTracksInCone,1,2);
2835                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 2, 1);
2836                             }
2837                           if((rjet->Pt()>20.)&&(rjet->Pt()<30.))
2838                             {
2839                               fPtInPerpCon->Fill(fPerpCone->Pt(),2,2);
2840                               fTotTracksCone->Fill(fTotTracksInCone,2,2);
2841                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 2, 2);
2842                             }
2843                           if((rjet->Pt()>30.)&&(rjet->Pt()<40.))
2844                             {
2845                               fPtInPerpCon->Fill(fPerpCone->Pt(),3,2);
2846                               fTotTracksCone->Fill(fTotTracksInCone,3,2);
2847                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 2, 3);
2848                             }
2849                           if((rjet->Pt()>40.)&&(rjet->Pt()<50.))
2850                             {
2851                               fPtInPerpCon->Fill(fPerpCone->Pt(),4,2);
2852                               fTotTracksCone->Fill(fTotTracksInCone,4,2);
2853                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 2, 4);
2854                             }
2855                           if((rjet->Pt()>50.)&&(rjet->Pt()<60.))
2856                             {
2857                               fPtInPerpCon->Fill(fPerpCone->Pt(),5,2);
2858                               fTotTracksCone->Fill(fTotTracksInCone,5,2);
2859                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 2, 5);
2860                             }
2861                           if((rjet->Pt()>60.)&&(rjet->Pt()<70.))
2862                             {
2863                               fPtInPerpCon->Fill(fPerpCone->Pt(),6,2);
2864                               fTotTracksCone->Fill(fTotTracksInCone,6,2);
2865                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 2, 6);
2866                             }
2867                           if((rjet->Pt()>70.)&&(rjet->Pt()<80.))
2868                             {
2869                               fPtInPerpCon->Fill(fPerpCone->Pt(),7,2);
2870                               fTotTracksCone->Fill(fTotTracksInCone,7,2);
2871                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 2, 7);
2872                             }
2873                           if((rjet->Pt()>80.)&&(rjet->Pt()<90.))
2874                             {
2875                               fPtInPerpCon->Fill(fPerpCone->Pt(),8,2);
2876                               fTotTracksCone->Fill(fTotTracksInCone,8,2);
2877                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 2, 8);
2878                             }
2879                         }
2880                     }
2881                   if(fEventCent>=10&&fEventCent<20.)
2882                     {
2883                       fNChTrRecECorr->Fill(ntxrecoRecalc,pTbs,3); //filling second centrality bin
2884                       fJetPtCentPbPbRaw->Fill(rjet->Pt(),3); 
2885                       fJetPtCentPbPbCorr->Fill(pTbs,3); 
2886                       if(!IsEqualRel(fCurrentJetMinPtNT90Recalc, 7000.))
2887                         {
2888                           fNChTrRecPerpECorr->Fill(nTRecAboveThresholdPerpRecalc,pTbs,3); //filling second centrality bin
2889                           fMinTrackPtInNTXRecalc->Fill(fMinTrackPtInNTXR,rjet->Pt(),3); 
2890                           fMinTrackPtInNTXh[1]->Fill(fMinTrackPtInNTX,rjet->Pt(),3); 
2891                           if((rjet->Pt()>10.)&&(rjet->Pt()<20.))
2892                             {
2893                               fPtInPerpCon->Fill(fPerpCone->Pt(),1,3);
2894                               fTotTracksCone->Fill(fTotTracksInCone,1,3);
2895                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 3, 1);
2896                             }
2897                           if((rjet->Pt()>20.)&&(rjet->Pt()<30.))
2898                             {
2899                               fPtInPerpCon->Fill(fPerpCone->Pt(),2,3);
2900                               fTotTracksCone->Fill(fTotTracksInCone,2,3);
2901                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 3, 2);
2902                             }
2903                           if((rjet->Pt()>30.)&&(rjet->Pt()<40.))
2904                             {
2905                               fPtInPerpCon->Fill(fPerpCone->Pt(),3,3);
2906                               fTotTracksCone->Fill(fTotTracksInCone,3,3);
2907                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 3, 3);
2908                             }
2909                           if((rjet->Pt()>40.)&&(rjet->Pt()<50.))
2910                             {
2911                               fPtInPerpCon->Fill(fPerpCone->Pt(),4,3);
2912                               fTotTracksCone->Fill(fTotTracksInCone,4,3);
2913                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 3, 4);
2914                             }
2915                           if((rjet->Pt()>50.)&&(rjet->Pt()<60.))
2916                             {
2917                               fPtInPerpCon->Fill(fPerpCone->Pt(),5,3);
2918                               fTotTracksCone->Fill(fTotTracksInCone,5,3);
2919                               FillPerpConeHisto(fPtDistInPerpConeRaw, tracksAOD, fAOD, 3, 5);