AliAODEvent::GetHeader() returns AliVHeader
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / FourierDecomposition / AliDhcTask.cxx
1 // Dihadron correlations task - simple task to read ESD or AOD input,
2 // calculate same- and mixed-event correlations, and fill THnSparse
3 // output. -A. Adare, Apr 2011
4
5 #include <TAxis.h>
6 #include <TCanvas.h>
7 #include <TChain.h>
8 #include <TFormula.h>
9 #include <TH1.h>
10 #include <TH2.h>
11 #include <TH3.h>
12 #include <THn.h>
13 #include <TProfile2D.h>
14 #include <TROOT.h>
15 #include <TTree.h>
16 #include "AliAnalysisUtils.h"
17 #include "AliAODEvent.h"
18 #include "AliAODInputHandler.h"
19 #include "AliAnalysisManager.h"
20 #include "AliAnalysisTask.h"
21 #include "AliCentrality.h"
22 #include "AliDhcTask.h"
23 #include "AliESDEvent.h"
24 #include "AliESDInputHandler.h"
25 #include "AliESDMuonTrack.h"
26 #include "AliESDtrackCuts.h"
27 #include "AliPool.h"
28 #include "AliVParticle.h"
29
30 using std::cout;
31 using std::endl;
32
33 ClassImp(AliDhcTask)
34
35 //________________________________________________________________________
36 AliDhcTask::AliDhcTask()
37 : AliAnalysisTaskSE(), fVerbosity(0), fEtaMax(1), fZVtxMax(10), fPtMin(0.25), fPtMax(15),
38   fTrackDepth(1000), fPoolSize(200), fTracksName(), fDoWeights(kFALSE), fFillMuons(kFALSE),
39   fPtTACrit(kTRUE), fAllTAHists(kFALSE), fMixInEtaT(kFALSE),
40   fEtaTLo(-1.0), fEtaTHi(1.0), fEtaALo(-1.0), fEtaAHi(1.0), fOmitFirstEv(kTRUE),
41   fDoFillSame(kFALSE), fDoMassCut(kFALSE), fCheckVertex(kTRUE), fClassName(),
42   fCentMethod("V0M"), fNBdeta(20), fNBdphi(36), fTriggerMatch(kTRUE),
43   fBPtT(0x0), fBPtA(0x0), fBCent(0x0), fBZvtx(0x0),
44   fMixBCent(0x0), fMixBZvtx(0x0), fHEffT(0x0), fHEffA(0x0),
45   fESD(0x0), fAOD(0x0), fOutputList(0x0), fHEvt(0x0), fHTrk(0x0), fHPoolReady(0x0),
46   fHPtAss(0x0), fHPtTrg(0x0), fHPtTrgEvt(0x0),
47   fHPtTrgNorm1S(0x0), fHPtTrgNorm1M(0x0), fHPtTrgNorm2S(0x0), fHPtTrgNorm2M(0x0),
48   fHCent(0x0), fHZvtx(0x0), fNbins(0), fHSs(0x0), fHMs(0x0), fHPts(0x0), fHSMass(0x0), fHMMass(0x0),
49   fHQATp(0x0), fHQAAp(0x0), fHQATpCorr(0x0), fHQAApCorr(0x0),
50   fHQATm(0x0), fHQAAm(0x0), fHQATmCorr(0x0), fHQAAmCorr(0x0),
51   fHPtCentT(0x0), fHPtCentA(0x0), fIndex(0x0),
52   fCentrality(99), fZVertex(99), fEsdTPCOnly(0), fPoolMgr(0),
53   fUtils(0x0)
54 {
55   // Constructor
56 }
57
58 //________________________________________________________________________
59 AliDhcTask::AliDhcTask(const char *name, Bool_t def) 
60 : AliAnalysisTaskSE(name), fVerbosity(0), fEtaMax(1), fZVtxMax(10), fPtMin(0.25), fPtMax(15),
61   fTrackDepth(1000), fPoolSize(200), fTracksName(), fDoWeights(kFALSE), fFillMuons(kFALSE),
62   fPtTACrit(kTRUE), fAllTAHists(kFALSE), fMixInEtaT(kFALSE),
63   fEtaTLo(-1.0), fEtaTHi(1.0), fEtaALo(-1.0), fEtaAHi(1.0), fOmitFirstEv(kTRUE),
64   fDoFillSame(kFALSE), fDoMassCut(kFALSE), fCheckVertex(kTRUE), fClassName(),
65   fCentMethod("V0M"), fNBdeta(20), fNBdphi(36), fTriggerMatch(kTRUE),
66   fBPtT(0x0), fBPtA(0x0), fBCent(0x0), fBZvtx(0x0),
67   fMixBCent(0x0), fMixBZvtx(0x0), fHEffT(0x0), fHEffA(0x0),
68   fESD(0x0), fAOD(0x0), fOutputList(0x0), fHEvt(0x0), fHTrk(0x0), fHPoolReady(0x0),
69   fHPtAss(0x0), fHPtTrg(0x0), fHPtTrgEvt(0x0),
70   fHPtTrgNorm1S(0x0), fHPtTrgNorm1M(0x0), fHPtTrgNorm2S(0x0), fHPtTrgNorm2M(0x0),
71   fHCent(0x0), fHZvtx(0x0), fNbins(0), fHSs(0x0), fHMs(0x0), fHPts(0x0), fHSMass(0x0), fHMMass(0x0),
72   fHQATp(0x0), fHQAAp(0x0), fHQATpCorr(0x0), fHQAApCorr(0x0),
73   fHQATm(0x0), fHQAAm(0x0), fHQATmCorr(0x0), fHQAAmCorr(0x0),
74   fHPtCentT(0x0), fHPtCentA(0x0), fIndex(0x0),
75   fCentrality(99), fZVertex(99), fEsdTPCOnly(0), fPoolMgr(0),
76   fUtils(0x0)
77 {
78   // Constructor
79
80   // Define input and output slots here
81   // Input slot #0 works with a TChain
82   DefineInput(0, TChain::Class());
83   // Output slot #0 id reserved by the base class for AOD
84   // Output slot #1 writes into a TH1 container
85   DefineOutput(1, TList::Class());
86
87   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,Tracks "
88                "AOD:header,tracks,vertices,";
89
90   if (def) {
91     Double_t ptt[4] = {0.25, 1.0, 2.0, 15.0};
92     fBPtT  = new TAxis(3,ptt); 
93     Double_t pta[4] = {0.25, 1.0, 2.0, 15.0};
94     fBPtA  = new TAxis(3,pta); 
95     Double_t cent[2] = {-100.0, 100.0};
96     fBCent = new TAxis(1,cent);
97     Double_t zvtx[2] = {-10, 10};
98     fBZvtx = new TAxis(1,zvtx);
99     Double_t centmix[2] = {-100.0, 100.0};
100     fMixBCent = new TAxis(1,centmix);
101     Double_t zvtxmix[9] = {-10,-6,-4,-2,0,2,4,6,10};
102     fMixBZvtx = new TAxis(8,zvtxmix);
103   }
104 }
105
106 //________________________________________________________________________
107 void AliDhcTask::UserCreateOutputObjects()
108 {
109   // Create histograms
110   // Called once (per slave on PROOF!)
111   PrintDhcSettings();
112
113   fOutputList = new TList();
114   fOutputList->SetOwner(1);
115
116   fUtils = new AliAnalysisUtils();
117
118   BookHistos();
119   InitEventMixer(); 
120   PostData(1, fOutputList);
121 }
122
123 //________________________________________________________________________
124 void AliDhcTask::PrintDhcSettings()
125 {
126   AliInfo(Form("Dhc Task %s settings",fName.Data()));
127   AliInfo(Form(" centrality estimator %s", fCentMethod.Data()));
128   AliInfo(Form(" using tracks named %s", fTracksName.Data()));
129   AliInfo(Form(" efficiency correct triggers? %d", fHEffT!=0));
130   if (fHEffT!=0) {
131     AliInfo(Form(" %d dimensions (t)", fHEffT->GetNdimensions()));
132   }
133   AliInfo(Form(" efficiency correct associates? %d", fHEffA!=0));
134   if (fHEffA!=0) {
135     AliInfo(Form(" %d dimensions (a)", fHEffA->GetNdimensions()));
136   }
137   AliInfo(Form(" fill muons? %d", fFillMuons));
138   AliInfo(Form(" use pTT > pTA criterion? %d", fPtTACrit));
139   AliInfo(Form(" create all pTT, pTA hists? %d", fAllTAHists));
140   AliInfo(Form(" Mix in eta_T bins instead of z_vertex? %d", fMixInEtaT));
141   AliInfo(Form(" trigger eta range %f .. %f", fEtaTLo, fEtaTHi));
142   AliInfo(Form(" associate eta range %f .. %f", fEtaALo, fEtaAHi));
143   AliInfo(Form(" fill same event in any case? %d", fDoFillSame));
144   AliInfo(Form(" do invariant mass cut? %d", fDoMassCut));
145   AliInfo(Form(" omit first event? %d\n", fOmitFirstEv));
146   AliInfo(Form(" check the vertex? %d\n", fCheckVertex));
147 }
148
149 //________________________________________________________________________
150 void AliDhcTask::BookHistos()
151 {
152   // Book histograms.
153
154   if (fVerbosity > 1) {
155     AliInfo(Form("Number of pt(t) bins: %d", fBPtT->GetNbins()));
156     for (Int_t i=1; i<=fBPtT->GetNbins(); i++) {
157       AliInfo(Form("pt(t) bin %d, %f to %f", i, fBPtT->GetBinLowEdge(i), fBPtT->GetBinUpEdge(i)));
158     }
159     AliInfo(Form("Number of pt(a) bins: %d", fBPtA->GetNbins()));
160     for (Int_t i=1; i<=fBPtA->GetNbins(); i++) {
161       AliInfo(Form("pt(a) bin %d, %f to %f", i, fBPtA->GetBinLowEdge(i), fBPtA->GetBinUpEdge(i)));
162     }
163   }
164
165   Int_t nPtAssc=fBPtA->GetNbins();
166   Int_t nPtTrig=fBPtT->GetNbins();
167   Int_t nCent=fBCent->GetNbins();
168   Int_t nZvtx=fBZvtx->GetNbins();
169   Double_t ptt[nPtTrig+1];
170   ptt[0] = fBPtT->GetBinLowEdge(1);
171   for (Int_t i=1; i<=nPtTrig; i++) {
172     ptt[i] = fBPtT->GetBinUpEdge(i);
173   }
174   Double_t pta[nPtAssc+1];
175   pta[0] = fBPtA->GetBinLowEdge(1);
176   for (Int_t i=1; i<=nPtAssc; i++) {
177     pta[i] = fBPtA->GetBinUpEdge(i);
178   }
179   Double_t cent[nCent+1];
180   cent[0] = fBCent->GetBinLowEdge(1);
181   for (Int_t i=1; i<=nCent; i++) {
182     cent[i] = fBCent->GetBinUpEdge(i);
183   }
184   Double_t zvtx[nZvtx+1];
185   zvtx[0] = fBZvtx->GetBinLowEdge(1);
186   for (Int_t i=1; i<=nZvtx; i++) {
187     zvtx[i] = fBZvtx->GetBinUpEdge(i);
188   }
189   
190   fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", nZvtx, zvtx, nCent, cent);
191   fOutputList->Add(fHEvt);
192   fHTrk = new TH2F("fHTrk", "Track-level variables",
193                    100, 0, TMath::TwoPi(), 100, -fEtaMax, +fEtaMax);
194   fOutputList->Add(fHTrk);
195   
196   // Centrality mixing axis
197   Int_t nCentBins=fMixBCent->GetNbins();
198   Double_t centBins[nCentBins+1];
199   centBins[0] = fMixBCent->GetBinLowEdge(1);
200   for (Int_t i=1; i<=nCentBins; i++) {
201     centBins[i] = fMixBCent->GetBinUpEdge(i);
202   }
203   // Z-vertex mixing axis
204   Int_t nZvtxBins=fMixBZvtx->GetNbins();
205   Double_t zvtxbin[nZvtxBins+1];
206   zvtxbin[0] = fMixBZvtx->GetBinLowEdge(1);
207   for (Int_t i=1; i<=nZvtxBins; i++) {
208     zvtxbin[i] = fMixBZvtx->GetBinUpEdge(i);
209   }
210   fHPoolReady = new TH2F("fHPoolReady","mixing started", nZvtxBins, zvtxbin, nCentBins, centBins);
211   fOutputList->Add(fHPoolReady);
212   
213   fHPtAss = new TH1F("fHPtAss","PtAssoc;P_{T} (GeV/c) [GeV/c]",nPtAssc,pta);
214   fOutputList->Add(fHPtAss);
215   fHPtTrg = new TH1F("fHPtTrg","PtTrig;P_{T} (GeV/c) [GeV/c]",nPtTrig,ptt);
216   fOutputList->Add(fHPtTrg);
217   fHPtTrgEvt = new TH1F("fHPtTrgEvt","PtTrig;P_{T} (GeV/c) [GeV/c]",nPtTrig,ptt);
218   fHPtTrgNorm1S = new TH3F("fHPtTrgNorm1S","PtTrig;P_{T} (GeV/c) [GeV/c];centrality;z_{vtx}",nPtTrig,ptt,nCent,cent,nZvtx,zvtx);
219   fOutputList->Add(fHPtTrgNorm1S);
220   fHPtTrgNorm1M = (TH3*) fHPtTrgNorm1S->Clone("fHPtTrgNorm1M");
221   fOutputList->Add(fHPtTrgNorm1M);
222   fHPtTrgNorm2S = (TH3*) fHPtTrgNorm1S->Clone("fHPtTrgNorm2S");
223   fOutputList->Add(fHPtTrgNorm2S);
224   fHPtTrgNorm2M = (TH3*) fHPtTrgNorm1S->Clone("fHPtTrgNorm2M");
225   fOutputList->Add(fHPtTrgNorm2M);
226   
227   fHCent = new TH1F("fHCent","Cent;bins",nCent,cent);
228   fOutputList->Add(fHCent);
229   fHZvtx = new TH1F("fHZvtx","Zvertex;bins",nZvtx,zvtx);
230   fOutputList->Add(fHZvtx);
231   
232   fHQATp = new TH3F("fHQATp","QA trigger;p_{T} (GeV/c);#eta;#phi",
233                     100,0.0,10.0,
234                     40,fEtaTLo,fEtaTHi,
235                     36,0.0,TMath::TwoPi());
236   fOutputList->Add(fHQATp);
237   fHQAAp = new TH3F("fHQAAp","QA associated;p_{T} (GeV/c);#eta;#phi",
238                     100,0.0,10.0,
239                     40,fEtaALo,fEtaAHi,
240                     36,0.0,TMath::TwoPi());
241   fOutputList->Add(fHQAAp);
242   fHQATpCorr = (TH3 *) fHQATp->Clone("fHQATpCorr");
243   fOutputList->Add(fHQATpCorr);
244   fHQAApCorr = (TH3 *) fHQAAp->Clone("fHQAApCorr");
245   fOutputList->Add(fHQAApCorr);
246   fHQATm = (TH3 *) fHQATp->Clone("fHQATm");
247   fOutputList->Add(fHQATm);
248   fHQAAm = (TH3 *) fHQAAp->Clone("fHQAAm");
249   fOutputList->Add(fHQAAm);
250   fHQATmCorr = (TH3 *) fHQATm->Clone("fHQATmCorr");
251   fOutputList->Add(fHQATmCorr);
252   fHQAAmCorr = (TH3 *) fHQAAm->Clone("fHQAAmCorr");
253   fOutputList->Add(fHQAAmCorr);
254
255   fHPtCentT = new TH2F("fHPtCentT",Form("trigger particles;p_{T} (GeV/c);centrality (%s)",fCentMethod.Data()),
256                        100,0.0,10.0,
257                        100,cent[0],cent[nCent]);
258   fOutputList->Add(fHPtCentT);
259   fHPtCentA = new TH2F("fHPtCentA",Form("associated particles;p_{T} (GeV/c);centrality (%s)",fCentMethod.Data()),
260                        100,0.0,10.0,
261                        100,cent[0],cent[nCent]);
262   fOutputList->Add(fHPtCentA);
263
264   fNbins  = nPtTrig*nPtAssc*nCent*nZvtx;
265   fHSs    = new TH2*[fNbins];
266   fHMs    = new TH2*[fNbins];
267   fHPts   = new TH1*[fNbins];
268   fHSMass = new TH1*[fNbins];
269   fHMMass = new TH1*[fNbins];
270
271   fIndex = new TFormula("GlobIndex","(t-1)*[0]*[1]*[2]+(z-1)*[0]*[1]+(x-1)*[0]+(y-1)+0*[4]");
272   fIndex->SetParameters(nPtTrig,nPtAssc,nZvtx,nCent);
273   fIndex->SetParNames("NTrigBins","NAssocBins", "NZvertexBins", "NCentBins");
274   fOutputList->Add(fIndex);
275   
276   Int_t count = 0;
277   for (Int_t c=1; c<=nCent; ++c) {
278     for (Int_t z=1; z<=nZvtx; ++z) {
279       for (Int_t t=1; t<=nPtTrig; ++t) {
280         for (Int_t a=1; a<=nPtAssc; ++a) {
281           fHSs[count]    = 0;
282           fHMs[count]    = 0;
283           fHPts[count]   = 0;
284           fHSMass[count] = 0;
285           fHMMass[count] = 0;
286           if ((a>t)&&!fAllTAHists) {
287             ++count;
288             continue;
289           }
290           if (z==1 && t==1 && a==1) {
291             TString title(Form("cen=%d (%.1f to %.1f)",
292                                c, fBCent->GetBinLowEdge(c), fBCent->GetBinUpEdge(c)));
293             fHSMass[count] = new TH1F(Form("hSMass%d",count), Form("Mass Same Event %s;m (GeV)",title.Data()), 10000, 0,10);
294             fOutputList->Add(fHSMass[count]);
295             fHMMass[count] = new TH1F(Form("hMMass%d",count), Form("Mass Mixed Event %s;m (GeV)",title.Data()), 10000, 0,10);
296             fOutputList->Add(fHMMass[count]);
297           }
298           if (t==1 && a==1) {
299             TString title(Form("cen=%d (%.1f to %.1f), zVtx=%d (%.1f to %.1f)",
300                                c, fBCent->GetBinLowEdge(c), fBCent->GetBinUpEdge(c),
301                                z, fBZvtx->GetBinLowEdge(z), fBZvtx->GetBinUpEdge(z)));
302             fHPts[count] = new TH1F(Form("hPt%d",count), Form("Ptdist %s;p_{T} (GeV/c)",title.Data()), 200,0,20);
303             fOutputList->Add(fHPts[count]);
304           }
305           TString title(Form("cen=%d (%.1f to %.1f), zVtx=%d (%.1f to %.1f), trig=%d (%.1f to %.1f), assoc=%d (%.1f to %.1f)",
306                              c, fBCent->GetBinLowEdge(c), fBCent->GetBinUpEdge(c),
307                              z, fBZvtx->GetBinLowEdge(z), fBZvtx->GetBinUpEdge(z),
308                              t, fBPtT->GetBinLowEdge(t),  fBPtT->GetBinUpEdge(t),
309                              a, fBPtA->GetBinLowEdge(a),  fBPtA->GetBinUpEdge(a)));
310           fHSs[count] = new TH2F(Form("hS%d",count), Form("Signal %s;#Delta #varphi;#Delta #eta",title.Data()),
311                                  fNBdphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),fNBdeta,-2*fEtaMax,2*fEtaMax);
312           fHMs[count] = new TH2F(Form("hM%d",count), Form("Mixed %s;#Delta #varphi;#Delta #eta",title.Data()),
313                                  fNBdphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),fNBdeta,-2*fEtaMax,2*fEtaMax);
314           fOutputList->Add(fHSs[count]);
315           fOutputList->Add(fHMs[count]);
316           Int_t tcount = (Int_t)fIndex->Eval(t,a,z,c);
317           if (fVerbosity>5)
318             cout << count << " " << tcount << ": " << title << endl;
319           if (count != tcount)
320             AliFatal(Form("Index mismatch: %d %d", count, tcount));
321           ++count;
322         }
323       }
324     }
325   }
326 }
327
328 //________________________________________________________________________
329 void AliDhcTask::SetAnaMode(Int_t iAna)
330 {
331   if (iAna==kHH) {
332     SetFillMuons(kFALSE);
333     fEtaTLo = -1.6;
334     fEtaTHi = +1.6;
335     fEtaALo = -1.6;
336     fEtaAHi = +1.6;
337   } else if (iAna==kMuH) {
338     SetFillMuons(kTRUE);
339     fEtaTLo = -5.0;
340     fEtaTHi = -2.0;
341     fEtaALo = -1.6;
342     fEtaAHi = +1.6;
343   } else if (iAna==kHMu) {
344     SetFillMuons(kTRUE);
345     fEtaTLo = -1.6;
346     fEtaTHi = +1.6;
347     fEtaALo = -5.0;
348     fEtaAHi = -2.0;
349   } else if (iAna==kMuMu) {
350     SetFillMuons(kTRUE);
351     fEtaTLo = -5.0;
352     fEtaTHi = -2.0;
353     fEtaALo = -5.0;
354     fEtaAHi = -2.0;
355   } else if (iAna==kPSide) {
356     SetFillMuons(kFALSE);
357     fEtaTLo = -1.6;
358     fEtaTHi = -0.465;
359     fEtaALo = -1.6;
360     fEtaAHi = +1.6;
361   } else if (iAna==kASide) {
362     SetFillMuons(kFALSE);
363     fEtaTLo = +0.465;
364     fEtaTHi = +1.6;
365     fEtaALo = -1.6;
366     fEtaAHi = +1.6;
367   } else {
368     AliInfo(Form("Unrecognized analysis option: %d", iAna));
369   }
370 }
371
372 //________________________________________________________________________
373 void AliDhcTask::InitEventMixer()
374 {
375   // The effective pool size in events is set by trackDepth, so more
376   // low-mult events are required to maintain the threshold than
377   // high-mult events. Centrality pools are indep. of data histogram
378   // binning, no need to match.
379
380   // Centrality pools
381   Int_t nCentBins=fMixBCent->GetNbins();
382   Double_t centBins[nCentBins+1];
383   centBins[0] = fMixBCent->GetBinLowEdge(1);
384   for (Int_t i=1; i<=nCentBins; i++) {
385     centBins[i] = fMixBCent->GetBinUpEdge(i);
386   }
387  
388   // Z-vertex pools
389   Int_t nZvtxBins=fMixBZvtx->GetNbins();
390   Double_t zvtxbin[nZvtxBins+1];
391   zvtxbin[0] = fMixBZvtx->GetBinLowEdge(1);
392   for (Int_t i=1; i<=nZvtxBins; i++) {
393     zvtxbin[i] = fMixBZvtx->GetBinUpEdge(i);
394   }
395
396   fPoolMgr = new AliEvtPoolManager();
397   fPoolMgr->SetTargetTrackDepth(fTrackDepth);
398   if (fVerbosity>4)
399     fPoolMgr->SetDebug(1);
400   fPoolMgr->InitEventPools(fPoolSize, nCentBins, centBins, nZvtxBins, zvtxbin);
401 }
402
403 //________________________________________________________________________
404 void AliDhcTask::UserExec(Option_t *) 
405 {
406   // Main loop, called for each event.
407
408   LoadBranches();
409
410   if (fOmitFirstEv) {
411     if (fUtils->IsFirstEventInChunk(InputEvent())) 
412       return;
413   }
414
415   // Get event pointers, check for signs of life
416   Int_t dType = -1;       // Will be set to kESD or kAOD.
417   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
418   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
419   if (fESD)
420     dType = kESD;
421   else if (fAOD)
422     dType = kAOD;
423   else {
424     AliError("Neither fESD nor fAOD available");
425     return;
426   }
427
428   // select specific trigger classes?
429   if (fClassName.Length()>0) {
430     TString strFiredClass;
431     if (fESD)
432       strFiredClass = fESD->GetFiredTriggerClasses();
433     else
434       strFiredClass = fAOD->GetFiredTriggerClasses();
435     
436     if (fVerbosity > 10) {
437       AliInfo(Form("Trigger class selection: This event has classes %s", strFiredClass.Data()));
438       AliInfo(Form("Trigger class selection: selecting classes %s", fClassName.Data()));
439     }
440
441     TObjArray *arrClass = fClassName.Tokenize(",");
442     Int_t nClass = arrClass->GetEntries();
443     
444     TString strOneClass;
445     Bool_t bAccEvent = kFALSE;
446     for (Int_t iClass=0; iClass<nClass; iClass++) {
447       strOneClass = arrClass->At(iClass)->GetName();
448       if (strFiredClass.Contains(strOneClass))
449         bAccEvent = kTRUE;
450     }
451     
452     if (!bAccEvent)
453       return;
454     
455     if (fVerbosity > 10) {
456       AliInfo(Form("After trigger class selection: This event has classes %s", strFiredClass.Data()));
457     }
458   }
459
460   Bool_t mcgen = 0;
461   if (fTracksName.Contains("Gen"))
462     mcgen = 1;
463
464   // Centrality, vertex, other event variables...
465   if (mcgen) {
466     fZVertex = 0;
467     TList *list = InputEvent()->GetList();
468     TClonesArray *tcaTracks = 0;
469     if (list) 
470       tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
471     if (tcaTracks) 
472       fCentrality = tcaTracks->GetEntries();
473   } else {
474     if (dType == kESD) {
475       const AliESDVertex* vertex = fESD->GetPrimaryVertex();
476       fZVertex = vertex->GetZ();
477       if (fCheckVertex && !VertexOk()) {
478         if (fVerbosity > 1)
479           AliInfo(Form("Event REJECTED (ESD vertex not OK). z = %.1f", fZVertex));
480         return;
481       }
482       if(fESD->GetCentrality()) {
483         fCentrality = 
484           fESD->GetCentrality()->GetCentralityPercentile(fCentMethod);
485       }
486     } else if (dType == kAOD) {
487       const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
488       fZVertex = vertex->GetZ();
489       if (fCheckVertex && !VertexOk()) {
490         if (fVerbosity > 1)
491           Info("Exec()", "Event REJECTED (AOD vertex not OK). z = %.1f", fZVertex);
492         return;
493       }
494       AliAODHeader * header = dynamic_cast<AliAODHeader*>(fAOD);
495       if(!header) AliFatal("Not a standard AOD");
496
497       const AliCentrality *aodCent = header->GetCentralityP();
498       if (aodCent) {
499         fCentrality = aodCent->GetCentralityPercentile(fCentMethod);
500       }
501     }
502   }
503
504   // Fill event histogram
505   fHEvt->Fill(fZVertex, fCentrality);
506   if (fCentrality > fBCent->GetXmax() || fCentrality < fBCent->GetXmin()) {
507     if (fVerbosity > 1)
508       AliInfo(Form("Event REJECTED (centrality out of range). fCentrality = %.1f", fCentrality));
509     return;
510   }
511
512   // Get pool containing tracks from other events like this one
513   AliEvtPool* pool = fPoolMgr->GetEventPool(fCentrality, fZVertex);
514   if (!pool) {
515     AliWarning(Form("No pool found. Centrality %f, ZVertex %f", 
516                     fCentrality, fZVertex));
517     return;
518   }
519
520   // Get array of selected tracks
521   MiniEvent* sTracks = new MiniEvent(0); // deleted by pool mgr.
522   if (dType == kESD) {
523     GetESDTracks(sTracks);
524   } else if (dType == kAOD) {
525     GetAODTracks(sTracks);
526   }
527
528   if (sTracks->size()==0) {
529     AliWarning(Form("Track array empty"));
530     delete sTracks;
531     return;
532   }
533
534   if (!pool->IsReady()) {
535     pool->UpdatePool(sTracks);
536     if (fDoFillSame) { // fill same event right away if requested
537       Correlate(*sTracks, *sTracks, kSameEvt);  
538     }
539     PostData(1, fOutputList);
540     return;
541   }
542
543   if (pool->IsFirstReady()) {
544     fHPoolReady->Fill(fZVertex,fCentrality);
545     // recover events missed before the pool is ready
546     Int_t nEvs = pool->GetCurrentNEvents();
547     if (nEvs>1) {
548       for (Int_t i=0; i<nEvs; ++i) {
549         MiniEvent* evI = pool->GetEvent(i);
550         for (Int_t j=0; j<nEvs; ++j) {
551           MiniEvent* evJ = pool->GetEvent(j);
552           if (i==j) {
553             if (!fDoFillSame) { // do not fill the same events twice
554               Correlate(*evI, *evJ, kSameEvt);
555             }
556           } else {
557             Correlate(*evI, *evJ, kDiffEvt);
558           }
559         }
560       }
561     }
562   } else { /* standard case: same event, then mix*/
563     Correlate(*sTracks, *sTracks, kSameEvt);  
564     Int_t nMix = pool->GetCurrentNEvents();
565     for (Int_t jMix=0; jMix<nMix; ++jMix) {
566       MiniEvent* bgTracks = pool->GetEvent(jMix);
567       Correlate(*sTracks, *bgTracks, kDiffEvt);
568     }
569   }
570
571   pool->UpdatePool(sTracks);
572   PostData(1, fOutputList);
573 }
574
575 //________________________________________________________________________
576 void AliDhcTask::GetESDTracks(MiniEvent* miniEvt)
577 {
578   // Loop twice: 1. Count sel. tracks. 2. Fill vector.
579
580   if (fTracksName.IsNull()) {
581     const AliESDVertex *vtxSPD = fESD->GetPrimaryVertexSPD();
582     if (!vtxSPD)
583       return;
584     
585     Int_t nTrax = fESD->GetNumberOfTracks();
586     if (fVerbosity > 2)
587       AliInfo(Form("%d tracks in event",nTrax));
588
589     // Loop 1.
590     Int_t nSelTrax = 0;
591     TObjArray arr(nTrax);
592     arr.SetOwner(1);
593
594     if (!fEsdTPCOnly) {
595       fEsdTPCOnly = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
596       //fEsdTPCOnly->SetMinNClustersTPC(70);
597       fEsdTPCOnly->SetMinNCrossedRowsTPC(70);
598       fEsdTPCOnly->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
599     }
600
601     for (Int_t i = 0; i < nTrax; ++i) {
602       AliESDtrack* esdtrack = fESD->GetTrack(i);
603       if (!esdtrack) {
604         AliError(Form("Couldn't get ESD track %d\n", i));
605         continue;
606       }
607       Bool_t trkOK = fEsdTPCOnly->AcceptTrack(esdtrack);
608       if (!trkOK)
609         continue;
610       Double_t pt = esdtrack->Pt();
611       Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
612       if (!ptOK)
613         continue;
614       Double_t eta = esdtrack->Eta();
615       if (TMath::Abs(eta) > fEtaMax)
616         continue;
617
618       // create a tpc only track
619       AliESDtrack *newtrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
620       if(!newtrack)
621         continue;
622       if (newtrack->Pt()<=0) {
623         delete newtrack;
624         continue;
625       }
626
627       AliExternalTrackParam exParam;
628       Bool_t relate = newtrack->RelateToVertexTPC(vtxSPD,fESD->GetMagneticField(),kVeryBig,&exParam);
629       if (!relate) {
630         delete newtrack;
631         continue;
632       }
633
634       // set the constraint parameters to the track
635       newtrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
636
637       pt = newtrack->Pt();
638       ptOK = pt >= fPtMin && pt < fPtMax;
639       if (!ptOK) {
640         delete newtrack;
641         continue;
642       }
643       eta  = esdtrack->Eta();
644       if (TMath::Abs(eta) > fEtaMax) {
645         delete newtrack;
646         continue;
647       }
648       arr.Add(newtrack);
649       nSelTrax++;
650     }
651
652     for(Int_t itrack = 0; itrack < nSelTrax; itrack++) {
653       AliVTrack *esdtrack = static_cast<AliESDtrack*>(arr.At(itrack));
654       if(!esdtrack) {
655         AliError(Form("ERROR: Could not retrieve esdtrack %d",itrack));
656         continue;
657       }
658       Double_t pt = esdtrack->Pt();
659       Double_t eta  = esdtrack->Eta();
660       Double_t phi  = esdtrack->Phi();
661       Int_t    sign = esdtrack->Charge() > 0 ? 1 : -1;
662       if (IsTrigger(eta,pt)||IsAssociated(eta,pt))
663         miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
664     }
665   } else {
666     TList *list = InputEvent()->GetList();
667     TClonesArray *tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
668
669     if (!tcaTracks){
670       AliError("Ptr to tcaTracks zero");
671       return;
672     }
673
674     const Int_t ntracks = tcaTracks->GetEntries();
675     Int_t nGoodTracks = 0;
676     // count good tracks
677     for (Int_t itrack = 0; itrack < ntracks; itrack++) {
678       AliVTrack *vtrack = static_cast<AliVTrack*>(tcaTracks->At(itrack));
679       if (!vtrack) {
680         AliError(Form("ERROR: Could not retrieve vtrack %d",itrack));
681         continue;
682       }
683       Double_t pt   = vtrack->Pt();
684       Double_t eta  = vtrack->Eta();
685       if (IsTrigger(eta,pt)||IsAssociated(eta,pt))
686         nGoodTracks++;
687     }
688     if (miniEvt)
689       miniEvt->reserve(nGoodTracks);
690     else {
691       AliError("Ptr to miniEvt zero");
692       return;
693     }
694     // fill good tracks into minievent
695     for (Int_t itrack = 0; itrack < ntracks; itrack++) {
696       AliVTrack *vtrack = static_cast<AliVTrack*>(tcaTracks->At(itrack));
697       if (!vtrack) {
698         AliError(Form("ERROR: Could not retrieve vtrack %d",itrack));
699         continue;
700       }
701       Double_t pt   = vtrack->Pt();
702       Double_t eta  = vtrack->Eta();
703       Double_t phi  = vtrack->Phi();
704       Int_t    sign = vtrack->Charge() > 0 ? 1 : -1;
705       if (IsTrigger(eta,pt)||IsAssociated(eta,pt))
706         miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
707     }
708   }
709   
710   if (fFillMuons) {
711     // count good muons
712     Int_t nGoodMuons = 0;
713     for (Int_t iMu = 0; iMu<fESD->GetNumberOfMuonTracks(); iMu++) {
714       AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iMu);
715       if (muonTrack) {
716         if (!IsGoodMUONtrack(*muonTrack))
717           continue;
718         Double_t ptMu   = muonTrack->Pt();
719         Double_t etaMu  = muonTrack->Eta();
720         if ((IsTrigger(etaMu,ptMu)||IsAssociated(etaMu,ptMu)) && (IsGoodMUONtrack(*muonTrack)))
721           nGoodMuons++;
722       }
723     }
724     miniEvt->reserve(miniEvt->size()+nGoodMuons);
725     // fill them into the mini event
726     for (Int_t iMu = 0; iMu<fESD->GetNumberOfMuonTracks(); iMu++) {
727       AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iMu);
728       if (muonTrack) {
729         if (!IsGoodMUONtrack(*muonTrack)) 
730           continue;
731         Double_t ptMu   = muonTrack->Pt();
732         Double_t etaMu  = muonTrack->Eta();
733         Double_t phiMu  = muonTrack->Phi();
734         Int_t    signMu = muonTrack->Charge() > 0 ? 1 : -1;
735         if (IsTrigger(etaMu,ptMu)||IsAssociated(etaMu,ptMu))
736           miniEvt->push_back(AliMiniTrack(ptMu, etaMu, phiMu, signMu));
737       }
738     }
739   }
740 }
741
742 //________________________________________________________________________
743 void AliDhcTask::GetAODTracks(MiniEvent* miniEvt)
744 {
745   // Loop twice: 1. Count sel. tracks. 2. Fill vector.
746
747   if (fTracksName.IsNull()) {
748     Int_t nTrax = fAOD->GetNumberOfTracks();
749     Int_t nSelTrax = 0;
750
751     if (fVerbosity > 2)
752       AliInfo(Form("%d tracks in event",nTrax));
753
754     // Loop 1.
755     for (Int_t i = 0; i < nTrax; ++i) {
756       AliAODTrack* aodtrack = fAOD->GetTrack(i);
757       if (!aodtrack) {
758         AliError(Form("Couldn't get AOD track %d\n", i));
759         continue;
760       }
761       // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
762       UInt_t tpcOnly = 1 << 7;
763       Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
764       if (!trkOK)
765         continue;
766       Double_t pt = aodtrack->Pt();
767       Double_t eta = aodtrack->Eta();
768       if (IsTrigger(eta,pt)||IsAssociated(eta,pt))
769         nSelTrax++;
770     }
771
772     if (miniEvt)
773       miniEvt->reserve(nSelTrax);
774     else {
775       AliError("!miniEvt");
776       return;
777     }
778   
779     // Loop 2.  
780     for (Int_t i = 0; i < nTrax; ++i) {
781       AliAODTrack* aodtrack = fAOD->GetTrack(i);
782       if (!aodtrack) {
783         AliError(Form("Couldn't get AOD track %d\n", i));
784         continue;
785       }
786     
787       // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
788       UInt_t tpcOnly = 1 << 7;
789       Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
790       if (!trkOK)
791         continue;
792       Double_t pt = aodtrack->Pt();
793       Double_t eta  = aodtrack->Eta();
794       Double_t phi  = aodtrack->Phi();
795       Int_t    sign = aodtrack->Charge() > 0 ? 1 : -1;
796       if (IsTrigger(eta,pt)||IsAssociated(eta,pt))
797         miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
798     }
799   } else {
800     TList *list = InputEvent()->GetList();
801     TClonesArray *tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
802
803     if (!tcaTracks){
804       AliError("Ptr to tcaTracks zero");
805       return;
806     }
807
808     const Int_t ntracks = tcaTracks->GetEntries();
809     Int_t nGoodTracks = 0;
810     // count good tracks
811     for (Int_t itrack = 0; itrack < ntracks; itrack++) {
812       AliVTrack *vtrack = static_cast<AliVTrack*>(tcaTracks->At(itrack));
813       if (!vtrack) {
814         AliError(Form("ERROR: Could not retrieve vtrack %d",itrack));
815         continue;
816       }
817       Double_t pt   = vtrack->Pt();
818       Double_t eta  = vtrack->Eta();
819       if (IsTrigger(eta,pt)||IsAssociated(eta,pt))
820         nGoodTracks++;
821     }
822     if (miniEvt)
823       miniEvt->reserve(nGoodTracks);
824     else {
825       AliError("Ptr to miniEvt zero");
826       return;
827     }
828     // fill good tracks into minievent
829     for (Int_t itrack = 0; itrack < ntracks; itrack++) {
830       AliVTrack *vtrack = static_cast<AliVTrack*>(tcaTracks->At(itrack));
831       if (!vtrack) {
832         AliError(Form("ERROR: Could not retrieve vtrack %d",itrack));
833         continue;
834       }
835       Double_t pt   = vtrack->Pt();
836       Double_t eta  = vtrack->Eta();
837       Double_t phi  = vtrack->Phi();
838       Int_t    sign = vtrack->Charge() > 0 ? 1 : -1;
839       if (IsTrigger(eta,pt)||IsAssociated(eta,pt))
840         miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
841     }
842   }
843
844   if (fFillMuons) {
845     // count good muons
846     Int_t nGoodMuons = 0;
847     for (Int_t iMu = 0; iMu<fAOD->GetNumberOfTracks(); iMu++) {
848       AliAODTrack* muonTrack = fAOD->GetTrack(iMu);
849       if (muonTrack) {
850         if (!IsGoodMUONtrack(*muonTrack))
851           continue;
852         Double_t ptMu   = muonTrack->Pt();
853         Double_t etaMu  = muonTrack->Eta();
854         if ((IsTrigger(etaMu,ptMu)||IsAssociated(etaMu,ptMu)) && (IsGoodMUONtrack(*muonTrack)))
855           nGoodMuons++;
856       }
857     }
858     miniEvt->reserve(miniEvt->size()+nGoodMuons);
859     // fill them into the mini event
860     for (Int_t iMu = 0; iMu<fAOD->GetNumberOfTracks(); iMu++) {
861       AliAODTrack* muonTrack = fAOD->GetTrack(iMu);
862       if (muonTrack) {
863         if (!IsGoodMUONtrack(*muonTrack)) 
864           continue;
865         Double_t ptMu   = muonTrack->Pt();
866         Double_t etaMu  = muonTrack->Eta();
867         Double_t phiMu  = muonTrack->Phi();
868         Int_t    signMu = muonTrack->Charge() > 0 ? 1 : -1;
869         if (IsTrigger(etaMu,ptMu)||IsAssociated(etaMu,ptMu))
870           miniEvt->push_back(AliMiniTrack(ptMu, etaMu, phiMu, signMu));
871       }
872     }
873   }
874 }
875
876 //________________________________________________________________________
877 Double_t AliDhcTask::DeltaPhi(Double_t phia, Double_t phib,
878                               Double_t rangeMin, Double_t rangeMax) const
879 {
880   Double_t dphi = -999;
881   Double_t pi = TMath::Pi();
882   
883   if (phia < 0)         phia += 2*pi;
884   else if (phia > 2*pi) phia -= 2*pi;
885   if (phib < 0)         phib += 2*pi;
886   else if (phib > 2*pi) phib -= 2*pi;
887   dphi = phib - phia;
888   if (dphi < rangeMin)      dphi += 2*pi;
889   else if (dphi > rangeMax) dphi -= 2*pi;
890   
891   return dphi;
892 }
893
894 //________________________________________________________________________
895 Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, Int_t pairing)
896 {
897   // Triggered angular correlations. If pairing is kSameEvt, particles
898   // within evt1 are correlated. If kDiffEvt, correlate triggers from
899   // evt1 with partners from evt2.
900
901   Int_t cbin = fHCent->FindBin(fCentrality);
902   if (fHCent->IsBinOverflow(cbin) ||
903       fHCent->IsBinUnderflow(cbin))
904     return 0;
905
906   Int_t zbin = fHZvtx->FindBin(fZVertex);
907   if (fHZvtx->IsBinOverflow(zbin) ||
908       fHZvtx->IsBinUnderflow(zbin))
909     return 0;
910
911   Int_t iMax = evt1.size();
912   Int_t jMax = evt2.size();
913
914   TH2  **hist = fHMs;
915   if (pairing == kSameEvt) {
916     hist = fHSs;
917     fHCent->Fill(fCentrality);
918     fHZvtx->Fill(fZVertex);
919   }
920
921   Int_t nZvtx = fHZvtx->GetNbinsX();
922   Int_t nPtTrig = fHPtTrg->GetNbinsX();
923   Int_t nPtAssc = fHPtAss->GetNbinsX();
924
925   Int_t globIndex = (cbin-1)*nZvtx*nPtTrig*nPtAssc+(zbin-1)*nPtTrig*nPtAssc;
926   Int_t ptindex   = (Int_t)fIndex->Eval(1,1,zbin,cbin);
927   Int_t mindex    = (Int_t)fIndex->Eval(1,1,1,cbin);
928
929
930   fHPtTrgEvt->Reset();
931   for (Int_t i=0; i<iMax; ++i) {
932     const AliMiniTrack &a(evt1.at(i));
933     Float_t pta  = a.Pt();
934     fHPtTrgEvt->Fill(pta);
935     if (pairing == kSameEvt) {
936       fHPts[ptindex]->Fill(pta);
937     }
938   }
939
940   Bool_t bCountTrg; // only count the trigger if an associated particle is found
941
942   for (Int_t i=0; i<iMax; ++i) {
943     // Trigger particles
944     const AliMiniTrack &a(evt1.at(i));
945
946     Float_t pta  = a.Pt();
947     Float_t etaa = a.Eta();
948     Float_t phia = a.Phi();
949     Int_t   sgna = a.Sign();
950     
951     // brief intermezzo in the trigger particle loop: do associated particle QA here in order to avoid double counting
952     if (pairing == kSameEvt) {
953       if (IsAssociated(etaa,pta)) {
954         Double_t aQAWght = 1.0;
955         if (fHEffA) {
956           const Int_t nEffDimA = fHEffA->GetNdimensions();
957           Int_t effBinA[nEffDimA];
958           effBinA[0] = fHEffA->GetAxis(0)->FindBin(etaa);
959           effBinA[1] = fHEffA->GetAxis(1)->FindBin(pta);
960           effBinA[2] = fHEffA->GetAxis(2)->FindBin(fCentrality);
961           effBinA[3] = fHEffA->GetAxis(3)->FindBin(fZVertex);
962           if (nEffDimA>4) {
963             effBinA[4] = fHEffA->GetAxis(4)->FindBin(phia);
964             if (nEffDimA>5) {
965               effBinA[5] = fHEffA->GetAxis(5)->FindBin(sgna);
966             }
967           }
968           aQAWght = fHEffA->GetBinContent(effBinA);
969         }
970         // fill every associated particle once unweighted, once weighted
971         if (sgna>0.0) {
972           fHQAAp->Fill(pta,etaa,phia);
973           fHQAApCorr->Fill(pta,etaa,phia,aQAWght);
974         } else {
975           fHQAAm->Fill(pta,etaa,phia);
976           fHQAAmCorr->Fill(pta,etaa,phia,aQAWght);
977         }
978         fHPtCentA->Fill(pta,fCentrality);
979       }
980     }
981     
982     // back to triggers
983     if (!IsTrigger(etaa,pta))
984       continue;
985
986     Int_t abin = fHPtAss->FindBin(pta);
987
988     // efficiency weighting
989     Double_t effWtT = 1.0;
990     if (fHEffT) {
991       // trigger particle
992       const Int_t nEffDimT = fHEffT->GetNdimensions();
993       Int_t effBinT[nEffDimT];
994       effBinT[0] = fHEffT->GetAxis(0)->FindBin(etaa);
995       effBinT[1] = fHEffT->GetAxis(1)->FindBin(pta);
996       effBinT[2] = fHEffT->GetAxis(2)->FindBin(fCentrality);
997       effBinT[3] = fHEffT->GetAxis(3)->FindBin(fZVertex);
998       if (nEffDimT>4) {
999         effBinT[4] = fHEffT->GetAxis(4)->FindBin(phia);
1000         if (nEffDimT>5) {
1001           effBinT[5] = fHEffT->GetAxis(5)->FindBin(sgna);
1002         }
1003       }
1004       effWtT = fHEffT->GetBinContent(effBinT);
1005     }
1006     
1007     if (pairing == kSameEvt) {
1008       fHTrk->Fill(phia,etaa);
1009       if (sgna>0.0) {
1010         fHQATp->Fill(pta,etaa,phia);
1011         fHQATpCorr->Fill(pta,etaa,phia,effWtT);
1012       } else {
1013         fHQATm->Fill(pta,etaa,phia);
1014         fHQATmCorr->Fill(pta,etaa,phia,effWtT);
1015       }
1016       fHPtCentT->Fill(pta,fCentrality);
1017       fHPtTrg->Fill(pta);
1018       fHPtTrgNorm1S->Fill(pta,fCentrality,fZVertex,effWtT);
1019     } else {
1020       fHPtTrgNorm1M->Fill(pta,fCentrality,fZVertex,effWtT);
1021     }
1022     bCountTrg = kFALSE;
1023
1024     for (Int_t j=0; j<jMax; ++j) {
1025       // Associated particles
1026       if (pairing == kSameEvt && i==j)
1027         continue;
1028
1029       const AliMiniTrack &b(evt2.at(j));
1030       
1031       Float_t ptb  = b.Pt();
1032       Float_t etab = b.Eta();
1033       Float_t phib = b.Phi();
1034       Int_t   sgnb = b.Sign();
1035       if (fPtTACrit&&(pta < ptb)) {
1036         continue;
1037       }
1038
1039       if (!IsAssociated(etab,ptb))
1040         continue;
1041
1042       Int_t bbin = fHPtAss->FindBin(ptb);
1043
1044       Float_t dphi = DeltaPhi(phia, phib);
1045       Float_t deta = etaa - etab;
1046       // invariant mass under mu mu assumption
1047       Float_t muMass2 = 0.1056583715*0.1056583715;
1048       Float_t ea2  = muMass2 + pta*pta*TMath::CosH(etaa)*TMath::CosH(etaa);
1049       Float_t eb2  = muMass2 + ptb*ptb*TMath::CosH(etab)*TMath::CosH(etab);
1050       Float_t papb = pta*TMath::Cos(phia)*ptb*TMath::Cos(phib) +
1051                      pta*TMath::Sin(phia)*ptb*TMath::Sin(phib) +
1052                      pta*TMath::SinH(etaa)*ptb*TMath::SinH(etab);
1053       Float_t mass = TMath::Sqrt(2.0*(muMass2 + TMath::Sqrt(ea2*eb2) - papb));
1054       Int_t q2 = sgna + sgnb;
1055       if ((q2==0) && fDoMassCut) {
1056         if (mass>3.0 && mass<3.2)
1057           continue;
1058         if (mass>9.2&&mass<9.8)
1059           continue;
1060       }
1061
1062       Int_t index = globIndex+(abin-1)*nPtAssc+(bbin-1);
1063       Double_t weight = 1.0;
1064       // number of trigger weight event-by-event (CMS method)
1065       if (fDoWeights) {
1066         Double_t nTrgs = fHPtTrgEvt->GetBinContent(abin);
1067         if (nTrgs==0.0) {
1068           AliError(Form("Trying to work with number of triggers weight = %f",nTrgs));
1069           return 0;
1070         }
1071         weight *= 1./nTrgs;
1072       }
1073
1074       // efficiency weighting
1075       if (fHEffT) {
1076         // trigger particle
1077         weight *= effWtT;
1078       }
1079
1080       if (fHEffA) {
1081         // associated particle
1082         const Int_t nEffDimA = fHEffA->GetNdimensions();
1083         Int_t effBinA[nEffDimA];
1084         effBinA[0] = fHEffA->GetAxis(0)->FindBin(etab);
1085         effBinA[1] = fHEffA->GetAxis(1)->FindBin(ptb);
1086         effBinA[2] = fHEffA->GetAxis(2)->FindBin(fCentrality);
1087         effBinA[3] = fHEffA->GetAxis(3)->FindBin(fZVertex);
1088         if (nEffDimA>4) {
1089           effBinA[4] = fHEffA->GetAxis(4)->FindBin(phib);
1090           if (nEffDimA>5) {
1091             effBinA[5] = fHEffA->GetAxis(5)->FindBin(sgnb);
1092           }
1093         }
1094         weight *= fHEffA->GetBinContent(effBinA);
1095       }
1096
1097       if (hist[index]) { // check if this histogram exists, relevant in the fPtTACrit==kFALSE case
1098         hist[index]->Fill(dphi,deta,weight);
1099         bCountTrg = kTRUE;
1100         if (pairing == kSameEvt) {
1101           fHPtAss->Fill(ptb); // fill every associated particle every time it is used
1102           if (q2==0)
1103             fHSMass[mindex]->Fill(mass);
1104         } else {
1105           if (q2==0)
1106             fHMMass[mindex]->Fill(mass);
1107         }
1108       }
1109     }
1110     if (bCountTrg) {
1111       if (pairing == kSameEvt) {
1112         fHPtTrgNorm2S->Fill(pta,fCentrality,fZVertex,effWtT);
1113       } else {
1114         fHPtTrgNorm2M->Fill(pta,fCentrality,fZVertex,effWtT);
1115       }
1116     }
1117   }
1118
1119   return 1;
1120 }
1121
1122 //________________________________________________________________________
1123 void AliDhcTask::Terminate(Option_t *) 
1124 {
1125   // Draw result to the screen
1126   // Called once at the end of the query
1127 }
1128
1129 //________________________________________________________________________
1130 Bool_t AliDhcTask::VertexOk() const
1131 {
1132   // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
1133   
1134   Int_t nContributors  = 0;
1135   Double_t zVertex     = 999;
1136   TString name("");
1137
1138   Int_t runno = InputEvent()->GetRunNumber();
1139   if (runno>=176326 && runno<=197692) { // year 12 and 13
1140     if (!fUtils->IsVertexSelected2013pA(InputEvent())) 
1141       return 0;
1142   }
1143
1144   if (fESD) {
1145     const AliESDVertex* vtx = fESD->GetPrimaryVertex();
1146     if (!vtx)
1147       return 0;
1148     nContributors = vtx->GetNContributors();
1149     zVertex       = vtx->GetZ();
1150     name          = vtx->GetName();
1151   } else {
1152     if (fAOD->GetNumberOfVertices() < 1)
1153       return 0;
1154     const AliAODVertex* vtx = fAOD->GetPrimaryVertex();
1155     nContributors = vtx->GetNContributors();
1156     zVertex       = vtx->GetZ();
1157     name          = vtx->GetName();
1158   }
1159   
1160   // Reject if TPC-only vertex
1161   if (name.CompareTo("TPCVertex")==0)
1162     return kFALSE;
1163   
1164   // Check # contributors and range...
1165   if( nContributors < 1 || TMath::Abs(zVertex) > fZVtxMax ) {
1166     return kFALSE;
1167   }
1168   
1169   return kTRUE;
1170 }
1171
1172 //________________________________________________________________________
1173 Bool_t AliDhcTask::IsGoodMUONtrack(AliESDMuonTrack &track)
1174 {
1175   // Applying track cuts for MUON tracks
1176
1177   if (!track.ContainTrackerData()) 
1178     return kFALSE;
1179   Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
1180   if ((thetaTrackAbsEnd < 2.) || (thetaTrackAbsEnd > 10.)) 
1181     return kFALSE;
1182   Double_t eta = track.Eta();
1183   if ((eta < -4.) || (eta > -2.5))
1184     return kFALSE;
1185   if (fTriggerMatch) {
1186     if (!track.ContainTriggerData()) 
1187       return kFALSE;
1188     if (track.GetMatchTrigger() < 0.5) 
1189       return kFALSE;
1190   }
1191   return kTRUE;
1192 }
1193
1194 //________________________________________________________________________
1195 Bool_t AliDhcTask::IsGoodMUONtrack(AliAODTrack &track)
1196 {
1197   // Applying track cuts for MUON tracks
1198
1199   if (!track.IsMuonTrack()) 
1200     return kFALSE;
1201   Double_t dThetaAbs = TMath::ATan(track.GetRAtAbsorberEnd()/505.)* TMath::RadToDeg();
1202   if ((dThetaAbs<2.) || (dThetaAbs>10.)) 
1203     return kFALSE;
1204   Double_t dEta = track.Eta();
1205   if ((dEta<-4.) || (dEta>-2.5)) 
1206     return kFALSE;
1207   if (fTriggerMatch) {
1208     if (track.GetMatchTrigger()<0.5) 
1209       return kFALSE;
1210   }
1211   return kTRUE;
1212 }
1213
1214 //________________________________________________________________________
1215 Bool_t AliDhcTask::IsTrigger(Double_t eta, Double_t pt)
1216 {
1217   // is in trigger eta,pt range?
1218   Int_t bin = fHPtTrg->FindBin(pt);
1219   if (fHPtTrg->IsBinOverflow(bin) || fHPtTrg->IsBinUnderflow(bin))
1220     return kFALSE;
1221   
1222   if (eta<fEtaTLo || eta>fEtaTHi)
1223     return kFALSE;
1224   
1225   return kTRUE;
1226 }
1227
1228 //________________________________________________________________________
1229 Bool_t AliDhcTask::IsAssociated(Double_t eta, Double_t pt)
1230 {
1231   // is in associated eta,pt range?
1232   Int_t bin = fHPtAss->FindBin(pt);
1233   if (fHPtAss->IsBinOverflow(bin) || fHPtAss->IsBinUnderflow(bin))
1234     return kFALSE;
1235   
1236   if (eta<fEtaALo || eta>fEtaAHi)
1237     return kFALSE;
1238   
1239   return kTRUE;
1240 }
1241
1242 //________________________________________________________________________
1243 AliDhcTask::~AliDhcTask()
1244 {
1245   //Destructor
1246   if (fPoolMgr) 
1247     delete fPoolMgr;
1248 }