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