AliPhysics  7f2a7c4 (7f2a7c4)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Cumulants.cxx
Go to the documentation of this file.
1 #include <Riostream.h>
2 #include <TComplex.h>
3 #include <TF1.h>
4 #include <TH1D.h>
5 #include <TLorentzVector.h>
6 #include <TMath.h>
7 #include <TParticlePDG.h>
8 #include <TProfile.h>
9 #include <TRandom.h>
10 #include <TRandom3.h>
11 #include <TStopwatch.h>
12 #include "Cumulants.h"
13 
14 CPart::CPart(const TParticle &p) : fPt(p.Pt()), fEta(p.Eta()), fPhi(TVector2::Phi_0_2pi(p.Phi()))
15 {
16  TParticlePDG *pdg = p.GetPDG(1);
17  fCharge = pdg->Charge();
18 }
19 
20 Cumulants::Cumulants(const char *name, Int_t mbins, Int_t minM) :
21  TNamed(name,""), fCumMBins(mbins), fMinM(minM), fEtaMin(-1), fEtaMax(1), fPtMin(0.3), fPtMax(3), fDoEtaGap(0),
22  fDoCharge(1), fDoQC(0), fDoQC44(0), fMaxNL4(50), fEGminNL4(0), fDoQC43(0), fMaxNL3(100), fEGminNL3(0),
23  fDoQC4withEG(0), fDoDebug(0), fDoPrint(0)
24 {
25  fList = new TList;
26  fList->SetName(name);
27  fParts.resize(9999);
28  fHists[0] = new TH1D("fCumMnocut",";mult",fCumMBins,0,fCumMBins);
29  fList->Add(fHists[0]);
30  fHists[1] = new TH1D("fCumMall",";mult",fCumMBins,0,fCumMBins);
31  fList->Add(fHists[1]);
32  fHists[2] = new TH1D("fCumMin",";mult",fCumMBins,0,fCumMBins);
33  fList->Add(fHists[2]);
34  fHists[3] = new TH1D("fCumM",";mult",fCumMBins,0,fCumMBins);
35  fList->Add(fHists[3]);
36  Info("Cumulants","Standalone version 1.0 started");
37 }
38 
40 {
41  // Add (and enable) QC4 measurements with one eta gap (for 1-3 and 2-4)
42  fDoQC4withEG = 1;
43 
44  Double_t eg=TMath::Abs(etagap);
45  fEGQCCuts.push_back(eg);
46  Int_t n = fEGQCCuts.size()-1;
47  Int_t hindex=400+4*n;
48  fHists[hindex] = new TH1D(Form("fCumMnegeta%d",n),Form(";mult (#eta<-%.1f)",eg),fCumMBins,0,fCumMBins);
49  fList->Add(fHists[hindex++]);
50  fHists[hindex] = new TH1D(Form("fCumMposeta%d",n),Form(";mult (#eta>+%.1f)",eg),fCumMBins,0,fCumMBins);
51  fList->Add(fHists[hindex++]);
52  fHists[hindex] = new TProfile(Form("fCumQC2wEG%d",n),Form(";M;qc2 (|#eta|>%.1f)",eg),fCumMBins,0,fCumMBins);
53  fList->Add(fHists[hindex++]);
54  fHists[hindex] = new TProfile(Form("fCumQC4wEG%d",n),Form(";M;qc4 (|#eta|>%.1f)",eg),fCumMBins,0,fCumMBins);
55  fList->Add(fHists[hindex++]);
56 }
57 
59 {
60  // Enable eta gap measurements
61  Double_t etacuts[] = {0.1,0.2,0.3,0.4,0.5,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.5};
62  const Int_t N = 14;
63  fEGCuts.resize(N);
64  fEGC2.resize(N);
65  fEGC3.resize(N);
66  fEGS2.resize(N);
67  fEGS3.resize(N);
68  fEGCounts.resize(N);
69  Int_t hindex=100;
70  for (Int_t i=0;i<N;++i) {
71  fEGCuts[i]=etacuts[i];
72  fHists[hindex] = new TProfile(Form("fEtaGapC2%02d",Int_t(etacuts[i]*10)),Form(";M;#LTcos(2#phi_{1}-2#phi_{2})#GT (|#eta|>%.1f)",etacuts[i]),fCumMBins,0,fCumMBins);
73  fList->Add(fHists[hindex]);
74  ++hindex;
75  fHists[hindex] = new TProfile(Form("fEtaGapC3%02d",Int_t(etacuts[i]*10)),Form(";M;#LTcos(3#phi_{1}-3#phi_{2})#GT (|#eta|>%.1f)",etacuts[i]),fCumMBins,0,fCumMBins);
76  fList->Add(fHists[hindex]);
77  ++hindex;
78  fHists[hindex] = new TProfile(Form("fEtaGapS2%02d",Int_t(etacuts[i]*10)),Form(";M;#LTsin(2#phi_{1}-2#phi_{2})#GT (|#eta|>%.1f)",etacuts[i]),fCumMBins,0,fCumMBins);
79  fList->Add(fHists[hindex]);
80  ++hindex;
81  fHists[hindex] = new TProfile(Form("fEtaGapS3%02d",Int_t(etacuts[i]*10)),Form(";M;#LTsin(3#phi_{1}-3#phi_{2})#GT (|#eta|>%.1f)",etacuts[i]),fCumMBins,0,fCumMBins);
82  fList->Add(fHists[hindex]);
83  ++hindex;
84  }
85  fDoEtaGap = 1;
86 }
87 
89 {
90  // Enable QC measurements
91  fDoQC = 1;
92  Int_t hindex=200;
93  fHists[hindex] = new TProfile("fCumQ2r",";M;q2r",fCumMBins,0,fCumMBins);
94  fList->Add(fHists[hindex++]);
95  fHists[hindex] = new TProfile("fCumQ2i",";M;q2i",fCumMBins,0,fCumMBins);
96  fList->Add(fHists[hindex++]);
97  fHists[hindex] = new TProfile("fCumQ3r",";M;q3r",fCumMBins,0,fCumMBins);
98  fList->Add(fHists[hindex++]);
99  fHists[hindex] = new TProfile("fCumQ3i",";M;q3i",fCumMBins,0,fCumMBins);
100  fList->Add(fHists[hindex++]);
101  fHists[hindex] = new TProfile("fCumQ4r",";M;q4r",fCumMBins,0,fCumMBins);
102  fList->Add(fHists[hindex++]);
103  fHists[hindex] = new TProfile("fCumQ4i",";M;q4i",fCumMBins,0,fCumMBins);
104  fList->Add(fHists[hindex++]);
105  fHists[hindex] = new TProfile("fCumQ6r",";M;q6r",fCumMBins,0,fCumMBins);
106  fList->Add(fHists[hindex++]);
107  fHists[hindex] = new TProfile("fCumQ6i",";M;q6i",fCumMBins,0,fCumMBins);
108  fList->Add(fHists[hindex++]);
109  hindex=230;
110  fHists[hindex] = new TProfile("fCumQC2",";M;qc2",fCumMBins,0,fCumMBins);
111  fList->Add(fHists[hindex++]);
112  fHists[hindex] = new TProfile("fCumQC4",";M;qc4",fCumMBins,0,fCumMBins);
113  fList->Add(fHists[hindex++]);
114  fHists[hindex] = new TProfile("fCumQC6",";M;qc6",fCumMBins,0,fCumMBins);
115  fList->Add(fHists[hindex++]);
116  fHists[hindex] = new TProfile("fCum3QC2",";M;qc2",fCumMBins,0,fCumMBins);
117  fList->Add(fHists[hindex++]);
118  fHists[hindex] = new TProfile("fCum3QC4",";M;qc4",fCumMBins,0,fCumMBins);
119  if (fDoCharge) {
120  hindex=260;
121  fHists[hindex] = new TH1D("fCumMpp",";mult (pos)",fCumMBins,0,fCumMBins);
122  fList->Add(fHists[hindex++]);
123  fHists[hindex] = new TH1D("fCumMnn",";mult (neg)",fCumMBins,0,fCumMBins);
124  fList->Add(fHists[hindex++]);
125  fHists[hindex] = new TProfile("fCumQC2pp",";M;qc2 (pos)",fCumMBins,0,fCumMBins);
126  fList->Add(fHists[hindex++]);
127  fHists[hindex] = new TProfile("fCum3QC2pp",";M;qc2 (pos)",fCumMBins,0,fCumMBins);
128  fList->Add(fHists[hindex++]);
129  fHists[hindex] = new TProfile("fCumQC4pp",";M;qc4 (pos)",fCumMBins,0,fCumMBins);
130  fList->Add(fHists[hindex++]);
131  fHists[hindex] = new TProfile("fCum3QC4pp",";M;qc4 (pos)",fCumMBins,0,fCumMBins);
132  fList->Add(fHists[hindex++]);
133  fHists[hindex] = new TProfile("fCumQC2nn",";M;qc2 (neg)",fCumMBins,0,fCumMBins);
134  fList->Add(fHists[hindex++]);
135  fHists[hindex] = new TProfile("fCumQC4nn",";M;qc4 (neg)",fCumMBins,0,fCumMBins);
136  fList->Add(fHists[hindex++]);
137  fHists[hindex] = new TProfile("fCum3QC2nn",";M;qc2 (neg)",fCumMBins,0,fCumMBins);
138  fList->Add(fHists[hindex++]);
139  fHists[hindex] = new TProfile("fCum3QC4nn",";M;qc4 (neg)",fCumMBins,0,fCumMBins);
140  fList->Add(fHists[hindex++]);
141  hindex=280;
142  fHists[hindex] = new TProfile("fCumQC2ss",";M;qc2 (ss)",fCumMBins,0,fCumMBins);
143  fList->Add(fHists[hindex++]);
144  fHists[hindex] = new TProfile("fCumQC4ss",";M;qc4 (ss)",fCumMBins,0,fCumMBins);
145  fList->Add(fHists[hindex++]);
146  fHists[hindex] = new TProfile("fCum3QC2ss",";M;qc2 (ss)",fCumMBins,0,fCumMBins);
147  fList->Add(fHists[hindex++]);
148  fHists[hindex] = new TProfile("fCum3QC4ss",";M;qc4 (ss)",fCumMBins,0,fCumMBins);
149  fList->Add(fHists[hindex++]);
150  }
151 }
152 
154 {
155  // Enable QC measurements with 4NL
156  fDoQC44 = 1;
157  fMaxNL4 = mn;
158  fEGminNL4 = etamin;
159  Int_t hindex=300;
160  for (Int_t i=0;i<6;++i) {
161  Double_t etagap=(Double_t)i/10;
162  fHists[hindex] = new TProfile(Form("fQC4NL42%d",i),Form(";M;QC4 (|#eta|>%.1f)",etagap),fCumMBins,0,fCumMBins);
163  fList->Add(fHists[hindex]);
164  ++hindex;
165  }
166 }
167 
169 {
170  // Enable QC measurements with 3NL
171  fDoQC43 = 1;
172  fMaxNL3 = mn;
173  fEGminNL3 = etamin;
174  Int_t hindex=350;
175  for (Int_t i=0;i<6;++i) {
176  Double_t etagap=(Double_t)i/10;
177  fHists[hindex] = new TProfile(Form("fQC3NL42%d",i),Form(";M;QC4 (|#eta|>%.1f)",etagap),fCumMBins,0,fCumMBins);
178  fList->Add(fHists[hindex]);
179  ++hindex;
180  }
181 }
182 
184 {
185  // Enable QC4 measurements with one eta gap (for 1-3 and 2-4)
186  for (Double_t eg=0.05;eg<=0.5;eg+=0.05)
187  AddQC4withEG(eg);
188 }
189 
191 {
192  // Run all enabled methods
193  if (fM<fMinM)
194  return;
195  fHists[3]->Fill(fM);
196 
197  TStopwatch w;
198  if (fDoEtaGap) {
199  w.Start();
200  RunEG();
201  w.Stop();
202  }
203  if (fDoQC) {
204  w.Start();
205  RunQC();
206  w.Stop();
207  }
208  if (fDoQC44) {
209  w.Start();
210  RunQC4with4NL();
211  w.Stop();
212  }
213  if (fDoQC43) {
214  w.Start();
215  RunQC4with3NL();
216  w.Stop();
217  }
218  if (fDoQC4withEG) {
219  w.Start();
220  RunQC4withEG();
221  w.Stop();
222  }
223 }
224 
226 {
227  // Run eta gap two particle method
228  const UInt_t N=fEGCuts.size();
229  fEGC2.assign(N,0);
230  fEGC3.assign(N,0);
231  fEGS2.assign(N,0);
232  fEGS3.assign(N,0);
233  fEGCounts.assign(N,0);
234  for (Int_t i=0; i<fM; ++i) {
235  CPart &p1 = fParts.at(i);
236  for (Int_t j=i+1; j<fM; ++j) {
237  CPart &p2 = fParts.at(j);
238  Double_t dphi=p1.Phi()-p2.Phi();
239  Double_t deta=TMath::Abs(p1.Eta()-p2.Eta());
240  Double_t c2 = TMath::Cos(2*dphi);
241  Double_t s2 = TMath::Sin(2*dphi);
242  Double_t c3 = TMath::Cos(3*dphi);
243  Double_t s3 = TMath::Sin(3*dphi);
244  for (UInt_t k=0;k<N;++k) {
245  if (deta>fEGCuts.at(k)) {
246  fEGC2[k] += c2;
247  fEGC3[k] += c3;
248  fEGS2[k] += s2;
249  fEGS3[k] += s3;
250  fEGCounts[k] += 1;
251  } else
252  break;
253  }
254  }
255  }
256  Int_t hind=100;
257  for (UInt_t i=0;i<N;++i) {
258  if (fEGCounts[i]<=0)
259  continue;
260  fEGC2[i]/=fEGCounts[i];
261  fHists[hind++]->Fill(fM,fEGC2[i]);
262  fEGC3[i]/=fEGCounts[i];
263  fHists[hind++]->Fill(fM,fEGC3[i]);
264  fEGS2[i]/=fEGCounts[i];
265  fHists[hind++]->Fill(fM,fEGS2[i]);
266  fEGS3[i]/=fEGCounts[i];
267  fHists[hind++]->Fill(fM,fEGS3[i]);
268  }
269 }
270 
272 {
273  // Run QC method
274  fQC[2]=0;fQC[3]=0;fQC[4]=0;fQC[5]=0;fQC[6]=0;
275  for (Int_t i=0; i<fM; ++i) {
276  CPart &p = fParts.at(i);
277  Double_t phi=p.Phi();
278  for (Int_t k=2; k<7; ++k) {
279  fQC[k] += TComplex(TMath::Cos(k*phi),TMath::Sin(k*phi));
280  }
281  }
282 
283  Double_t Q22 = fQC[2].Rho2();
284  Double_t Q32 = fQC[3].Rho2();
285  Double_t Q42 = fQC[4].Rho2();
286  Double_t Q52 = fQC[5].Rho2();
287  Double_t Q62 = fQC[6].Rho2();
288  Double_t Q32re = (fQC[6]*TComplex::Power(TComplex::Conjugate(fQC[3]),2)).Re();
289  Double_t Q42re = TComplex(fQC[4]*TComplex::Power(TComplex::Conjugate(fQC[2]),2)).Re();
290  Double_t Q6are = TComplex(fQC[4]*fQC[2]*TComplex::Power(TComplex::Conjugate(fQC[2]),3)).Re();
291  Double_t Q6bre = TComplex(fQC[6]*TComplex::Power(TComplex::Conjugate(fQC[2]),3)).Re();
292  Double_t Q6cre = TComplex(fQC[6]*TComplex::Conjugate(fQC[4])*TComplex::Conjugate(fQC[2])).Re();
293 
294  Int_t M = fM;
295  Int_t M2 = M*(M-1);
296  Int_t M4 = M*(M-1)*(M-2)*(M-3);
297  if (M>1) {
298  Int_t hindex=200;
299  fHists[hindex++]->Fill(M,fQC[2].Re()/M);
300  fHists[hindex++]->Fill(M,fQC[2].Im()/M);
301  fHists[hindex++]->Fill(M,fQC[3].Re()/M);
302  fHists[hindex++]->Fill(M,fQC[3].Im()/M);
303  fHists[hindex++]->Fill(M,fQC[4].Re()/M);
304  fHists[hindex++]->Fill(M,fQC[4].Im()/M);
305  fHists[hindex++]->Fill(M,fQC[6].Re()/M);
306  fHists[hindex++]->Fill(M,fQC[6].Im()/M);
307  hindex=230;
308  fHists[hindex++]->Fill(M,(Q22-M)/M2);
309  if (M>3) {
310  Double_t qc42tmp = (Q22*Q22+Q42-2*Q42re-4*(M-2)*Q22+2*M*(M-3));
311  fHists[hindex]->Fill(M,qc42tmp/M4);
312  }
313  hindex++;
314  if (M>5) {
315  Double_t qc6tmp = Q22*Q22*Q22 + 9*Q42*Q22 - 6*Q6are
316  + 4*Q6bre - 12*Q6cre
317  + 18*(M-4)*Q42re + 4*Q62
318  - 9*(M-4)*Q22*Q22 - 9*(M-4)*Q42
319  + 18*(M-2)*(M-5)*Q22
320  - 6*M*(M-4)*(M-5);
321  fHists[hindex]->Fill(M,qc6tmp/M/(M-1)/(M-2)/(M-3)/(M-4)/(M-5));
322  }
323  hindex++;
324  fHists[hindex++]->Fill(M,(Q32-M)/M2);
325  if (M>3) {
326  Double_t qc43tmp = (Q32*Q32+Q62-2*Q32re-4*(M-2)*Q32+2*M*(M-3));
327  fHists[hindex]->Fill(M,qc43tmp/M4);
328  }
329  }
330 
331  if (!fDoCharge)
332  return;
333 
334  Int_t Mpp=0,Mnn=0;
335  TComplex qcpp[7],qcnn[7];
336  qcpp[2]=0;qcpp[3]=0;qcpp[4]=0,qcpp[5]=0,qcpp[6]=0;
337  qcnn[2]=0;qcnn[3]=0;qcnn[4]=0,qcnn[5]=0,qcnn[6]=0;
338  for (Int_t i=0; i<fM; ++i) {
339  CPart &p = fParts.at(i);
340  Double_t phi=p.Phi();
341  Short_t c=p.Charge();
342  if (c==0)
343  continue;
344  if (c>0) {
345  ++Mpp;
346  for (Int_t k=2; k<5; ++k)
347  qcpp[k] += TComplex(TMath::Cos(k*phi),TMath::Sin(k*phi));
348  } else {
349  ++Mnn;
350  for (Int_t k=2; k<5; ++k)
351  qcnn[k] += TComplex(TMath::Cos(k*phi),TMath::Sin(k*phi));
352  }
353  }
354 
355  Double_t Q22pp = qcpp[2].Rho2();
356  Double_t Q32pp = qcpp[3].Rho2();
357  Double_t Q42pp = qcpp[4].Rho2();
358  Double_t Q62pp = qcpp[6].Rho2();
359  Double_t Q32ppre = (qcpp[6]*TComplex::Power(TComplex::Conjugate(qcpp[3]),2)).Re();
360  Double_t Q42ppre = (qcpp[4]*TComplex::Power(TComplex::Conjugate(qcpp[2]),2)).Re();
361  Double_t Q22nn = qcnn[2].Rho2();
362  Double_t Q32nn = qcnn[3].Rho2();
363  Double_t Q42nn = qcnn[4].Rho2();
364  Double_t Q62nn = qcnn[6].Rho2();
365  Double_t Q32nnre = (qcnn[6]*TComplex::Power(TComplex::Conjugate(qcnn[3]),2)).Re();
366  Double_t Q42nnre = (qcnn[4]*TComplex::Power(TComplex::Conjugate(qcnn[2]),2)).Re();
367 
368  Int_t hindex=260;
369  const Int_t sind=280;
370  fHists[260]->Fill(Mpp);
371  fHists[261]->Fill(Mnn);
372  Int_t Mpp2 = Mpp*(Mpp-1);
373  Int_t Mpp4 = Mpp*(Mpp-1)*(Mpp-2)*(Mpp-3);
374  Int_t Mnn2 = Mnn*(Mnn-1);
375  Int_t Mnn4 = Mnn*(Mnn-1)*(Mnn-2)*(Mnn-3);
376  if (Mpp>1) {
377  fHists[262]->Fill(Mpp,(Q22pp-Mpp)/Mpp2);
378  fHists[280]->Fill(M,(Q22pp-Mpp)/Mpp2);
379  fHists[263]->Fill(Mpp,(Q32pp-Mpp)/Mpp2);
380  fHists[282]->Fill(M,(Q32pp-Mpp)/Mpp2);
381  if (Mpp>3) {
382  Double_t qc4tmp = (Q22pp*Q22pp+Q42pp-2*Q42ppre-4*(Mpp-2)*Q22pp+2*Mpp*(Mpp-3));
383  fHists[264]->Fill(Mpp,qc4tmp/Mpp4);
384  fHists[281]->Fill(M,qc4tmp/Mpp4);
385  Double_t qc43tmp = (Q32pp*Q32pp+Q62pp-2*Q32ppre-4*(Mpp-2)*Q32pp+2*Mpp*(Mpp-3));
386  fHists[265]->Fill(Mpp,qc43tmp/Mpp4);
387  fHists[283]->Fill(M,qc43tmp/Mpp4);
388  }
389  }
390  if (Mnn>1) {
391  fHists[266]->Fill(Mnn,(Q22nn-Mnn)/Mnn2);
392  fHists[280]->Fill(M,(Q22nn-Mnn)/Mnn2);
393  fHists[267]->Fill(Mnn,(Q32nn-Mnn)/Mnn2);
394  fHists[282]->Fill(M,(Q32nn-Mnn)/Mnn2);
395  if (Mnn>3) {
396  Double_t qc4tmp = (Q22nn*Q22nn+Q42nn-2*Q42nnre-4*(Mnn-2)*Q22nn+2*Mnn*(Mnn-3));
397  fHists[268]->Fill(Mnn,qc4tmp/Mnn4);
398  fHists[281]->Fill(M,qc4tmp/Mnn4);
399  Double_t qc43tmp = (Q32nn*Q32nn+Q62nn-2*Q32nnre-4*(Mnn-2)*Q32nn+2*Mnn*(Mnn-3));
400  fHists[269]->Fill(Mnn,qc43tmp/Mnn4);
401  fHists[283]->Fill(M,qc43tmp/Mnn4);
402  }
403  }
404 }
405 
407 {
408  // Run QC method with 4 nested loops
409  if (fM>fMaxNL4)
410  return;
411 
412  Double_t er[6] = {0};
413  Double_t nc[6] = {0};
414 
415  for (Int_t i1=0; i1<fM; ++i1) {
416  CPart &p1 = fParts.at(i1);
417  Double_t eta1=p1.Eta();
418  Double_t phi1=p1.Phi();
419  for (Int_t i2=0; i2<fM; ++i2) {
420  if (i2==i1)
421  continue;
422  CPart &p2 = fParts.at(i2);
423  Double_t eta2=p2.Eta();
424  Double_t eta12=TMath::Abs(eta1-eta2);
425  if (eta12<fEGminNL4)
426  continue;
427  Double_t phi2=p2.Phi();
428  for (Int_t i3=0; i3<fM; ++i3) {
429  if (i3==i2)
430  continue;
431  if (i3==i1)
432  continue;
433  CPart &p3 = fParts.at(i3);
434  Double_t eta3=p3.Eta();
435  Double_t eta13=TMath::Abs(eta1-eta3);
436  if (eta13<fEGminNL4)
437  continue;
438  Double_t eta23=TMath::Abs(eta2-eta3);
439  if (eta23<fEGminNL4)
440  continue;
441  Double_t phi3=p3.Phi();
442  for (Int_t i4=0; i4<fM; ++i4) {
443  if (i4==i3)
444  continue;
445  if (i4==i2)
446  continue;
447  if (i4==i1)
448  continue;
449  CPart &p4 = fParts.at(i4);
450  Double_t eta4=p4.Eta();
451  Double_t phi4=p4.Phi();
452  Double_t eta14=TMath::Abs(eta1-eta4);
453  Double_t eta24=TMath::Abs(eta2-eta4);
454  Double_t eta34=TMath::Abs(eta3-eta4);
455  if (eta14<fEGminNL4)
456  continue;
457  if (eta24<fEGminNL4)
458  continue;
459  if (eta23<fEGminNL4)
460  continue;
461  Double_t arg=(phi1+phi2-phi3-phi4);
462  Double_t val = TMath::Cos(2*arg);;
463  for (Int_t eg=0;eg<6;++eg) {
464  Double_t etagap=Double_t(eg)/10;
465  if ((eta12>etagap)&&(eta13>etagap)&&(eta14>etagap)&&
466  (eta23>etagap)&&(eta24>etagap)&&(eta34>etagap)) {
467  er[eg] += val;
468  nc[eg] += 1;
469  }
470  }
471  }
472  }
473  }
474  }
475 
476  Int_t hindex=300;
477  for (Int_t eg=0;eg<6;++eg) {
478  Int_t n=nc[eg];
479  if (n<=0)
480  continue;
481  Double_t val = er[eg]/n;
482  fHists[hindex+eg]->Fill(fM,val);
483  }
484 }
485 
487 {
488  // Run QC method with 3 nested loops (from Dhevan's thesis with 7.19 corrected)
489  if (fM>fMaxNL3)
490  return;
491 
492  TComplex nq2[6],nq3[6],nq4[6];
493  Int_t np[6]={0}, ns[6]={0};
494  for (Int_t i1=0; i1<fM; ++i1) {
495  CPart &p1 = fParts.at(i1);
496  Double_t eta1=p1.Eta();
497  Double_t phi1=p1.Phi();
498  for (Int_t i2=0; i2<fM; ++i2) {
499  if (i1==i2)
500  continue;
501  CPart &p2 = fParts.at(i2);
502  Double_t eta2=p2.Eta();
503  Double_t eta12=TMath::Abs(eta1-eta2);
504  if (eta12<fEGminNL3)
505  continue;
506  Double_t phi2=p2.Phi();
507  Double_t phi12=phi1-phi2;
508  TComplex v2(TMath::Cos(2*phi12),TMath::Sin(2*phi12));
509  TComplex v4(1+TMath::Cos(4*phi12),TMath::Sin(4*phi12));
510  for (Int_t eg=0;eg<6;++eg) {
511  Double_t etagap=Double_t(eg)/10;
512  if (eta12>etagap) {
513  nq2[eg] += v2;
514  nq4[eg] += v4;
515  np[eg] += 1;
516  ns[eg] += 2;
517  }
518  }
519  for (Int_t i3=0; i3<fM; ++i3) {
520  if (i1==i3)
521  continue;
522  if (i2==i3)
523  continue;
524  CPart &p3 = fParts.at(i3);
525  Double_t eta3=p3.Eta();
526  Double_t phi3=p3.Phi();
527  Double_t eta13=TMath::Abs(eta1-eta3);
528  if (eta13<fEGminNL3)
529  continue;
530  Double_t eta23=TMath::Abs(eta2-eta3);
531  if (eta23<fEGminNL3)
532  continue;
533  Double_t phi13=phi1-phi3;
534  Double_t phi23=phi2-phi3;
535  Double_t t1=2*TMath::Cos(2*(phi12+phi13))+2*TMath::Cos(2*(phi12-phi13));
536  Double_t t2=2*TMath::Sin(2*(phi12+phi13))+2*TMath::Sin(2*(phi12-phi13));
537  TComplex val(t1,t2);
538  for (Int_t eg=0;eg<6;++eg) {
539  Double_t etagap=Double_t(eg)/10;
540  if ((eta13>etagap)&&(eta23>etagap)) {
541  nq3[eg] += val;
542  ns[eg] += 4;
543  }
544  }
545  }
546  }
547  }
548 
549  Int_t hindex=350;
550  for (Int_t eg=0;eg<6;++eg) {
551  Int_t n=np[eg]*np[eg]-ns[eg];
552  if (n<=0)
553  continue;
554  Double_t val=TComplex(nq2[eg]*nq2[eg]-nq4[eg]-nq3[eg]).Re()/n;
555  fHists[hindex+eg]->Fill(fM,val);
556  }
557 }
558 
559 void Cumulants::RunQC4withEG(Double_t etagap, Int_t &nn, Int_t &np, Double_t &val2, Double_t &val4)
560 {
561  // Run QC4 with one eta gap (between 1-3, and 2-4, from https://aliceinfo.cern.ch/Notes/node/554)
562 
563  const Double_t el=-TMath::Abs(etagap);
564  const Double_t eh=+TMath::Abs(etagap);
565  TComplex cp[5][5]={0};
566  TComplex cm[5][5]={0};
567  nn=0,np=0;
568  for (Int_t i=0; i<fM; ++i) {
569  CPart &p = fParts.at(i);
570  Double_t phi=p.Phi();
571  Double_t eta=p.Eta();
572  if (eta<el) {
573  ++nn;
574  for(Int_t iharm=0; iharm<5; ++iharm) {
575  for(Int_t ipow=0; ipow<5; ++ipow) {
576  cm[iharm][ipow] += TComplex(TMath::Power(TMath::Cos(iharm*phi),ipow),TMath::Power(TMath::Sin(iharm*phi),ipow));
577  }
578  }
579  } else if (eta>eh) {
580  ++np;
581  for(Int_t iharm=0; iharm<5; ++iharm) {
582  for(Int_t ipow=0; ipow<5; ++ipow) {
583  cp[iharm][ipow] += TComplex(TMath::Power(TMath::Cos(iharm*phi),ipow),TMath::Power(TMath::Sin(iharm*phi),ipow));
584  }
585  }
586  }
587  }
588 
589  if (nn>0&&np>0)
590  val2 = TComplex(cp[2][1]*TComplex::Conjugate(cm[2][1])).Re()/nn/np;
591 
592  TComplex Dn4Gap(cp[0][1]*cp[0][1]*cm[0][1]*cm[0][1] - cp[0][2]*cm[0][1]*cm[0][1] - cp[0][1]*cp[0][1]*cm[0][2] + cp[0][2]*cm[0][2]);
593  TComplex v24Gap(cp[2][1]*cp[2][1]*TComplex::Conjugate(cm[2][1])*TComplex::Conjugate(cm[2][1])
594  - cp[4][2]*TComplex::Conjugate(cm[2][1])*TComplex::Conjugate(cm[2][1])
595  - cp[2][1]*cp[2][1]*TComplex::Conjugate(cm[4][2]) + cp[4][2]*TComplex::Conjugate(cm[4][2]));
596  if (Dn4Gap.Re()>0)
597  val4 = v24Gap.Re()/Dn4Gap.Re();
598 }
599 
601 {
602  // Run QC4 with two eta gaps (between 1-3, and 2-4, from https://aliceinfo.cern.ch/Notes/node/554)
603 
604  Int_t nn=0,np=0;
605  Double_t val2=0,val4=0;
606 
607  Int_t n=fEGQCCuts.size();
608  for (Int_t i=0;i<n;++i) {
609  Double_t etagap = fEGQCCuts.at(i);
610  RunQC4withEG(etagap,nn,np,val2,val4);
611  Int_t hindex=400+4*i;
612  fHists[hindex++]->Fill(nn);
613  fHists[hindex++]->Fill(np);
614  if ((nn<2)||(np<2))
615  return;
616  fHists[hindex++]->Fill(fM,val2);
617  fHists[hindex++]->Fill(fM,val4);
618  }
619 }
620 
621 void Cumulants::SetTracks(TObjArray &trks, Bool_t doKinCuts)
622 {
623  // Set tracks from array, and do kinematic cuts if required.
624  fM = 0;
625  Int_t Mall = 0;
626  Int_t Mnoc = 0;
627 
628  const Int_t en = trks.GetEntries();
629  for (Int_t i=0;i<en;++i) {
630  CPart part;
631  AliVParticle *pav = dynamic_cast<AliVParticle*>(trks.At(i));
632  if (pav) {
633  part = CPart(*pav);
634  } else {
635  TParticle *par = dynamic_cast<TParticle*>(trks.At(i));
636  if (par) {
637  part = CPart(*par);
638  } else {
639  ::Error("Cumulants::SetTracks","Type of particle not known!");
640  continue;
641  }
642  }
643 
644  ++Mnoc;
645  if (doKinCuts) {
646  if (part.Eta()<fEtaMin)
647  continue;
648  if (part.Eta()>fEtaMax)
649  continue;
650  ++Mall;
651  if (part.Pt()<fPtMin)
652  continue;
653  if (part.Pt()>fPtMax)
654  continue;
655  }
656  fParts[fM]=part;
657  ++fM;
658  }
659  fHists[0]->Fill(Mnoc);
660  fHists[1]->Fill(Mall);
661  fHists[2]->Fill(fM);
662 }
663 
664 #if 0 /*from Alex*/
665 void v224pPb()
666 {
667  TProfile* hqqqqG = new TProfile("hqqqqG", ";multiplicity;Q_{1}Q_{2}Q_{3}^{*}Q_{4}^{*}/(M_{1}M_{2}M_{3}M_{4})", 30, 0, 150);
668 
669  TProfile* hqqG13 = new TProfile("hqqG13", ";multiplicity;Q_{1}Q_{2}^{*}/(M_{1}M_{2})", 30, 0, 150);
670  TProfile* hqqG24 = new TProfile("hqqG24", ";multiplicity;Q_{1}Q_{2}^{*}/(M_{1}M_{2})", 30, 0, 150);
671  TProfile* hqqG14 = new TProfile("hqqG14", ";multiplicity;Q_{1}Q_{2}^{*}/(M_{1}M_{2})", 30, 0, 150);
672  TProfile* hqqG23 = new TProfile("hqqG23", ";multiplicity;Q_{1}Q_{2}^{*}/(M_{1}M_{2})", 30, 0, 150);
673  TProfile* hqqG34 = new TProfile("hqqG34", ";multiplicity;Q_{1}Q_{2}^{*}/(M_{1}M_{2})", 30, 0, 150);
674 
675  for(Int_t n = 0; n < nEvents; n++) {
676  inTree->GetEntry(n);
677 
678  Int_t nTracks = trackArray->GetEntries();
679  Double_t QxRnd1g = 0, QyRnd1g = 0;
680  Double_t QxRnd2g = 0, QyRnd2g = 0;
681  Double_t QxRnd3g = 0, QyRnd3g = 0;
682  Double_t QxRnd4g = 0, QyRnd4g = 0;
683  Int_t mult1g = 0, mult2g = 0, mult3g = 0, mult4g = 0;
684 
685  for(Int_t i = 0; i < nTracks; i++) {
686  Track* trk = (Track*)trackArray->At(i);
687 
688  if ((trk->eta>-0.8) && (trk->eta<-0.55)){
689  QxRnd1g += TMath::Cos(2*trk->phi);
690  QyRnd1g += TMath::Sin(2*trk->phi);
691  mult1g++;
692  }
693 
694  if ((trk->eta>-0.35) && (trk->eta<-0.1)){
695  QxRnd2g += TMath::Cos(2*trk->phi);
696  QyRnd2g += TMath::Sin(2*trk->phi);
697  mult2g++;
698  }
699 
700  if ((trk->eta>0.1) && (trk->eta<0.35)){
701  QxRnd3g += TMath::Cos(2*trk->phi);
702  QyRnd3g += TMath::Sin(2*trk->phi);
703  mult3g++;
704  }
705 
706  if ((trk->eta>0.55) && (trk->eta<0.8)){
707  QxRnd4g += TMath::Cos(2*trk->phi);
708  QyRnd4g += TMath::Sin(2*trk->phi);
709  mult4g++;
710  }
711 
712  }
713 
714  if (mult1g != 0 && mult2g != 0 && mult3g != 0 && mult4g != 0){
715 
716  Double_t qqqqg = QxRnd1g*QxRnd2g*QxRnd3g*QxRnd4g - QxRnd1g*QxRnd2g*QyRnd3g*QyRnd4g + QxRnd1g*QyRnd2g*QxRnd3g*QyRnd4g + QxRnd1g*QyRnd2g*QyRnd3g*QxRnd4g + QyRnd1g*QxRnd2g*QxRnd3g*QyRnd4g + QyRnd1g*QxRnd2g*QyRnd3g*QxRnd4g - QyRnd1g*QyRnd2g*QxRnd3g*QxRnd4g + QyRnd1g*QyRnd2g*QyRnd3g*QyRnd4g;
717  hqqqqG->Fill(mult, qqqqg/Double_t(mult1g*mult2g*mult3g*mult4g));
718 
719  hqqG13->Fill(mult, (QxRnd1g*QxRnd3g + QyRnd1g*QyRnd3g)/Double_t(mult1g*mult3g));
720  hqqG24->Fill(mult, (QxRnd2g*QxRnd4g + QyRnd2g*QyRnd4g)/Double_t(mult2g*mult4g));
721  hqqG14->Fill(mult, (QxRnd1g*QxRnd4g + QyRnd1g*QyRnd4g)/Double_t(mult1g*mult4g));
722  hqqG23->Fill(mult, (QxRnd2g*QxRnd3g + QyRnd2g*QyRnd3g)/Double_t(mult2g*mult3g));
723  hqqG34->Fill(mult, (QxRnd3g*QxRnd4g + QyRnd3g*QyRnd4g)/Double_t(mult3g*mult4g));
724  }
725 }
726 #endif
Int_t pdg
Double_t fEGminNL4
Definition: Cumulants.h:51
Int_t fMaxNL3
Definition: Cumulants.h:53
Double_t fPtMin
Definition: Cumulants.h:44
std::vector< Double_t > fEGQCCuts
Definition: Cumulants.h:59
double Double_t
Definition: External.C:58
Double_t Eta() const
Definition: Cumulants.h:81
Bool_t fDoQC
Definition: Cumulants.h:48
Int_t fMaxNL4
Definition: Cumulants.h:50
void EnableQC4with4NL(Int_t mn=50, Double_t etamin=0.0)
Definition: Cumulants.cxx:153
void AddQC4withEG(Double_t etagap)
Definition: Cumulants.cxx:39
TList * fList
vector with particles
Definition: Cumulants.h:68
void EnableEG()
Definition: Cumulants.cxx:58
void RunQC()
Definition: Cumulants.cxx:271
TH1 * fHists[999]
list with histograms
Definition: Cumulants.h:69
std::vector< CPart > fParts
number of particles
Definition: Cumulants.h:67
Bool_t fDoQC4withEG
Definition: Cumulants.h:55
void RunAll()
Definition: Cumulants.cxx:190
Int_t fMinM
Definition: Cumulants.h:41
Int_t fCumMBins
Definition: Cumulants.h:40
Double_t Pt() const
Definition: Cumulants.h:80
TCanvas * c
Definition: TestFitELoss.C:172
std::vector< Double_t > fEGC3
eta gap c2
Definition: Cumulants.h:61
std::vector< Double_t > fEGCuts
Definition: Cumulants.h:58
Bool_t fDoCharge
Definition: Cumulants.h:47
std::vector< Double_t > fEGS2
eta gap c3
Definition: Cumulants.h:62
Double_t fPtMax
Definition: Cumulants.h:45
const Double_t etamin
void RunEG()
Definition: Cumulants.cxx:225
int Int_t
Definition: External.C:63
TComplex fQC[7]
eta counts per gap
Definition: Cumulants.h:65
unsigned int UInt_t
Definition: External.C:33
void EnableQC()
Definition: Cumulants.cxx:88
void SetTracks(TObjArray &trks, Bool_t doKinCuts=1)
Definition: Cumulants.cxx:621
std::vector< Double_t > fEGS3
eta gap s2
Definition: Cumulants.h:63
Definition: External.C:212
std::vector< Double_t > fEGC2
Definition: Cumulants.h:60
Double_t fEGminNL3
Definition: Cumulants.h:54
Double_t fEtaMin
Definition: Cumulants.h:42
Short_t fCharge
Definition: Cumulants.h:88
void RunQC4with3NL()
Definition: Cumulants.cxx:486
short Short_t
Definition: External.C:23
Bool_t fDoQC44
Definition: Cumulants.h:49
std::vector< Int_t > fEGCounts
eta gap s3
Definition: Cumulants.h:64
void EnableQC4withEG()
Definition: Cumulants.cxx:183
Float_t nEvents[nProd]
Bool_t fDoEtaGap
Definition: Cumulants.h:46
void RunQC4with4NL()
Definition: Cumulants.cxx:406
CPart()
Definition: Cumulants.h:75
bool Bool_t
Definition: External.C:53
Double_t Phi() const
Definition: Cumulants.h:82
void RunQC4withEG()
Definition: Cumulants.cxx:600
Int_t fM
QC values.
Definition: Cumulants.h:66
Cumulants(const char *name="cumulants", Int_t mbins=350, Int_t minM=10)
Definition: Cumulants.cxx:20
void EnableQC4with3NL(Int_t mn=100, Double_t etamin=0.0)
Definition: Cumulants.cxx:168
Short_t Charge() const
Definition: Cumulants.h:83
Double_t fEtaMax
Definition: Cumulants.h:43
Bool_t fDoQC43
Definition: Cumulants.h:52