Added analysis task for * Event by event PID fluctuation analysis
[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliAnalysisTaskPIDFluctuation.cxx
1 #include "AliAnalysisTaskPIDFluctuation.h"
2 #include "AliInputEventHandler.h"
3 #include "AliAnalysisManager.h"
4 #include "AliESDEvent.h"
5 #include "AliESDtrack.h"
6 #include "AliAODEvent.h"
7 #include "AliAODTrack.h"
8 #include "AliCentrality.h"
9 #include "AliPIDResponse.h"
10 #include "AliESDpid.h"
11 #include "AliAODpidUtil.h"
12 #include "TFile.h"
13 #include "TKey.h"
14 #include "TList.h"
15 #include "TH1F.h"
16 #include "TH2F.h"
17 #include "TH3F.h"
18 #include "TCanvas.h"
19 #include "AliESDtrackCuts.h"
20
21 /* 
22  * Event by event PID fluctuation analysis
23  * author: Roberto Preghenella (R+)
24  * email:  preghenella@bo.infn.it
25  *
26  */
27
28 ClassImp(AliAnalysisTaskPIDFluctuation)
29
30 //________________________________________________________________________
31
32 const Char_t *AliAnalysisTaskPIDFluctuation::fgkEventCounterName[AliAnalysisTaskPIDFluctuation::kNEventCounters] = {
33   "AllEvents",
34   "PhysicsSelection",
35   "PrimayVertex",
36   "PrimayVertexSPD",
37   "VertexAccepted",
38   "GoodCentrality",
39   "AcceptedEvents"
40 };
41
42 const Char_t *AliAnalysisTaskPIDFluctuation::fgkEventCounterTitle[AliAnalysisTaskPIDFluctuation::kNEventCounters] = {
43   "all events",
44   "physics selection",
45   "primary vertex",
46   "primary vertex SPD",
47   "vertex accepted",
48   "good centrality",
49   "accepted events"
50 };
51
52 const Char_t *AliAnalysisTaskPIDFluctuation::fgkSparseDataName[AliAnalysisTaskPIDFluctuation::kNSparseData] = {
53   "centV0M",
54   "centTRK",
55   "Nch",
56   "Nplus",
57   "Nminus",
58   "Npi",
59   "Npiplus",
60   "Npiminus",
61   "Nka",
62   "Nkaplus",
63   "Nkaminus",
64   "Npr",
65   "Nprplus",
66   "Nprminus"
67 };
68
69 const Char_t *AliAnalysisTaskPIDFluctuation::fgkSparseDataTitle[AliAnalysisTaskPIDFluctuation::kNSparseData] = {
70   "centrality percentile (V0M)",
71   "centrality percentile (TRK)",
72   "N_{charged}",
73   "N_{+}",
74   "N_{-}",
75   "N_{#pi}",
76   "N_{#pi^{+}}",
77   "N_{#pi^{-}}",
78   "N_{K}",
79   "N_{K^{+}}",
80   "N_{K^{-}}",
81   "N_{p+#bar{p}}",
82   "N_{p}",
83   "N_{#bar{p}}"
84 };
85
86 //________________________________________________________________________
87
88 AliAnalysisTaskPIDFluctuation::AliAnalysisTaskPIDFluctuation(const Char_t *name) : 
89   AliAnalysisTaskSE(name),
90   fPIDMethod(kTPCTOF),
91   fESDtrackCuts(NULL),
92   fAODfilterBit(AliAODTrack::kTrkGlobal),
93   fEtaMin(-0.8),
94   fEtaMax(0.8),
95   fPtMin(0.3),
96   fPtMax(1.5),
97   fPID(NULL),
98   fHistoList(NULL)
99 {
100   
101   /*
102    * default constructor
103    */
104
105   DefineOutput(1, TList::Class());
106 }
107
108 //________________________________________________________________________
109
110 AliAnalysisTaskPIDFluctuation::~AliAnalysisTaskPIDFluctuation()
111 {
112   
113   /*
114    * default destructor
115    */
116
117   if (fPID) delete fPID;
118   if (fHistoList) delete fHistoList;
119   
120 }
121
122 //________________________________________________________________________
123
124 void AliAnalysisTaskPIDFluctuation::UserCreateOutputObjects() 
125 {
126   
127   /*
128    * user create output objects
129    */
130
131   fHistoList = new TList();
132   fHistoList->SetOwner(kTRUE);
133
134   fHistoEventCounter = new TH1F("hHistoEventCounter", "", kNEventCounters, 0, kNEventCounters);
135   for (Int_t ievc = 0; ievc < kNEventCounters; ievc++)
136     fHistoEventCounter->GetXaxis()->SetBinLabel(ievc + 1, fgkEventCounterTitle[ievc]);
137   fHistoList->Add(fHistoEventCounter);
138
139   fHistoAcceptedTracks = new TH2F("hHistoAcceptedTracks", ";p_{T} (GeV/c)", 10, 0., 100., 50, 0., 5.);
140   fHistoList->Add(fHistoAcceptedTracks);
141   
142   fHistoTOFMatchedTracks = new TH2F("hHistoTOFMatchedTracks", ";p_{T} (GeV/c)", 10, 0., 100., 50, 0., 5.);
143   fHistoList->Add(fHistoTOFMatchedTracks);
144   
145   fHistoTPCdEdx = new TH3F("hHistoTPCdEdx", ";centrality percentile;p_{TPCin} (GeV/c);dE/dx (au.)", 10, 0., 100., 50, 0., 5., 500, 0., 500.);
146   fHistoList->Add(fHistoTPCdEdx);
147   
148   fHistoTPCdEdx_inclusive = new TH3F("hHistoTPCdEdx_inclusive", ";centrality percentile;p_{TPCin} (GeV/c);dE/dx (au.)", 10, 0., 100., 50, 0., 5., 500, 0., 500.);
149   fHistoList->Add(fHistoTPCdEdx_inclusive);
150   
151   fHistoTOFbeta = new TH3F(Form("hHistoTOFbeta"), ";centrality percentile;p (GeV/c);v/c", 10, 0., 100., 50, 0., 5., 500, 0.1, 1.1);
152   fHistoList->Add(fHistoTOFbeta);
153   
154   /* loop over species */
155   for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
156     
157     fHistoTPCdEdx_selected[ipart] = new TH3F(Form("hHistoTPCdEdx_selected_%s", AliPID::ParticleName(ipart)), ";centrality percentile;p_{TPCin} (GeV/c);dE/dx (au.)", 10, 0., 100., 50, 0., 5., 500, 0., 500.);
158     fHistoList->Add(fHistoTPCdEdx_selected[ipart]);
159     
160     fHistoTOFbeta_selected[ipart] = new TH3F(Form("hHistoTOFbeta_selected_%s", AliPID::ParticleName(ipart)), ";centrality percentile;p (GeV/c);v/c", 10, 0., 100., 50, 0., 5., 500, 0.1, 1.1);
161     fHistoList->Add(fHistoTOFbeta_selected[ipart]);
162     
163     fHistoNSigmaTPC[ipart] = new TH3F(Form("hHistoNSigmaTPC_%s", AliPID::ParticleName(ipart)), Form(";centrality percentile;p_{T} (GeV/c); N_{#sigma-TPC}^{%s}", AliPID::ParticleLatexName(ipart)), 10, 0., 100., 50, 0., 5., 200, -10., 10.);
164     fHistoList->Add(fHistoNSigmaTPC[ipart]);
165     
166     fHistoNSigmaTOF[ipart] = new TH3F(Form("hHistoNSigmaTOF_%s", AliPID::ParticleName(ipart)), Form(";centrality percentile;p_{T} (GeV/c); N_{#sigma-TOF}^{%s}", AliPID::ParticleLatexName(ipart)), 10, 0., 100., 50, 0., 5., 200, -10., 10.);
167     fHistoList->Add(fHistoNSigmaTOF[ipart]);
168     
169   } /* end of loop over species */
170
171   Int_t fgSparseDataBins[kNSparseData] = {
172     20, 
173     20, 
174     5000, 
175     2500, 
176     2500,
177     3000,
178     1500,
179     1500,
180     1000,
181     500,
182     500,
183     500,
184     250,
185     250
186   };
187   Double_t fgSparseDataMin[kNSparseData] = {
188     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.
189   };
190   Double_t fgSparseDataMax[kNSparseData] = {
191     100., 100., 5000., 2500., 2500., 3000., 1500., 1500., 1000., 500., 500., 500., 250., 250.
192   };
193
194
195   fHistoCorrelation = new THnSparseI("hHistoCorrelation", "", kNSparseData, fgSparseDataBins, fgSparseDataMin, fgSparseDataMax);
196   for (Int_t iaxis = 0; iaxis < kNSparseData; iaxis++)
197     fHistoCorrelation->GetAxis(iaxis)->SetTitle(fgkSparseDataTitle[iaxis]);
198   fHistoList->Add(fHistoCorrelation);
199     
200   PostData(1, fHistoList);
201 }
202
203 //________________________________________________________________________
204
205 void AliAnalysisTaskPIDFluctuation::UserExec(Option_t *) 
206 {
207
208   /*
209    * user exec
210    */
211
212   /* get ESD event */
213   AliVEvent *event = InputEvent();
214   if (!event) return;
215
216   /* accept event */
217   if (!AcceptEvent(event)) return;
218
219   /* get centrality object and centrality */
220   AliCentrality *centrality = event->GetCentrality();
221   Float_t cent_v0m = centrality->GetCentralityPercentileUnchecked("V0M");
222   Float_t cent_trk = centrality->GetCentralityPercentileUnchecked("TRK");
223
224   /* init PID */
225   InitPID(event);
226
227   Bool_t pidFlag[AliPID::kSPECIES];
228   Int_t icharge;
229   Int_t chargedCounts[2];
230   Int_t pidCounts[AliPID::kSPECIES][2];
231   for (Int_t i = 0; i < 2; i++) {
232     chargedCounts[i] = 0;
233     for (Int_t ii = 0; ii < 5; ii++) {
234       pidCounts[ii][i] = 0;
235     }
236   }
237
238   /* loop over tracks */
239   for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
240
241     /* get track */
242     AliVTrack *track = (AliVTrack *)event->GetTrack(itrk);
243     if (!track) continue;
244
245     /* accept track */
246     if (!AcceptTrack(track)) continue;
247
248     /* get charge */
249     icharge = track->Charge() > 0 ? 0 : 1;
250     chargedCounts[icharge]++;
251
252     /* make PID */
253     MakePID(track, pidFlag, cent_v0m);
254     for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++)
255       if (pidFlag[ipart])
256         pidCounts[ipart][icharge]++;
257     
258   }
259
260   /* fill histogram */
261   Double_t vsparse[kNSparseData];
262   vsparse[kCent_V0M] = cent_v0m;
263   vsparse[kCent_TRK] = cent_trk;
264   vsparse[kNch] = chargedCounts[0] + chargedCounts[1];
265   vsparse[kNch_plus] = chargedCounts[0];
266   vsparse[kNch_minus] = chargedCounts[1];
267   vsparse[kNpi] = pidCounts[AliPID::kPion][0] + pidCounts[AliPID::kPion][1];
268   vsparse[kNpi_plus] = pidCounts[AliPID::kPion][0];
269   vsparse[kNpi_minus] = pidCounts[AliPID::kPion][1];
270   vsparse[kNka] = pidCounts[AliPID::kKaon][0] + pidCounts[AliPID::kKaon][1];
271   vsparse[kNka_plus] = pidCounts[AliPID::kKaon][0];
272   vsparse[kNka_minus] = pidCounts[AliPID::kKaon][1];
273   vsparse[kNpr] = pidCounts[AliPID::kProton][0] + pidCounts[AliPID::kProton][1];
274   vsparse[kNpr_plus] = pidCounts[AliPID::kProton][0];
275   vsparse[kNpr_minus] = pidCounts[AliPID::kProton][1];
276   fHistoCorrelation->Fill(vsparse);
277
278   PostData(1, fHistoList);
279 }      
280
281 //___________________________________________________________
282
283 Bool_t
284 AliAnalysisTaskPIDFluctuation::AcceptEvent(AliVEvent *event) const
285 {
286   /*
287    * accept event
288    */
289
290   /* fill histo event counter */
291   fHistoEventCounter->Fill(kAllEvents);
292   
293   /* check physics selection */
294   if (!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB)) return kFALSE;
295   fHistoEventCounter->Fill(kPhysicsSelection);
296
297   /* check primary vertex */
298   const AliVVertex *vertex = event->GetPrimaryVertex();
299   if (vertex->GetNContributors() < 1) {
300     /* get ESD vertex SPD */
301     if (event->InheritsFrom("AliESDEvent")) {
302       AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(event);
303       vertex = esdevent->GetPrimaryVertexSPD();
304     }
305     /* get AOD vertex SPD */
306     else if (event->InheritsFrom("AliAODEvent")) {
307       AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(event);
308       vertex = aodevent->GetPrimaryVertexSPD();
309     }
310     if (vertex->GetNContributors() < 1) return kFALSE;
311     fHistoEventCounter->Fill(kPrimaryVertexSPD);
312   }
313   fHistoEventCounter->Fill(kPrimaryVertex);
314
315   /* check vertex position */
316   if (TMath::Abs(vertex->GetZ()) > 10.) return kFALSE;
317   fHistoEventCounter->Fill(kVertexAccepted);
318
319   /* get centrality object and check quality */
320   AliCentrality *centrality = event->GetCentrality();
321   if (centrality->GetQuality() != 0) return kFALSE;
322   fHistoEventCounter->Fill(kGoodCentrality);
323   
324   /* event accepted */
325   fHistoEventCounter->Fill(kAcceptedEvents);
326   return kTRUE;
327 }
328
329 //___________________________________________________________
330
331 Bool_t
332 AliAnalysisTaskPIDFluctuation::AcceptTrack(AliVTrack *track) const
333 {
334   /*
335    * accept track
336    */
337
338   /* check ESD track cuts */
339   if (track->InheritsFrom("AliESDtrack")) {
340     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
341     if (!fESDtrackCuts->AcceptTrack(esdtrack)) return kFALSE;
342   }
343   /* check AOD filter bit */
344   else if (track->InheritsFrom("AliAODTrack")) {
345     AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
346     if (!aodtrack->TestFilterBit(fAODfilterBit)) return kFALSE;
347   }
348
349   /* check eta range */
350   if (track->Eta() < fEtaMin ||
351       track->Eta() > fEtaMax) return kFALSE;
352   /* check pt range */
353   if (track->Pt() < fPtMin ||
354       track->Pt() > fPtMax) return kFALSE;
355
356   /* track accepted */
357   return kTRUE;
358 }
359
360 //___________________________________________________________
361
362 Bool_t
363 AliAnalysisTaskPIDFluctuation::HasTPCPID(AliVTrack *track) const
364 {
365   /*
366    * has TPC PID
367    */
368   
369   /* check PID signal */
370   if (track->GetTPCsignal() <= 0. ||
371       track->GetTPCsignalN() == 0) return kFALSE;
372   return kTRUE;
373   
374 }
375   
376 //___________________________________________________________
377
378 Bool_t
379 AliAnalysisTaskPIDFluctuation::HasTOFPID(AliVTrack *track) const
380 {
381   /*
382    * has TOF PID
383    */
384
385   /* check TOF matched track */
386   if (!(track->GetStatus() & AliESDtrack::kTOFout)||
387       !(track->GetStatus() & AliESDtrack::kTIME)) return kFALSE;
388   return kTRUE;
389
390 }
391
392 //___________________________________________________________
393
394 Double_t
395 AliAnalysisTaskPIDFluctuation::MakeTPCPID(AliVTrack *track, Double_t *nSigma) const
396 {
397   /*
398    * make TPC PID
399    * returns measured dEdx if PID available, otherwise -1.
400    * fills nSigma array with TPC nsigmas for e, mu, pi, K, p
401    */
402   
403   /* check TPC PID */
404   if (!HasTPCPID(track)) return -1.;
405
406   /* get TPC info */
407   Double_t ptpc = track->GetTPCmomentum();
408   Double_t dEdx = track->GetTPCsignal();
409   Int_t dEdxN = track->GetTPCsignalN();
410     
411   /* loop over particles */
412   for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
413     Double_t bethe = fPID->GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)ipart);
414     Double_t diff = dEdx - bethe;
415     Double_t sigma = fPID->GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)ipart);
416     nSigma[ipart] = diff / sigma;
417   }
418   return dEdx;
419
420 }
421
422 //___________________________________________________________
423
424 Double_t
425 AliAnalysisTaskPIDFluctuation::MakeTOFPID(AliVTrack *track, Double_t *nSigma) const
426 {
427   /*
428    * make TOF PID
429    * returns measured beta if PID available, otherwise -1.
430    * fills nSigma array with TOF nsigmas for e, mu, pi, K, p
431    */
432   
433   /* check TOF PID */
434   if (!HasTOFPID(track)) return -1.;
435
436   /* get TOF info */
437   Double_t p = track->P();
438   Double_t time = track->GetTOFsignal() - fPID->GetTOFResponse().GetStartTime(p);
439   Double_t timei[5];
440   track->GetIntegratedTimes(timei);
441   
442   /* loop over particles */
443   for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
444     Double_t timez = time - timei[ipart];
445     Double_t sigma = fPID->GetTOFResponse().GetExpectedSigma(p, timei[ipart], AliPID::ParticleMass(ipart));
446     nSigma[ipart] = timez / sigma;
447   }
448   
449   return timei[0] / time;
450
451 }
452
453 //___________________________________________________________
454
455 void
456 AliAnalysisTaskPIDFluctuation::MakePID(AliVTrack *track, Bool_t *pidFlag, Float_t centrality) const
457 {
458   /*
459    * make PID
460    * fills PID QA plots
461    * fills pidFlag array with PID flags for e, mu, pi, K, p
462    */
463   
464   /* cut definitions
465      (better put them as static variables so they can be changed from outside) */
466   Double_t fgTPCPIDmomcut[AliPID::kSPECIES] = {0., 0., 0.5, 0.5, 0.7};
467   Double_t fgTPCPIDsigmacut[AliPID::kSPECIES] = {0., 0., 2., 2., 2.};
468   Double_t fgTPCPIDcompcut[AliPID::kSPECIES] = {0., 0., 3., 3., 3.};
469   Double_t fgTOFPIDmomcut[AliPID::kSPECIES] = {0., 0., 1.5, 1.5, 2.0};
470   Double_t fgTOFPIDsigmacut[AliPID::kSPECIES] = {0., 0., 2., 2., 2.};
471
472   /* get momentum information */
473   Double_t p = track->P();
474   Double_t pt = track->Pt();
475   Double_t ptpc = track->GetTPCmomentum();
476
477   /* make pid and check if available */
478   Double_t nsigmaTPC[AliPID::kSPECIES];
479   Double_t nsigmaTOF[AliPID::kSPECIES];
480   Double_t dEdx = MakeTPCPID(track, nsigmaTPC);
481   Double_t beta = MakeTOFPID(track, nsigmaTOF);
482   Bool_t hasTPCPID = dEdx > 0.;
483   Bool_t hasTOFPID = beta > 0.;
484
485   /* check PID method */
486   if (fPIDMethod == kTPConly) hasTOFPID = kFALSE; // inhibit TOF PID
487   if (fPIDMethod == kTOFonly) hasTPCPID = kFALSE; // inhibit TPC PID
488
489   /* fill qa histos */
490   fHistoAcceptedTracks->Fill(centrality, pt);
491   if (hasTPCPID)
492     fHistoTPCdEdx->Fill(centrality, ptpc, dEdx);
493   if (hasTOFPID) {
494     fHistoTOFMatchedTracks->Fill(centrality, pt);
495     fHistoTOFbeta->Fill(centrality, p, beta);
496   }
497   if (hasTPCPID && !hasTOFPID)
498     fHistoTPCdEdx_inclusive->Fill(centrality, ptpc, dEdx);
499     
500   /* loop over species */
501   for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
502     
503     /* reset PID flag */
504     pidFlag[ipart] = kFALSE;
505     
506     /* fill qa histos */
507     if (hasTPCPID)
508       fHistoNSigmaTPC[ipart]->Fill(centrality, pt, nsigmaTPC[ipart]);
509     if (hasTOFPID)
510       fHistoNSigmaTOF[ipart]->Fill(centrality, pt, nsigmaTOF[ipart]);
511     
512     /* combined PID cuts */
513     if (hasTPCPID && hasTOFPID) {
514       if (pt < fgTOFPIDmomcut[ipart] &&
515           TMath::Abs(nsigmaTOF[ipart]) < fgTOFPIDsigmacut[ipart] &&
516           TMath::Abs(nsigmaTPC[ipart]) < fgTPCPIDcompcut[ipart]) {
517         fHistoTOFbeta_selected[ipart]->Fill(centrality, p, beta);
518         pidFlag[ipart] = kTRUE;
519       }
520     }
521     /* TPC-only PID cuts */
522     else if (hasTPCPID && !hasTOFPID) {
523       if (pt < fgTPCPIDmomcut[ipart] &&
524           TMath::Abs(nsigmaTPC[ipart]) < fgTPCPIDsigmacut[ipart]) { 
525         fHistoTPCdEdx_selected[ipart]->Fill(centrality, ptpc, dEdx);
526         pidFlag[ipart] = kTRUE;
527       }
528     }
529     /* TOF-only PID cuts */
530     else if (!hasTPCPID && hasTOFPID) {
531       if (pt < fgTOFPIDmomcut[ipart] &&
532           TMath::Abs(nsigmaTOF[ipart]) < fgTOFPIDsigmacut[ipart]) {
533         fHistoTOFbeta_selected[ipart]->Fill(centrality, p, beta);
534         pidFlag[ipart] = kTRUE;
535       }
536     }
537     
538   } /* end of loop over species */
539
540 }
541
542 //___________________________________________________________
543
544 void
545 AliAnalysisTaskPIDFluctuation::InitPID(AliVEvent *event)
546 {
547   /*
548    * init PID
549    */
550
551   /* create PID object if not there yet */
552   if (!fPID) {
553
554     /* instance object */
555     Bool_t mcFlag = kFALSE; /*** WARNING: check whether is MC ***/
556     if (event->InheritsFrom("AliESDEvent"))
557       fPID = new AliESDpid(mcFlag);
558     else if (event->InheritsFrom("AliAODEvent"))
559       fPID = new AliAODpidUtil(mcFlag);
560
561     /* set OADB path */
562     fPID->SetOADBPath("$ALICE_ROOT/OADB");
563   }
564
565   /* init ESD PID */
566   Int_t recoPass = 2; /*** WARNING: need to set the recoPass somehow better ***/
567   fPID->InitialiseEvent(event, recoPass); /* warning: this apparently sets TOF time        
568                                               * resolution to some hardcoded value,     
569                                               * therefore we have to set correct   
570                                               * resolution value after this call */
571
572   /* set TOF response */
573   Double_t tofReso = 85.; /* ps */ /*** WARNING: need to set tofReso somehow better ***/
574   fPID->GetTOFResponse().SetTimeResolution(tofReso);
575   fPID->GetTOFResponse().SetTrackParameter(0, 0.007);
576   fPID->GetTOFResponse().SetTrackParameter(1, 0.007);
577   fPID->GetTOFResponse().SetTrackParameter(2, 0.0);
578   fPID->GetTOFResponse().SetTrackParameter(3, 30);
579   
580 }
581
582 //___________________________________________________________
583
584 void
585 AliAnalysisTaskPIDFluctuation::MeasureNuDyn(const Char_t *filename, Int_t i1, Int_t i2, Int_t centralityEstimator)
586 {
587
588   /*
589    * measure nu
590    */
591   
592   printf("MeasureNuDyn: running for %s vs. %s\n", fgkSparseDataName[i1], fgkSparseDataName[i2]);
593
594   /* get data */
595   TFile *filein = TFile::Open(filename);
596   /* output */
597   TFile *fileout = TFile::Open(Form("MeasureNuDyn_%s_%s.%s", fgkSparseDataName[i1], fgkSparseDataName[i2], filename), "RECREATE");
598
599   /* loop over available containers */
600   TList *keylist = filein->GetListOfKeys();
601   for (Int_t ikey = 0; ikey < keylist->GetEntries(); ikey++) {
602     
603     /* get key and check */
604     TKey *key = (TKey *)keylist->At(ikey);
605     TString contname = key->GetName();
606     if (!contname.BeginsWith("PIDFluctuation")) continue;
607
608     /* get data */
609     TList *list = (TList *)filein->Get(contname.Data());
610     THnSparse *hsparse = (THnSparse *)list->FindObject("hHistoCorrelation");
611     
612     /* create output directory and cd there */
613     fileout->mkdir(contname.Data());
614     fileout->cd(contname.Data());
615
616     /* loop over centralities */
617     Double_t centBin[kNCentralityBins + 1] = {
618       0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90.
619     };
620     TH1D *hNu = new TH1D("hNu", ";centrality percentile;#nu", kNCentralityBins, centBin);
621     TH1D *hNuStat = new TH1D("hNuStat", ";centrality percentile;#nu_{stat}", kNCentralityBins, centBin);
622     TH1D *hNuDyn = new TH1D("hNuDyn", ";centrality percentile;#nu_{dyn}", kNCentralityBins, centBin);
623     for (Int_t icent = 0; icent < kNCentralityBins; icent++) {
624       
625       /* select centrality range */
626       hsparse->GetAxis(centralityEstimator)->SetRangeUser(centBin[icent] + 0.001, centBin[icent + 1] - 0.001);
627       /* projections */
628       TH1D *hcent = hsparse->Projection(centralityEstimator);
629       TH1D *h1 = hsparse->Projection(i1);
630       TH1D *h2 = hsparse->Projection(i2);
631       TH2D *hcorr = hsparse->Projection(i2, i1);
632       TH1D *hnu = new TH1D("hnu", "", 2000, -1., 1.);
633       
634       Double_t n, a, b, nev, meanx, rmsx, meany, rmsy, meannu, rmsnu;
635       Double_t meanx_err, meany_err, meannu_err;
636       
637       /* compute mean values */
638       nev = 0.; meanx = 0.; meany = 0.;
639       for (Int_t ibinx = 0; ibinx < hcorr->GetNbinsX(); ibinx++)
640         for (Int_t ibiny = 0; ibiny < hcorr->GetNbinsY(); ibiny++) {
641           n = hcorr->GetBinContent(ibinx + 1, ibiny + 1);
642           if (n <= 0.) continue;
643           meanx += n * ibinx;
644           meany += n * ibiny;
645           nev  += n;
646         }
647       meanx /= nev;
648       meany /= nev;
649       //    printf("nev = %f, meanx = %f, meany = %f\n", nev, meanx, meany);
650       
651       /* compute RMS values */
652       nev = 0.; rmsx = 0.; rmsy = 0.;
653       for (Int_t ibinx = 0; ibinx < hcorr->GetNbinsX(); ibinx++)
654         for (Int_t ibiny = 0; ibiny < hcorr->GetNbinsY(); ibiny++) {
655           n = hcorr->GetBinContent(ibinx + 1, ibiny + 1);
656           if (n <= 0.) continue;
657           a = ibinx - meanx;
658           rmsx += n * a * a;
659           a = ibiny - meany;
660           rmsy += n * a * a;
661           nev  += n;
662         }
663       rmsx /= nev;
664       rmsx = TMath::Sqrt(rmsx);
665       rmsy /= nev;
666       rmsy = TMath::Sqrt(rmsy);
667       //    printf("nev = %f, rmsx = %f, rmsy = %f\n", nev, rmsx, rmsy);
668       meanx_err = rmsx / TMath::Sqrt(nev);
669       meany_err = rmsy / TMath::Sqrt(nev);
670       //    printf("nev = %f, meanx_err = %f, meany_err = %f\n", nev, meanx_err, meany_err);
671       
672       /* compute mean nu */
673       nev = 0.; meannu = 0.;
674       for (Int_t ibinx = 0; ibinx < hcorr->GetNbinsX(); ibinx++)
675         for (Int_t ibiny = 0; ibiny < hcorr->GetNbinsY(); ibiny++) {
676           n = hcorr->GetBinContent(ibinx + 1, ibiny + 1);
677           if (n <= 0.) continue;
678           a = ibinx / meanx - ibiny / meany;
679           meannu += n * a * a;
680           hnu->Fill(a * a, n);
681           nev  += n;
682         }
683       meannu /= nev;
684       //    printf("nev = %f, meannu = %f\n", nev, meannu);
685       
686       /* compute RMS nu */
687       nev = 0.; rmsnu = 0.;
688       for (Int_t ibinx = 0; ibinx < hcorr->GetNbinsX(); ibinx++)
689         for (Int_t ibiny = 0; ibiny < hcorr->GetNbinsY(); ibiny++) {
690           n = hcorr->GetBinContent(ibinx + 1, ibiny + 1);
691           if (n <= 0.) continue;
692           a = ibinx / meanx - ibiny / meany;
693           b = a * a - meannu;
694           rmsnu += n * b * b;
695           nev  += n;
696         }
697       rmsnu /= nev;
698       rmsnu = TMath::Sqrt(rmsnu);
699       //    printf("nev = %f, rmsnu = %f\n", nev, rmsnu);
700       meannu_err = rmsnu / TMath::Sqrt(nev);
701       //    printf("nev = %f, meannu_err = %f\n", nev, meannu_err);
702       
703       /* final calculations */
704       Double_t nu = meannu;
705       Double_t nu_err = meannu_err;
706       Double_t nu_stat = 1. / meanx + 1. / meany;
707       Double_t meanx4 = meanx * meanx * meanx * meanx;
708       Double_t meanx_err2 = meanx_err * meanx_err;
709       Double_t meany4 = meany * meany * meany * meany;
710       Double_t meany_err2 = meany_err * meany_err;
711       Double_t nu_stat_err = TMath::Sqrt(meanx_err2 / meanx4 + meany_err2 / meany4);
712       Double_t nu_dyn = nu - nu_stat;
713       Double_t nu_dyn_err = TMath::Sqrt(nu_err * nu_err + nu_stat_err * nu_stat_err);
714       
715       /* setup final plots */
716       hNu->SetBinContent(icent + 1, nu);
717       hNu->SetBinError(icent + 1, nu_err);
718       hNuStat->SetBinContent(icent + 1, nu_stat);
719       hNuStat->SetBinError(icent + 1, nu_stat_err);
720       hNuDyn->SetBinContent(icent + 1, nu_dyn);
721       hNuDyn->SetBinError(icent + 1, nu_dyn_err);
722       
723       hcent->Write(Form("hcent_cent%d", icent));
724       h1->Write(Form("h1_cent%d", icent));
725       h2->Write(Form("h2_cent%d", icent));
726       hcorr->Write(Form("hcorr_cent%d", icent));
727       hnu->Write(Form("hnu_cent%d", icent));
728       
729       /* clean-up */
730       delete hcent;
731       delete h1;
732       delete h2;
733       delete hcorr;
734       delete hnu;
735     }
736
737     hNu->Write();
738     hNuStat->Write();
739     hNuDyn->Write();
740   }
741
742   fileout->Close();
743   filein->Close();
744
745 }
746
747 //___________________________________________________________