up from Megan
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalJetSpectraMECpA.cxx
1 // $Id: AliAnalysisTaskEmcalJetSpectraMECpA.cxx 3010 2012-06-10 05:40:56Z loizides $
2 //
3 // Jet spectrum task.
4 //
5 // Author: R.Reed, M.Connors
6
7 #include "AliAnalysisTaskEmcalJetSpectraMECpA.h"
8
9 #include <TCanvas.h>
10 #include <TChain.h>
11 #include <TClonesArray.h>
12 #include <TH1F.h>
13 #include <TH2F.h>
14 #include <TH3F.h>
15 #include <THnSparse.h>
16 #include <TList.h>
17 #include <TLorentzVector.h>
18 #include <TParameter.h>
19 #include <TParticle.h>
20 #include <TTree.h>
21 #include <TVector3.h>
22
23 #include "AliAODEvent.h"
24 #include "AliAnalysisManager.h"
25 #include "AliAnalysisTask.h"
26 #include "AliCentrality.h"
27 #include "AliESDEvent.h"
28 #include "AliESDInputHandler.h"
29 #include "AliEmcalJet.h"
30 #include "AliVCluster.h"
31 #include "AliRhoParameter.h"
32 #include "AliEmcalParticle.h"
33
34 ClassImp(AliAnalysisTaskEmcalJetSpectraMECpA)
35
36 //________________________________________________________________________
37 AliAnalysisTaskEmcalJetSpectraMECpA::AliAnalysisTaskEmcalJetSpectraMECpA() : 
38   AliAnalysisTaskEmcalJet("spectra",kFALSE), 
39   fHistRhovsCent(0),
40   fHistNjetvsCent(0)
41 {
42   // Default constructor.
43   for (Int_t i = 0;i<6;++i){
44     fHistJetPtvsTrackPt[i]      = 0;
45     fHistRawJetPtvsTrackPt[i]   = 0;
46     fHistTrackPt[i]             = 0;
47     fHistEP0[i]                 = 0;
48     fHistEP0A[i]                = 0;
49     fHistEP0C[i]                = 0;
50     fHistEPAvsC[i]              = 0;
51     fHistJetPtvsdEP[i]          = 0;
52     fHistJetPtvsdEPBias[i]      = 0;
53     fHistRhovsEP[i]             = 0;
54     fHistJetPtEtaPhi[i]         = 0;
55
56   }
57   SetCentralityEstimator("V0A");
58   SetMakeGeneralHistograms(kTRUE);
59 }
60
61 //________________________________________________________________________
62 AliAnalysisTaskEmcalJetSpectraMECpA::AliAnalysisTaskEmcalJetSpectraMECpA(const char *name) :
63   AliAnalysisTaskEmcalJet(name,kTRUE),
64   fHistRhovsCent(0),
65   fHistNjetvsCent(0)
66  { 
67    for (Int_t i = 0;i<6;++i){
68     fHistJetPtvsTrackPt[i]      = 0;
69     fHistRawJetPtvsTrackPt[i]   = 0;
70     fHistTrackPt[i]             = 0;
71     fHistEP0[i]                 = 0;
72     fHistEP0A[i]                = 0;
73     fHistEP0C[i]                = 0;
74     fHistEPAvsC[i]              = 0;
75     fHistJetPtvsdEP[i]          = 0;
76     fHistJetPtvsdEPBias[i]      = 0;
77     fHistRhovsEP[i]             = 0;
78     fHistJetPtEtaPhi[i]         = 0;
79    }
80    SetCentralityEstimator("V0A");
81    SetMakeGeneralHistograms(kTRUE);
82  }
83
84 //________________________________________________________________________
85 void AliAnalysisTaskEmcalJetSpectraMECpA::UserCreateOutputObjects()
86 {
87   if (! fCreateHisto)
88     return;
89   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
90
91   fHistRhovsCent             = new TH2F("RhovsCent",              "RhovsCent",             100, 0.0, 100.0, 500, 0, 500);
92   fHistNjetvsCent            = new TH2F("NjetvsCent",             "NjetvsCent",            100, 0.0, 100.0, 100, 0, 100);
93
94   TString name;
95   TString title;
96   for (Int_t i = 0;i<6;++i){
97     name = TString(Form("JetPtvsTrackPt_%i",i));
98     title = TString(Form("Jet pT vs Leading Track pT cent bin %i",i));
99     fHistJetPtvsTrackPt[i] = new TH2F(name,title,1000,-500,500,100,0,100);
100     fOutput->Add(fHistJetPtvsTrackPt[i]);
101     name = TString(Form("RawJetPtvsTrackPt_%i",i));
102     title = TString(Form("Raw Jet pT vs Leading Track pT cent bin %i",i));
103     fHistRawJetPtvsTrackPt[i] = new TH2F(name,title,1000,-500,500,100,0,100);
104     fOutput->Add(fHistRawJetPtvsTrackPt[i]);
105     name = TString(Form("TrackPt_%i",i));
106     title = TString(Form("Track pT cent bin %i",i));
107     fHistTrackPt[i] = new TH1F(name,title,1000,0,200);
108     fOutput->Add(fHistTrackPt[i]);
109    
110     name = TString(Form("EP0_%i",i));
111     title = TString(Form("EP VZero cent bin %i",i));
112     fHistEP0[i] = new TH1F(name,title,100,-TMath::Pi(),TMath::Pi());
113     fOutput->Add(fHistEP0[i]);
114     name = TString(Form("EP0A_%i",i));
115     title = TString(Form("EP VZero cent bin %i",i));
116     fHistEP0A[i] = new TH1F(name,title,100,-TMath::Pi(),TMath::Pi());
117     fOutput->Add(fHistEP0A[i]);
118     name = TString(Form("EP0C_%i",i));
119     title = TString(Form("EP VZero cent bin %i",i));
120     fHistEP0C[i] = new TH1F(name,title,100,-TMath::Pi(),TMath::Pi());
121     fOutput->Add(fHistEP0C[i]);
122     name = TString(Form("EPAvsC_%i",i));
123     title = TString(Form("EP VZero cent bin %i",i));
124     fHistEPAvsC[i] = new TH2F(name,title,100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());
125     fOutput->Add(fHistEPAvsC[i]);
126     name = TString(Form("JetPtvsdEP_%i",i));
127     title = TString(Form("Jet pt vs dEP cent bin %i",i));
128     fHistJetPtvsdEP[i] = new TH2F(name,title,1000,-500,500,400,-2*TMath::Pi(),2*TMath::Pi());
129     fOutput->Add(fHistJetPtvsdEP[i]);
130     name = TString(Form("JetPtvsdEPBias_%i",i));
131     title = TString(Form("Bias Jet pt vs dEP cent bin %i",i));
132     fHistJetPtvsdEPBias[i] = new TH2F(name,title,1000,-500,500,400,-2*TMath::Pi(),2*TMath::Pi());
133     fOutput->Add(fHistJetPtvsdEPBias[i]);
134     name = TString(Form("JetPtvsEP_%i",i));
135     title = TString(Form("Jet pt vs EP cent bin %i",i));
136     fHistJetPtvsEP[i] = new TH2F(name,title,1000,-500,500,400,-2*TMath::Pi(),2*TMath::Pi());
137     fOutput->Add(fHistJetPtvsEP[i]);
138     name = TString(Form("JetPtvsEPBias_%i",i));
139     title = TString(Form("Bias Jet pt vs EP cent bin %i",i));
140     fHistJetPtvsEPBias[i] = new TH2F(name,title,1000,-500,500,400,-2*TMath::Pi(),2*TMath::Pi());
141     fOutput->Add(fHistJetPtvsEPBias[i]);
142     name = TString(Form("RhovsEP_%i",i));
143     title = TString(Form("Rho vs EP cent bin %i",i));
144     fHistRhovsEP[i] = new TH2F(name,title,500,0,500,400,-2*TMath::Pi(),2*TMath::Pi());
145     fOutput->Add(fHistRhovsEP[i]);
146
147     name = TString(Form("JetPtEtaPhi_%i",i));
148     title = TString(Form("JetPtEtaPhi_%i",i));
149     fHistJetPtEtaPhi[i] = new TH3F(name,title,300,-100,200,26,-0.8,0.8,400,-2*TMath::Pi(),2*TMath::Pi());
150     fOutput->Add(fHistJetPtEtaPhi[i]);
151
152
153   }
154
155   
156   fOutput->Add(fHistRhovsCent);
157   fOutput->Add(fHistNjetvsCent);
158   
159    PostData(1, fOutput);
160 }
161
162 //________________________________________________________________________
163
164 Int_t AliAnalysisTaskEmcalJetSpectraMECpA::GetCentBin(Double_t cent) const 
165 {
166   // Get centrality bin.
167
168   Int_t centbin = -1;
169   if (cent>=0 && cent<10)
170     centbin = 0;
171   else if (cent>=10 && cent<20)
172     centbin = 1;
173   else if (cent>=20 && cent<30)
174     centbin = 2;
175   else if (cent>=30 && cent<40)
176     centbin = 3;
177   else if (cent>=40 && cent<50)
178     centbin = 4;
179   else if (cent>=50 && cent<90)
180     centbin = 5;
181   return centbin;
182 }
183
184 //________________________________________________________________________
185
186 Float_t AliAnalysisTaskEmcalJetSpectraMECpA:: RelativePhi(Double_t mphi,Double_t vphi) const
187 {
188   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
189   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
190   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
191   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
192   double dphi = mphi-vphi;
193   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
194   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
195
196   return dphi;//dphi in [-Pi, Pi]                                                                                                    
197 }
198
199
200 //________________________________________________________________________
201 Bool_t AliAnalysisTaskEmcalJetSpectraMECpA::Run()
202 {
203
204
205   Int_t centbin = GetCentBin(fCent);
206   //for pp analyses we will just use the first centrality bin
207   if (centbin == -1)
208     centbin = 0;
209
210   if (!fTracks){
211     cout << "Tracks not found: " << fCent << endl;
212     return kTRUE;
213   }
214   
215   const Int_t nTrack = fTracks->GetEntriesFast();
216   for (int i = 0;i<nTrack;i++){
217     AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(i));
218     if (!track)
219       continue;
220     fHistTrackPt[centbin]->Fill(track->Pt());
221   }
222
223   fHistEP0[centbin]->Fill(fEPV0);
224   fHistEP0A[centbin]->Fill(fEPV0A);
225   fHistEP0C[centbin]->Fill(fEPV0C);
226   fHistEPAvsC[centbin]->Fill(fEPV0A,fEPV0C);
227   TString fRhoScaledName = fRhoName;
228   fRho = GetRhoFromEvent(fRhoScaledName);
229   fRhoVal = fRho->GetVal();
230   fHistRhovsCent->Fill(fCent,fRhoVal);
231   fHistRhovsEP[centbin]->Fill(fRhoVal,fEPV0);
232   const Int_t Njets = fJets->GetEntriesFast();
233
234   Int_t NjetAcc = 0;
235   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
236      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
237      if (!jet)
238        continue; 
239      if (jet->Area()==0)
240        continue;
241      if (jet->Pt()<0.1)
242        continue;
243      if (jet->MaxTrackPt()>100)
244        continue;
245      if (! AcceptJet(jet))
246        continue;
247      //  jets.push_back(jet);
248      NjetAcc++;
249      Double_t jetPt = -500;
250      jetPt = jet->Pt()-jet->Area()*fRhoVal;    
251      fHistJetPtvsTrackPt[centbin]->Fill(jetPt,jet->MaxTrackPt());
252      fHistRawJetPtvsTrackPt[centbin]->Fill(jet->Pt(),jet->MaxTrackPt());
253
254      fHistJetPtEtaPhi[centbin]->Fill(jet->Pt(),jet->Eta(),jet->Phi());
255
256      fHistJetPtvsdEP[centbin]->Fill(jetPt,RelativePhi((fEPV0+TMath::Pi()),jet->Phi()));
257      fHistJetPtvsEP[centbin]->Fill(jetPt,fEPV0);
258      if (jet->MaxTrackPt()>5.0){
259        fHistJetPtvsdEPBias[centbin]->Fill(jetPt,RelativePhi((fEPV0+TMath::Pi()),jet->Phi()));
260        fHistJetPtvsEPBias[centbin]->Fill(jetPt,fEPV0);
261      }
262   }
263   
264   fHistNjetvsCent->Fill(fCent,NjetAcc);
265   return kTRUE;
266 }      
267
268
269
270
271