AliPhysics  b3a51e4 (b3a51e4)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskCLQA.cxx
Go to the documentation of this file.
1 // $Id: $
2 //
3 // Constantin's Task
4 //
5 // Author: C.Loizides
6 
7 #include <complex>
8 
9 #include <TChain.h>
10 #include <TClonesArray.h>
11 #include <TDirectory.h>
12 #include <TFile.h>
13 #include <TH1F.h>
14 #include <TH2F.h>
15 #include <TH3F.h>
16 #include <TList.h>
17 #include <TLorentzVector.h>
18 #include <TNtuple.h>
19 #include <TNtupleD.h>
20 #include <TProfile.h>
21 #include <TTree.h>
22 #include "AliESDMuonTrack.h"
23 #include "AliAODEvent.h"
24 #include "AliAnalysisManager.h"
25 #include "AliAnalysisUtils.h"
26 #include "AliCentrality.h"
27 #include "AliEMCALGeoParams.h"
28 #include "AliEMCALGeometry.h"
29 #include "AliEMCALRecoUtils.h"
30 #include "AliESDEvent.h"
31 #include "AliEmcalJet.h"
32 #include "AliExternalTrackParam.h"
33 #include "AliInputEventHandler.h"
34 #include "AliLog.h"
35 #include "AliPicoTrack.h"
36 #include "AliTrackerBase.h"
37 #include "AliVCluster.h"
38 #include "AliVEventHandler.h"
39 #include "AliVParticle.h"
40 #include "AliVTrack.h"
41 #include "AliAnalysisTaskCLQA.h"
42 
44 
45 //________________________________________________________________________
48  fDo2013VertexCut(1),
49  fDoTracking(1), fDoMuonTracking(0), fDoCumulants(0), fDoCumNtuple(0),
50  fCumPtMin(0.3), fCumPtMax(5.0), fCumEtaMin(-1.0), fCumEtaMax(1.0), fCumMmin(15),
51  fCumMbins(250), fCentCL1In(0), fCentV0AIn(0),
52  fNtupCum(0), fNtupCumInfo(0), fNtupZdcInfo(0)
53 {
54  // Default constructor.
55 
56  for (Int_t i=0;i<1000;++i)
57  fHists[i] = 0;
58 }
59 
60 //________________________________________________________________________
62  AliAnalysisTaskEmcal(name, kTRUE),
63  fDo2013VertexCut(1),
64  fDoTracking(1), fDoMuonTracking(0), fDoCumulants(0), fDoCumNtuple(0),
65  fCumPtMin(0.3), fCumPtMax(5.0), fCumEtaMin(-1.0), fCumEtaMax(1.0), fCumMmin(15),
66  fCumMbins(250), fCentCL1In(0), fCentV0AIn(0),
67  fNtupCum(0), fNtupCumInfo(0), fNtupZdcInfo(0)
68 {
69  // Standard constructor.
70 
71  for (Int_t i=0;i<1000;++i)
72  fHists[i] = 0;
73 }
74 
75 //________________________________________________________________________
77 {
78  // Destructor
79 }
80 
81 //________________________________________________________________________
82 Double_t AliAnalysisTaskCLQA::DeltaPhi(Double_t phia, Double_t phib,
83  Double_t rangeMin, Double_t rangeMax) const
84 {
85  // Calculate Delta Phi.
86 
87  Double_t dphi = -999;
88  Double_t pi = TMath::Pi();
89 
90  if (phia < 0) phia += 2*pi;
91  else if (phia > 2*pi) phia -= 2*pi;
92  if (phib < 0) phib += 2*pi;
93  else if (phib > 2*pi) phib -= 2*pi;
94  dphi = phib - phia;
95  if (dphi < rangeMin) dphi += 2*pi;
96  else if (dphi > rangeMax) dphi -= 2*pi;
97 
98  return dphi;
99 }
100 
101 //________________________________________________________________________
103 {
104  // Fill histograms.
105 
106  AliVEvent *event = InputEvent();
107  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
108 
109  UInt_t trig = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
110  for (Int_t i=0;i<31;++i) {
111  if (trig & (1<<i))
112  fHists[0]->Fill(i);
113  }
114 
115  Double_t vz = event->GetPrimaryVertex()->GetZ();
116  fHists[1]->Fill(vz);
117 
118  Int_t run = event->GetRunNumber();
119  Int_t vzn = event->GetPrimaryVertex()->GetNContributors();
120  if ((vzn<1)&&(run>0))
121  return kTRUE;
122  fHists[2]->Fill(vz);
123 
124  if (TMath::Abs(vz)>10)
125  return kTRUE;
126 
127  if (fDo2013VertexCut) {
128  if ((run>=188356&&run<=188366) || (run>=195344&&run<=197388)) {
129  AliAnalysisUtils anau;
130  if (anau.IsFirstEventInChunk(event))
131  return kFALSE;
132  if (!anau.IsVertexSelected2013pA(event))
133  return kFALSE;
134  }
135  }
136 
137  // accepted events
138  fHists[9]->Fill(1,run);
139 
140  AliCentrality *cent = InputEvent()->GetCentrality();
141  Double_t v0acent = cent->GetCentralityPercentile("V0A");
142  fHists[10]->Fill(v0acent);
143  Double_t znacent = cent->GetCentralityPercentile("ZNA");
144  fHists[11]->Fill(znacent);
145  Double_t v0mcent = cent->GetCentralityPercentile("V0M");
146  fHists[12]->Fill(v0mcent);
147 
148  if (fDoTracking) {
149  const Int_t ntracks = fTracks->GetEntries();
150  if (fTracks) {
151  for (Int_t i=0; i<ntracks; ++i) {
152  AliVTrack *track = dynamic_cast<AliVTrack*>(fTracks->At(i));
153  if (!track)
154  continue;
155  if (track->Charge()==0)
156  continue;
157 
158  AliPicoTrack *picot = dynamic_cast<AliPicoTrack*>(track);
159  if (picot && picot->GetTrack())
160  track = picot->GetTrack();
161 
162  Double_t phi = track->Phi();
163  Double_t eta = track->Eta();
164  Double_t pt = track->Pt();
165  fHists[20]->Fill(phi,eta);
166  fHists[21]->Fill(phi,pt);
167  fHists[22]->Fill(eta,pt);
168  if (TMath::Abs(eta)<0.8) {
169  fHists[23]->Fill(pt,v0acent);
170  fHists[24]->Fill(pt,znacent);
171  }
172  Int_t ttype = AliPicoTrack::GetTrackType(track);
173  fHists[25+ttype]->Fill(phi,pt);
174  if (track->IsExtrapolatedToEMCAL()) {
175  Double_t dphi = TVector2::Phi_mpi_pi(phi-track->GetTrackPhiOnEMCal());
176  fHists[28]->Fill(dphi,pt);
177  } else {
178  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(track,440);
179  if (track->IsExtrapolatedToEMCAL()) {
180  Double_t dphi = TVector2::Phi_mpi_pi(phi-track->GetTrackPhiOnEMCal());
181  fHists[29]->Fill(dphi,pt);
182  }
183  }
184  if (track->IsEMCAL() && track->IsExtrapolatedToEMCAL()) {
185  Int_t id = track->GetEMCALcluster();
186  AliVCluster *clus = InputEvent()->GetCaloCluster(id);
187  if (id>=0&&clus) {
188  Float_t pos[3];
189  clus->GetPosition(pos);
190  TVector3 vpos(pos);
191  Double_t dphi = TVector2::Phi_mpi_pi(vpos.Phi()-track->GetTrackPhiOnEMCal());
192  fHists[30]->Fill(dphi,pt);
193  }
194  }
195  if (track->IsExtrapolatedToEMCAL()) {
196  Double_t phi1 = track->GetTrackPhiOnEMCal();
197  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(track,440);
198  Double_t phi2 = track->GetTrackPhiOnEMCal();
199  Double_t dphi = TVector2::Phi_mpi_pi(phi1-phi2);
200  fHists[31]->Fill(dphi,pt);
201  }
202  }
203  }
204  }
205 
206  if (fDoMuonTracking) {
207  AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
208  if (aod) {
209  for (Int_t iMu = 0; iMu<aod->GetNumberOfTracks(); ++iMu) {
210  AliAODTrack* muonTrack = static_cast<AliAODTrack*>(aod->GetTrack(iMu));
211  if (!muonTrack)
212  continue;
213  if (!muonTrack->IsMuonTrack())
214  continue;
215  Double_t dThetaAbs = TMath::ATan(muonTrack->GetRAtAbsorberEnd()/505.)* TMath::RadToDeg();
216  if ((dThetaAbs<2.) || (dThetaAbs>10.))
217  continue;
218  Double_t dEta = muonTrack->Eta();
219  if ((dEta<-4.) || (dEta>-2.5))
220  continue;
221  if (0) {
222  if (muonTrack->GetMatchTrigger()<0.5)
223  continue;
224  }
225  Double_t ptMu = muonTrack->Pt();
226  Double_t etaMu = muonTrack->Eta();
227  Double_t phiMu = muonTrack->Phi();
228  fHists[50]->Fill(phiMu,etaMu);
229  fHists[51]->Fill(phiMu,ptMu);
230  fHists[52]->Fill(etaMu,ptMu);
231  fHists[53]->Fill(ptMu,v0acent);
232  fHists[54]->Fill(ptMu,znacent);
233  }
234  } else {
235  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
236  if (esd) {
237  for (Int_t iMu = 0; iMu<esd->GetNumberOfMuonTracks(); ++iMu) {
238  AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iMu);
239  if (!muonTrack)
240  continue;
241  if (!muonTrack->ContainTrackerData())
242  continue;
243  Double_t thetaTrackAbsEnd = TMath::ATan(muonTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
244  if ((thetaTrackAbsEnd < 2.) || (thetaTrackAbsEnd > 10.))
245  continue;
246  Double_t eta = muonTrack->Eta();
247  if ((eta < -4.) || (eta > -2.5))
248  return kFALSE;
249  if (0) {
250  if (!muonTrack->ContainTriggerData())
251  continue;
252  if (muonTrack->GetMatchTrigger() < 0.5)
253  continue;
254  }
255  }
256  }
257  }
258  }
259 
260  return kTRUE;
261 }
262 
263 //________________________________________________________________________
265 {
266  // Retrieve event objects.
267 
269  return kFALSE;
270 
271  return kTRUE;
272 }
273 
274 //________________________________________________________________________
276 {
277  // Run various functions.
278 
280 
281  return kTRUE;
282 }
283 
284 //________________________________________________________________________
285 void AliAnalysisTaskCLQA::RunCumulants(Double_t Mmin, Double_t ptmin, Double_t ptmax, Double_t etamin, Double_t etamax)
286 {
287  // Run cumulant analysis.
288 
289  if (!fDoCumulants)
290  return;
291 
292  if (!fTracks)
293  return;
294 
295  Bool_t isMC = 0;
296  TString tname(fTracks->GetName());
297  if (tname.Contains("mc"))
298  isMC = 1;
299 
300  const Int_t ntracks = fTracks->GetEntries();
301  Int_t Mall=0,M=0,Mall2=0;
302  Double_t ptmaxall=0,ptsumall=0,pt2sumall=0,ptsumall2=0;
303  Double_t tsa00=0,tsa10=0,tsa11=0;
304  Double_t Q2r=0,Q2i=0;
305  Double_t Q3r=0,Q3i=0;
306  Double_t Q4r=0,Q4i=0;
307  Double_t Q6r=0,Q6i=0;
308  Double_t mpt=0,mpt2=0,ptmaxq=0;
309  Double_t ts00=0,ts10=0,ts11=0;
310  Double_t v0ach=0, v0cch=0;
311  Double_t cl1ch=0;
312 
313  for (Int_t i =0; i<ntracks; ++i) {
314  AliVParticle *track = dynamic_cast<AliVParticle*>(fTracks->At(i));
315  if (!track)
316  continue;
317  if (track->Charge()==0)
318  continue;
319  Double_t eta = track->Eta();
320  if ((eta<5.1)&&(eta>2.8))
321  ++v0ach;
322  else if ((eta>-3.7)&&(eta<-1.7))
323  ++v0cch;
324  if (TMath::Abs(eta)<1.4) {
325  ++cl1ch;
326  }
327  if ((eta<etamin) || (eta>etamax))
328  continue;
329  Double_t pt = track->Pt();
330  if (pt>ptmaxall)
331  ptmaxall = pt;
332  if (pt>1) {
333  ptsumall2 += pt;
334  ++Mall2;
335  }
336  ptsumall +=pt;
337  pt2sumall +=pt*pt;
338  Double_t px = track->Px();
339  Double_t py = track->Py();
340  tsa00 += px*px/pt;
341  tsa10 += px*py/pt;
342  tsa11 += py*py/pt;
343  ++Mall;
344  if ((pt<ptmin) || (pt>ptmax))
345  continue;
346  if (pt>ptmaxq)
347  ptmaxq = pt;
348  Double_t phi = track->Phi();
349  ++M;
350  mpt += pt;
351  mpt2 += pt*pt;
352  ts00 += px*px/pt;
353  ts10 += px*py/pt;
354  ts11 += py*py/pt;
355  Q2r += TMath::Cos(2*phi);
356  Q2i += TMath::Sin(2*phi);
357  Q3r += TMath::Cos(3*phi);
358  Q3i += TMath::Sin(3*phi);
359  Q4r += TMath::Cos(4*phi);
360  Q4i += TMath::Sin(4*phi);
361  Q6r += TMath::Cos(6*phi);
362  Q6i += TMath::Sin(6*phi);
363  }
364 
365  if (M<=1)
366  return;
367 
368  Int_t pmult=0;
369  Double_t v2g=0;
370  Double_t v3g=0;
371  Int_t pmult14=0;
372  Double_t v2g14=0;
373  Double_t v3g14=0;
374  Int_t pmult18=0;
375  Double_t v2g18=0;
376  Double_t v3g18=0;
377  for (Int_t i=0; i<ntracks; ++i) {
378  AliVParticle *track1 = dynamic_cast<AliVParticle*>(fTracks->At(i));
379  if (!track1)
380  continue;
381  if (track1->Charge()==0)
382  continue;
383  Double_t eta1 = track1->Eta();
384  if ((eta1<etamin) || (eta1>etamax))
385  continue;
386  Double_t pt1 = track1->Pt();
387  if ((pt1<ptmin) || (pt1>ptmax))
388  continue;
389  Double_t phi1 = track1->Phi();
390  for (Int_t j = i+1; j<ntracks; ++j) {
391  AliVParticle *track2 = dynamic_cast<AliVParticle*>(fTracks->At(j));
392  if (!track2)
393  continue;
394  if (track2->Charge()==0)
395  continue;
396  Double_t eta2 = track2->Eta();
397  if ((eta2<etamin) || (eta2>etamax))
398  continue;
399  Double_t pt2 = track2->Pt();
400  if ((pt2<ptmin) || (pt2>ptmax))
401  continue;
402  Double_t deta=TMath::Abs(eta2-eta1);
403  ((TH3*)fHists[128])->Fill(DeltaPhi(phi1,track2->Phi()),deta,M);
404  fHists[129]->Fill(M);
405  if(deta<1)
406  continue;
407  Double_t dphi=TVector2::Phi_0_2pi(phi1-track2->Phi());
408  pmult++;
409  v2g+=TMath::Cos(2*dphi);
410  v3g+=TMath::Cos(3*dphi);
411  if (deta>1.4) {
412  pmult14++;
413  v2g14+=TMath::Cos(2*dphi);
414  v3g14+=TMath::Cos(3*dphi);
415  }
416  if (deta>1.8) {
417  pmult18++;
418  v2g18+=TMath::Cos(2*dphi);
419  v3g18+=TMath::Cos(3*dphi);
420  }
421  }
422  }
423 
424  if (pmult>0) {
425  v2g/=pmult;
426  v3g/=pmult;
427  }
428  if (pmult14>0) {
429  v2g14/=pmult14;
430  v3g14/=pmult14;
431  }
432  if (pmult18>0) {
433  v2g18/=pmult18;
434  v3g18/=pmult18;
435  }
436 
437  std::complex<double> q2(Q2r,Q2i);
438  std::complex<double> q3(Q3r,Q3i);
439  std::complex<double> q4(Q4r,Q4i);
440  std::complex<double> q6(Q6r,Q6i);
441  Double_t Q22 = std::abs(q2)*std::abs(q2);
442  Double_t Q32 = std::abs(q3)*std::abs(q3);
443  Double_t Q42 = std::abs(q4)*std::abs(q4);
444  Double_t Q62 = std::abs(q6)*std::abs(q6);
445  Double_t Q32re = std::real(q6*std::conj(q3)*std::conj(q3));
446  Double_t Q42re = std::real(q4*std::conj(q2)*std::conj(q2));
447  Double_t Q6are = std::real(q4*q2*std::conj(q2)*std::conj(q2)*std::conj(q2));
448  Double_t Q6bre = std::real(q6*std::conj(q2)*std::conj(q2)*std::conj(q2));
449  Double_t Q6cre = std::real(q6*std::conj(q4)*std::conj(q2));
450 
451  Double_t tsall = -1;
452  Double_t tsax = (tsa00+tsa11)*(tsa00+tsa11)-4*(tsa00*tsa11-tsa10*tsa10);
453  if (tsax>=0) {
454  Double_t l1 = 0.5*(tsa00+tsa11+TMath::Sqrt(tsax))/ptsumall;
455  Double_t l2 = 0.5*(tsa00+tsa11-TMath::Sqrt(tsax))/ptsumall;
456  tsall = 2*l2/(l1+l2);
457  }
458 
459  Double_t ts = -1;
460  Double_t tsx = (ts00+ts11)*(ts00+ts11)-4*(ts00*ts11-ts10*ts10);
461  if (tsx>=0) {
462  Double_t l1 = 0.5*(ts00+ts11+TMath::Sqrt(tsx))/ptsumall;
463  Double_t l2 = 0.5*(ts00+ts11-TMath::Sqrt(tsx))/ptsumall;
464  ts = 2*l2/(l1+l2);
465  }
466 
467  if (isMC) {
468  fHists[106]->Fill(cl1ch);
469  fHists[107]->Fill(v0ach);
470  fHists[108]->Fill(v0cch);
471  AliCentrality *cent = InputEvent()->GetCentrality();
472  if (fCentCL1In) {
473  cent->SetCentralityCL1(100*fCentCL1In->GetBinContent(fCentCL1In->FindBin(cl1ch)));
474  cent->SetQuality(0);
475  }
476  if (fCentV0AIn) {
477  cent->SetCentralityV0A(100*fCentV0AIn->GetBinContent(fCentV0AIn->FindBin(v0ach)));
478  cent->SetQuality(0);
479  }
480  }
481 
482  AliAnalysisUtils anau;
483  AliVEvent *event = InputEvent();
484  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
485 
486  fNtupCumInfo->fTrig = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
487  fNtupCumInfo->fRun = event->GetRunNumber();
488  fNtupCumInfo->fVz = event->GetPrimaryVertex()->GetZ();
489  fNtupCumInfo->fIsFEC = anau.IsFirstEventInChunk(event);
490  fNtupCumInfo->fIsVSel = anau.IsVertexSelected2013pA(event);
491  fNtupCumInfo->fIsP = event->IsPileupFromSPD(3/*minContributors*/,
492  0.8/*minZdist*/,
493  3./*nSigmaZdist*/,
494  2./*nSigmaDiamXY*/,
495  5./*nSigmaDiamZ*/);
496 
497  fNtupCumInfo->fMall = Mall;
498  fNtupCumInfo->fMall2 = Mall2;
499  fNtupCumInfo->fPtMaxall = ptmaxall;
500  fNtupCumInfo->fMPtall = ptsumall/Mall;
501  fNtupCumInfo->fMPt2all = pt2sumall/Mall;
502  if (Mall2>0)
503  fNtupCumInfo->fMPtall2 = ptsumall2/Mall2;
504  else
505  fNtupCumInfo->fMPtall2 = -1;
506  fNtupCumInfo->fTSall = tsall;
507  fNtupCumInfo->fM = M;
508  fNtupCumInfo->fQ2abs = Q22;
509  fNtupCumInfo->fQ4abs = Q42;
510  fNtupCumInfo->fQ42re = Q42re;
511  fNtupCumInfo->fCos2phi = Q2r;
512  fNtupCumInfo->fSin2phi = Q2i;
513  fNtupCumInfo->fPtMax = ptmaxq;
514  fNtupCumInfo->fMPt = mpt/M;
515  fNtupCumInfo->fMPt2 = mpt2/M;
516  fNtupCumInfo->fTS = ts;
517 
518  if (isMC) {
519  fNtupCumInfo->fMV0M = v0ach + v0cch;
520  } else {
521  AliVVZERO *vzero = InputEvent()->GetVZEROData();
522  fNtupCumInfo->fMV0M = vzero->GetMTotV0A()+vzero->GetMTotV0C();
523  }
524 
525  AliCentrality *cent = InputEvent()->GetCentrality();
526  fNtupCumInfo->fCl1 = cent->GetCentralityPercentile("CL1");
527  fNtupCumInfo->fV0M = cent->GetCentralityPercentile("V0M");
528  fNtupCumInfo->fV0MEq = cent->GetCentralityPercentile("V0MEq");
529  fNtupCumInfo->fV0A = cent->GetCentralityPercentile("V0A");
530  fNtupCumInfo->fV0AEq = cent->GetCentralityPercentile("V0AEq");
531  fNtupCumInfo->fZNA = cent->GetCentralityPercentile("ZNA");
532 
533  AliVZDC *vZDC = InputEvent()->GetZDCData();
534  const Double_t *znaTowers = vZDC->GetZNATowerEnergy();
535  fNtupZdcInfo->fZna0 = znaTowers[0];
536  fNtupZdcInfo->fZna1 = znaTowers[1];
537  fNtupZdcInfo->fZna2 = znaTowers[2];
538  fNtupZdcInfo->fZna3 = znaTowers[3];
539  fNtupZdcInfo->fZna4 = znaTowers[4];
540 
541  if (fDoCumNtuple && (M>=Mmin)) {
542  fNtupCum->Fill();
543  }
544 
545  fHists[109]->Fill(fNtupCumInfo->fCl1);
546  fHists[110]->Fill(fNtupCumInfo->fV0A);
547  fHists[111]->Fill(fNtupCumInfo->fZNA);
548 
549  Bool_t fillCumHist = kTRUE;
550  if (fillCumHist) {
551  Int_t run = InputEvent()->GetRunNumber();
552  if (fDo2013VertexCut) {
553  if ((run>=188356&&run<=188366) || (run>=195344&&run<=197388)) {
554  if (anau.IsFirstEventInChunk(event))
555  fillCumHist = kFALSE;
556  if (!anau.IsVertexSelected2013pA(event))
557  fillCumHist = kFALSE;
558  }
559  }
560  Double_t vz = InputEvent()->GetPrimaryVertex()->GetZ();
561  if (TMath::Abs(vz)>10)
562  fillCumHist = kFALSE;
563  }
564  if (fillCumHist) {
565  AliVVZERO *vzero = InputEvent()->GetVZEROData();
566  Double_t v0a = vzero->GetMTotV0A();
567  Double_t v0c = vzero->GetMTotV0C();
568  Double_t v0m = vzero->GetMTotV0A()+vzero->GetMTotV0C();
569  fHists[112]->Fill(Mall);
570  fHists[113]->Fill(M);
571  fHists[117]->Fill(v0a,M);
572  fHists[118]->Fill(v0c,M);
573  fHists[119]->Fill(v0m,M);
574  if (M>1) {
575  fHists[114]->Fill(M,(Q22-M)/M/(M-1));
576  fHists[120]->Fill(M,(Q32-M)/M/(M-1));
577  fHists[122]->Fill(M,v2g);
578  fHists[123]->Fill(M,v3g);
579  fHists[124]->Fill(M,v2g14);
580  fHists[125]->Fill(M,v3g14);
581  fHists[126]->Fill(M,v2g18);
582  fHists[127]->Fill(M,v3g18);
583  }
584  if (M>3) {
585  Double_t qc4tmp = (Q22*Q22+Q42-2*Q42re-4*(M-2)*Q22+2*M*(M-3));
586  fHists[115]->Fill(M,qc4tmp/M/(M-1)/(M-2)/(M-3));
587  qc4tmp = (Q32*Q32+Q62-2*Q32re-4*(M-2)*Q32+2*M*(M-3));
588  fHists[121]->Fill(M,qc4tmp/M/(M-1)/(M-2)/(M-3));
589  }
590  if (M>5) {
591  Double_t qc6tmp = Q22*Q22*Q22 + 9*Q42*Q22 - 6*Q6are
592  + 4*Q6bre - 12*Q6cre
593  + 18*(M-4)*Q42re + 4*Q62
594  - 9*(M-4)*Q22*Q22 - 9*(M-4)*Q42
595  + 18*(M-2)*(M-5)*Q22
596  - 6*M*(M-4)*(M-5);
597  fHists[116]->Fill(M,qc6tmp/M/(M-1)/(M-2)/(M-3)/(M-4)/(M-5));
598  }
599  }
600 
601  if ((isMC) ||
602  ((TMath::Abs(fNtupCumInfo->fVz)<10) && !fNtupCumInfo->fIsFEC && fNtupCumInfo->fIsVSel)) {
603  for (Int_t i =0; i<ntracks; ++i) {
604  AliVParticle *track1 = dynamic_cast<AliVParticle*>(fTracks->At(i));
605  if (!track1)
606  continue;
607  Double_t phi1 = track1->Phi();
608  Double_t eta1 = track1->Eta();
609  Double_t pt1 = track1->Pt();
610  ((TH3*)fHists[103])->Fill(pt1,eta1,fNtupCumInfo->fCl1);
611  ((TH3*)fHists[104])->Fill(pt1,eta1,fNtupCumInfo->fV0A);
612  ((TH3*)fHists[105])->Fill(pt1,eta1,fNtupCumInfo->fZNA);
613  if ((eta1<etamin) || (eta1>etamax))
614  continue;
615  if ((pt1<ptmin) || (pt1>ptmax))
616  continue;
617  for (Int_t j =0; j<ntracks; ++j) {
618  AliVParticle *track2 = dynamic_cast<AliVParticle*>(fTracks->At(j));
619  if (!track2)
620  continue;
621  Double_t eta2 = track2->Eta();
622  if ((eta2<etamin) || (eta2>etamax))
623  continue;
624  Double_t pt2 = track2->Pt();
625  if ((pt2<ptmin) || (pt2>ptmax))
626  continue;
627  Double_t phi2 = track2->Phi();
628  Double_t deta = eta1-eta2;
629  Double_t dphi = phi1-phi2;
630  while (dphi<-TMath::Pi())
631  dphi+=TMath::TwoPi();
632  while (dphi>3*TMath::Pi()/2)
633  dphi-=TMath::TwoPi();
634  ((TH3*)fHists[100])->Fill(dphi,deta,fNtupCumInfo->fCl1);
635  ((TH3*)fHists[101])->Fill(dphi,deta,fNtupCumInfo->fV0A);
636  ((TH3*)fHists[102])->Fill(dphi,deta,fNtupCumInfo->fZNA);
637  }
638  }
639  }
640 }
641 
642 //________________________________________________________________________
643 void AliAnalysisTaskCLQA::SetCumParams(Double_t Mmin, Double_t ptmin, Double_t ptmax, Double_t etamin, Double_t etamax)
644 {
645  // Set parameters for cumulants.
646 
647  fCumMmin = Mmin;
648  fCumPtMin = ptmin;
649  fCumPtMax = ptmax;
650  fCumEtaMin = etamin;
651  fCumEtaMax = etamax;
652 }
653 
654 //________________________________________________________________________
656 {
657  // Create histograms
658 
660 
661  fHists[0] = new TH1D("fTrigBits",";bit",32,-0.5,31.5);
662  fOutput->Add(fHists[0]);
663  fHists[1] = new TH1D("fVertexZ",";vertex z (cm)",51,-25.5,25.5);
664  fOutput->Add(fHists[1]);
665  fHists[2] = new TH1D("fVertexZnc",";vertex z (cm)",51,-25.5,25.5);
666  fOutput->Add(fHists[2]);
667  fHists[9] = new TProfile("fAccepted","",1,0.5,1.5);
668  fOutput->Add(fHists[9]);
669  fHists[10] = new TH1D("fV0ACent",";percentile",20,0,100);
670  fOutput->Add(fHists[10]);
671  fHists[11] = new TH1D("fZNACent",";percentile",20,0,100);
672  fOutput->Add(fHists[11]);
673  fHists[12] = new TH1D("fV0MCent",";percentile",20,0,100);
674  fOutput->Add(fHists[12]);
675 
676  if (fDoTracking) {
677  fHists[20] = new TH2D("fPhiEtaTracks",";#phi;#eta",60,0,TMath::TwoPi(),20,-2,2);
678  fHists[20]->Sumw2();
679  fOutput->Add(fHists[20]);
680  fHists[21] = new TH2D("fPhiPtTracks",";#phi;p_{T} (GeV/c)",60,0,TMath::TwoPi(),40,0,20);
681  fHists[21]->Sumw2();
682  fOutput->Add(fHists[21]);
683  fHists[22] = new TH2D("fEtaPtTracks",";#eta;p_{T} (GeV/c)",20,-2,2,40,0,20);
684  fHists[22]->Sumw2();
685  fOutput->Add(fHists[22]);
686  fHists[23] = new TH2D("fPtV0ATracks",";#p_{T} (GeV/c);percentile",100,0,20,20,0,100);
687  fHists[23]->Sumw2();
688  fOutput->Add(fHists[23]);
689  fHists[24] = new TH2D("fPtZNATracks",";#p_{T} (GeV/c);percentile",100,0,20,20,0,100);
690  fHists[24]->Sumw2();
691  fOutput->Add(fHists[24]);
692  fHists[25] = new TH2D("fPhiPtTracks_type0",";#phi;p_{T} (GeV/c)",60,0,TMath::TwoPi(),40,0,20);
693  fHists[25]->Sumw2();
694  fOutput->Add(fHists[25]);
695  fHists[26] = new TH2D("fPhiPtTracks_type1",";#phi;p_{T} (GeV/c)",60,0,TMath::TwoPi(),40,0,20);
696  fHists[26]->Sumw2();
697  fOutput->Add(fHists[26]);
698  fHists[27] = new TH2D("fPhiPtTracks_type2",";#phi;p_{T} (GeV/c)",60,0,TMath::TwoPi(),40,0,20);
699  fHists[27]->Sumw2();
700  fOutput->Add(fHists[27]);
701  fHists[28] = new TH2D("fDPhiPtTracks",";#Delta#phi;p_{T} (GeV/c)",60,-TMath::Pi(),TMath::Pi(),40,0,20);
702  fHists[28]->Sumw2();
703  fOutput->Add(fHists[28]);
704  fHists[29] = new TH2D("fDPhiPtTracks2",";#Delta#phi;p_{T} (GeV/c)",60,-TMath::Pi(),TMath::Pi(),40,0,20);
705  fHists[29]->Sumw2();
706  fOutput->Add(fHists[29]);
707  fHists[30] = new TH2D("fDPhiPtClusTracks",";#Delta#phi;p_{T} (GeV/c)",60,-TMath::Pi(),TMath::Pi(),40,0,20);
708  fHists[30]->Sumw2();
709  fOutput->Add(fHists[30]);
710  fHists[31] = new TH2D("fDPhiPtReTracks",";#Delta#phi;p_{T} (GeV/c)",60,-TMath::Pi(),TMath::Pi(),40,0,20);
711  fHists[31]->Sumw2();
712  fOutput->Add(fHists[31]);
713  }
714 
715  if (fDoMuonTracking) {
716  fHists[50] = new TH2D("fPhiEtaMuonTracks",";#phi;#eta",60,0,TMath::TwoPi(),15,-4,-2.5);
717  fHists[50]->Sumw2();
718  fOutput->Add(fHists[50]);
719  fHists[51] = new TH2D("fPhiPtMuonTracks",";#phi;p_{T} (GeV/c)",60,0,TMath::TwoPi(),200,0,20);
720  fHists[51]->Sumw2();
721  fOutput->Add(fHists[51]);
722  fHists[52] = new TH2D("fEtaPtMuonTracks",";#eta;p_{T} (GeV/c)",15,-4,-2.5,200,0,20);
723  fHists[52]->Sumw2();
724  fOutput->Add(fHists[52]);
725  fHists[53] = new TH2D("fPtV0AMuonTracks",";#p_{T} (GeV/c);percentile",100,0,10,20,0,100);
726  fHists[53]->Sumw2();
727  fOutput->Add(fHists[53]);
728  fHists[54] = new TH2D("fPtZNAMuonTracks",";#p_{T} (GeV/c);percentile",100,0,10,20,0,100);
729  fHists[54]->Sumw2();
730  fOutput->Add(fHists[54]);
731  }
732 
733  if (fDoCumulants) {
734  fHists[100] = new TH3D("fCumPhiEtaCl1",";#Delta#phi;#Delta#eta",32,-TMath::Pi()/2,3*TMath::Pi()/2,60,-3,3,10,0,100);
735  fOutput->Add(fHists[100]);
736  fHists[101] = new TH3D("fCumPhiEtaV0A",";#Delta#phi;#Delta#eta",32,-TMath::Pi()/2,3*TMath::Pi()/2,60,-3,3,10,0,100);
737  fOutput->Add(fHists[101]);
738  fHists[102] = new TH3D("fCumPhiEtaZNA",";#Delta#phi;#Delta#eta",32,-TMath::Pi()/2,3*TMath::Pi()/2,60,-3,3,10,0,100);
739  fOutput->Add(fHists[102]);
740  fHists[103] = new TH3D("fCumPtEtaCl1",";p_{T} (GeV/c);#eta",100,0,25,20,-2,2,10,0,100);
741  fOutput->Add(fHists[103]);
742  fHists[104] = new TH3D("fCumPtEtaV0A",";p_{T} (GeV/c);#eta",100,0,25,20,-2,2,10,0,100);
743  fOutput->Add(fHists[104]);
744  fHists[105] = new TH3D("fCumPtEtaZNA",";p_{T} (GeV/c);#eta",100,0,25,20,-2,2,10,0,100);
745  fOutput->Add(fHists[105]);
746  fHists[106] = new TH1D("fCumCL1MC",";#tracks",fCumMbins,0,fCumMbins);
747  fOutput->Add(fHists[106]);
748  fHists[107] = new TH1D("fCumV0AMC",";#tracks",fCumMbins,0,fCumMbins);
749  fOutput->Add(fHists[107]);
750  fHists[108] = new TH1D("fCumV0CMC",";#tracks",fCumMbins,0,fCumMbins);
751  fOutput->Add(fHists[108]);
752  fHists[109] = new TH1D("fCumCl1Cent",";percentile",10,0,100);
753  fOutput->Add(fHists[109]);
754  fHists[110] = new TH1D("fCumV0ACent",";percentile",10,0,100);
755  fOutput->Add(fHists[110]);
756  fHists[111] = new TH1D("fCumZNACent",";percentile",10,0,100);
757  fOutput->Add(fHists[111]);
758  fHists[112] = new TH1D("fCumMall",";mult",fCumMbins,0,fCumMbins);
759  fOutput->Add(fHists[112]);
760  fHists[113] = new TH1D("fCumM",";mult",fCumMbins,0,fCumMbins);
761  fOutput->Add(fHists[113]);
762  fHists[114] = new TProfile("fCumQC2",";qc2",fCumMbins,0,fCumMbins);
763  fOutput->Add(fHists[114]);
764  fHists[115] = new TProfile("fCumQC4",";qc4",fCumMbins,0,fCumMbins);
765  fOutput->Add(fHists[115]);
766  fHists[116] = new TProfile("fCumQC6",";qc6",fCumMbins,0,fCumMbins);
767  fOutput->Add(fHists[116]);
768  Int_t v0bins=1000;
769  if (fCumMbins>1000)
770  v0bins=25000;
771  fHists[117] = new TH2D("fCumV0ACentVsM",";v0a;M",v0bins,0,v0bins,fCumMbins,0,fCumMbins);
772  fOutput->Add(fHists[117]);
773  fHists[118] = new TH2D("fCumV0CCentVsM",";v0c;M",v0bins,0,v0bins,fCumMbins,0,fCumMbins);
774  fOutput->Add(fHists[118]);
775  fHists[119] = new TH2D("fCumV0MCentVsM",";v0m;M",v0bins,0,v0bins,fCumMbins,0,fCumMbins);
776  fOutput->Add(fHists[119]);
777  fHists[120] = new TProfile("fCum3QC2",";qc2",fCumMbins,0,fCumMbins);
778  fOutput->Add(fHists[120]);
779  fHists[121] = new TProfile("fCum3QC4",";qc4",fCumMbins,0,fCumMbins);
780  fOutput->Add(fHists[121]);
781  fHists[122] = new TProfile("fEtaGapV2",";M",fCumMbins,0,fCumMbins);
782  fOutput->Add(fHists[122]);
783  fHists[123] = new TProfile("fEtaGapV3",";M",fCumMbins,0,fCumMbins);
784  fOutput->Add(fHists[123]);
785  fHists[124] = new TProfile("fEtaGapV214",";M",fCumMbins,0,fCumMbins);
786  fOutput->Add(fHists[124]);
787  fHists[125] = new TProfile("fEtaGapV314",";M",fCumMbins,0,fCumMbins);
788  fOutput->Add(fHists[125]);
789  fHists[126] = new TProfile("fEtaGapV218",";M",fCumMbins,0,fCumMbins);
790  fOutput->Add(fHists[126]);
791  fHists[127] = new TProfile("fEtaGapV318",";M",fCumMbins,0,fCumMbins);
792  fOutput->Add(fHists[127]);
793  fHists[128] = new TH3D("fDPhiDEtaTracks",";#Delta#phi;#Delta#eta;M",64,-0.5*TMath::Pi(),1.5*TMath::Pi(),60,-3,3,
794  fCumMbins/10,0,fCumMbins);
795  fHists[128]->Sumw2();
796  fOutput->Add(fHists[128]);
797  fHists[129] = new TH1D("fDPhiDEtaTrigs",";M",fCumMbins/10,0,fCumMbins);
798  fHists[129]->Sumw2();
799  fOutput->Add(fHists[129]);
800 
801  fNtupCum = new TTree("NtupCum", "Ntuple for cumulant analysis");
802  if (1) {
803  fNtupCum->SetDirectory(0);
804  } else {
805  TFile *f = OpenFile(1);
806  if (f) {
807  f->SetCompressionLevel(2);
808  fNtupCum->SetDirectory(f);
809  fNtupCum->SetAutoFlush(-4*1024*1024);
810  fNtupCum->SetAutoSave(0);
811  }
812  }
814  fNtupCum->Branch("cumulants", &fNtupCumInfo, 32*1024, 99);
816  fNtupCum->Branch("zdc", &fNtupZdcInfo, 32*1024, 99);
817  if (fDoCumNtuple)
818  fOutput->Add(fNtupCum);
819  }
820 
821  PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
822 }
Base task in the EMCAL framework.
ClassImp(AliAnalysisTaskCLQA) AliAnalysisTaskCLQA
Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rangeMin=-TMath::Pi()/2, Double_t rangeMax=3 *TMath::Pi()/2) const
pointers to histograms
AliNtupZdcInfo * fNtupZdcInfo
object holding cumulant results
const Double_t etamin
void SetCumParams(Double_t Mmin, Double_t ptmin, Double_t ptmax, Double_t etamin, Double_t etamax)
const Double_t ptmax
const Double_t ptmin
Byte_t GetTrackType() const
Definition: AliPicoTrack.h:42
virtual Bool_t RetrieveEventObjects()
AliVTrack * GetTrack() const
Definition: AliPicoTrack.h:59
Bool_t isMC
AliEmcalList * fOutput
!output list
TClonesArray * fTracks
!tracks
TH1 * fHists[1000]
object holding zdc info
const Double_t etamax
const Double_t pi
AliNtupCumInfo * fNtupCumInfo
ntuple for cumulant analysis
void RunCumulants(Double_t Mmin, Double_t ptmin, Double_t ptmax, Double_t etamin, Double_t etamax)