Added class AliForwarddNdetaTask to do the dN/deta
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwarddNdetaTask.cxx
1 //====================================================================
2 #include "AliForwarddNdetaTask.h"
3 #include <TMath.h>
4 #include <TH2D.h>
5 #include <TH1D.h>
6 #include <THStack.h>
7 #include <TList.h>
8 #include <AliAnalysisManager.h>
9 #include <AliAODEvent.h>
10 #include <AliAODHandler.h>
11 #include <AliAODInputHandler.h>
12 #include "AliForwardUtil.h"
13
14 //____________________________________________________________________
15 AliForwarddNdetaTask::AliForwarddNdetaTask()
16   : AliAnalysisTaskSE(), 
17     fSumForward(0),     //  Sum of histograms 
18     fSumForwardMC(0),   //  Sum of MC histograms (if any)
19     fSumPrimary(0),     //  Sum of primary histograms
20     fSumCentral(0),     //  Sum of central histograms
21     fCentral(0),        // Cache of central histogram
22     fPrimary(0),        // Cache of primary histogram
23     fSums(0),           // Container of sums 
24     fOutput(0),         // Container of outputs 
25     fTriggers(0),       // Histogram of triggers 
26     fVtxMin(0),         // Minimum v_z
27     fVtxMax(0),         // Maximum v_z
28     fTriggerMask(0),    // Trigger mask 
29     fRebin(0),          // Rebinning factor 
30     fCutEdges(false),
31     fSNNString(0),
32     fSysString(0)
33 {}
34
35 //____________________________________________________________________
36 AliForwarddNdetaTask::AliForwarddNdetaTask(const char* /* name */)
37   : AliAnalysisTaskSE("Forward"), 
38     fSumForward(0),     //  Sum of histograms 
39     fSumForwardMC(0),   //  Sum of MC histograms (if any)
40     fSumPrimary(0),     //  Sum of primary histograms
41     fSumCentral(0),     //  Sum of central histograms
42     fCentral(0),        // Cache of central histogram
43     fPrimary(0),        // Cache of primary histogram
44     fSums(0),           // Container of sums 
45     fOutput(0),         // Container of outputs 
46     fTriggers(0),       // Histogram of triggers 
47     fVtxMin(-10),       // Minimum v_z
48     fVtxMax(10),        // Maximum v_z
49     fTriggerMask(AliAODForwardMult::kInel), 
50     fRebin(5),          // Rebinning factor 
51     fCutEdges(false),
52     fSNNString(0),
53     fSysString(0)
54 {
55   // Output slot #1 writes into a TH1 container
56   DefineOutput(1, TList::Class()); 
57   DefineOutput(2, TList::Class()); 
58 }
59
60 //____________________________________________________________________
61 AliForwarddNdetaTask::AliForwarddNdetaTask(const AliForwarddNdetaTask& o)
62   : AliAnalysisTaskSE(o), 
63     fSumForward(o.fSumForward),         // TH2D* -  Sum of histograms 
64     fSumForwardMC(o.fSumForwardMC),             // TH2D* -  Sum of MC histograms (if any)
65     fSumPrimary(o.fSumPrimary),         // TH2D* -  Sum of primary histograms
66     fSumCentral(o.fSumCentral),         // TH2D* -  Sum of central histograms
67     fCentral(o.fCentral),       //! Cache of central histogram
68     fPrimary(o.fPrimary),       //! Cache of primary histogram
69     fSums(o.fSums),             // TList* - Container of sums 
70     fOutput(o.fOutput),         // TList* - Container of outputs 
71     fTriggers(o.fTriggers),             // TH1D* - Histogram of triggers 
72     fVtxMin(o.fVtxMin),         // Double_t - Minimum v_z
73     fVtxMax(o.fVtxMax),         // Double_t - Maximum v_z
74     fTriggerMask(o.fTriggerMask),               // Int_t - Trigger mask 
75     fRebin(o.fRebin),           // Int_t - Rebinning factor 
76     fCutEdges(o.fCutEdges),             // Bool_t - Whether to cut edges when rebinning
77     fSNNString(o.fSNNString),           // TNamed* - 
78     fSysString(o.fSysString)            // TNamed* - 
79 {}
80
81 //____________________________________________________________________
82 AliForwarddNdetaTask::~AliForwarddNdetaTask()
83 {
84   if (fSums) { 
85     fSums->Delete();
86     delete fSums;
87     fSums = 0;
88   }
89   if (fOutputs) { 
90     fOutputs->Delete();
91     delete fOutputs;
92     fOutputs = 0;
93   }
94   if (fCentral)    delete fCentral;
95   if (fPrimary)    delete fPrimary;
96 }
97
98 //________________________________________________________________________
99 void 
100 AliForwarddNdetaTask::SetTriggerMask(const char* mask)
101 {
102   UShort_t    trgMask = 0;
103   TString     trgs(mask);
104   trgs.ToUpper();
105   TObjString* trg;
106   TIter       next(trgs.Tokenize(" ,|"));
107   while ((trg = static_cast<TObjString*>(next()))) { 
108     TString s(trg->GetString());
109     if      (s.IsNull()) continue;
110     if      (s.CompareTo("INEL")  == 0) trgMask = AliAODForwardMult::kInel;
111     else if (s.CompareTo("INEL>0")== 0) trgMask = AliAODForwardMult::kInelGt0;
112     else if (s.CompareTo("NSD")   == 0) trgMask = AliAODForwardMult::kNSD;
113     else 
114       Warning("SetTriggerMask", "Unknown trigger %s", s.Data());
115   }
116   if (trgMask == 0) trgMask = 1;
117   SetTriggerMask(trgMask);
118 }
119
120 //________________________________________________________________________
121 void 
122 AliForwarddNdetaTask::UserCreateOutputObjects()
123 {
124   // Create histograms
125   // Called once (on the worker node)
126
127   fOutput = new TList;
128   fOutput->SetName(Form("%s_result", GetName()));
129   fOutput->SetOwner();
130
131   fSums = new TList;
132   fSums->SetName(Form("%s_sums", GetName()));
133   fSums->SetOwner();
134
135
136   fTriggers = new TH1D("triggers", "Number of triggers", 
137                        kAccepted, 1, kAccepted);
138   fTriggers->SetYTitle("# of events");
139   fTriggers->GetXaxis()->SetBinLabel(kAll,         "All events");
140   fTriggers->GetXaxis()->SetBinLabel(kB,           "w/B trigger");
141   fTriggers->GetXaxis()->SetBinLabel(kA,           "w/A trigger");
142   fTriggers->GetXaxis()->SetBinLabel(kC,           "w/C trigger");
143   fTriggers->GetXaxis()->SetBinLabel(kE,           "w/E trigger");
144   fTriggers->GetXaxis()->SetBinLabel(kMB,          "w/Collision trigger");
145   fTriggers->GetXaxis()->SetBinLabel(kWithVertex,  "w/Vertex");
146   fTriggers->GetXaxis()->SetBinLabel(kWithTrigger, "w/Selected trigger");
147   fTriggers->GetXaxis()->SetBinLabel(kAccepted,    "Accepted by cut");
148   fTriggers->GetXaxis()->SetNdivisions(kAccepted, false);
149   fTriggers->SetFillColor(kRed+1);
150   fTriggers->SetFillStyle(3001);
151   fTriggers->SetStats(0);
152   fSums->Add(fTriggers);
153
154   fSNNString = new TNamed("sNN", "");
155   fSysString = new TNamed("sys", "");
156   fSums->Add(fSNNString);
157   fSums->Add(fSysString);
158
159   // Check that we have an AOD input handler 
160   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
161   AliAODInputHandler* ah = 
162     dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler());
163   if (!ah) AliFatal("No AOD input handler set in analysis manager");
164
165   // Post data for ALL output slots >0 here, to get at least an empty histogram
166   PostData(1, fSums); 
167 }
168
169 //____________________________________________________________________
170 TH2D*
171 AliForwarddNdetaTask::CloneHist(TH2D* in, const char* name) 
172 {
173   if (!in) return 0;
174   TH2D* ret = static_cast<TH2D*>(in->Clone(name));
175   ret->SetDirectory(0);
176   ret->Sumw2();
177   fSums->Add(ret);
178
179   return ret;
180 }
181
182 //____________________________________________________________________
183 void 
184 AliForwarddNdetaTask::UserExec(Option_t *) 
185 {
186   // Main loop
187   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
188   if (!aod) {
189     AliError("Cannot get the AOD event");
190     return;
191   }  
192
193   // Get objects from the event structure 
194   TObject* oForward   = aod->FindListObject("Forward");
195   TObject* oForwardMC = aod->FindListObject("ForwardMC");
196   TObject* oPrimary   = aod->FindListObject("primary");
197   TObject* oCentral   = aod->FindListObject("Central");
198
199   // We should have a forward object at least 
200   if (!oForward) { 
201     AliWarning("No Forward object found AOD");
202     return;
203   }
204
205   // Cast to good types 
206   AliAODForwardMult* forward   = static_cast<AliAODForwardMult*>(oForward);
207   AliAODForwardMult* forwardMC = static_cast<AliAODForwardMult*>(oForwardMC);
208   TH2D*              primary   = static_cast<TH2D*>(oPrimary);
209   TH2D*              central   = static_cast<TH2D*>(oCentral);
210
211   static bool first = true;
212   if (first) { 
213     UShort_t sNN = forward->GetSNN();
214     UShort_t sys = forward->GetSystem();
215     fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
216     fSNNString->SetUniqueID(sNN);
217     fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
218     fSysString->SetUniqueID(sys);
219     first = false;
220   }
221
222   // Create our sum histograms 
223   if (!fSumForward) fSumForward = CloneHist(&forward->GetHistogram(),"forward");
224   if (!fSumForwardMC && forwardMC) 
225     fSumForwardMC = CloneHist(&forwardMC->GetHistogram(),"forwardMC");
226   if (!fSumPrimary && primary) fSumPrimary = CloneHist(primary, "primary");
227   if (!fSumCentral && central) fSumCentral = CloneHist(central, "central");
228
229   // Add contribtion from MC 
230   if (primary) fSumPrimary->Add(primary);
231   
232   // Count event 
233   fTriggers->AddBinContent(kAll);
234   if (forward->IsTriggerBits(AliAODForwardMult::kB)) 
235     fTriggers->AddBinContent(kB);
236   if (forward->IsTriggerBits(AliAODForwardMult::kA)) 
237     fTriggers->AddBinContent(kA);
238   if (forward->IsTriggerBits(AliAODForwardMult::kC)) 
239     fTriggers->AddBinContent(kC);
240   if (forward->IsTriggerBits(AliAODForwardMult::kE)) 
241     fTriggers->AddBinContent(kE);
242   if (forward->IsTriggerBits(AliAODForwardMult::kInel)) 
243     fTriggers->AddBinContent(kMB);
244
245   // Check if we have an event of interest. 
246   if (!forward->IsTriggerBits(fTriggerMask)) return;
247   fTriggers->AddBinContent(kWithTrigger);
248
249   // Check that we have a valid vertex
250   if (!forward->HasIpZ()) return;
251   fTriggers->AddBinContent(kWithVertex);
252
253   // Check that vertex is within cuts 
254   if (!forward->InRange(fVtxMin, fVtxMax)) return;
255   fTriggers->AddBinContent(kAccepted);
256
257   // Add contribution 
258   fSumForward->Add(&(forward->GetHistogram()));
259   if (fSumForwardMC) fSumForwardMC->Add(&(forwardMC->GetHistogram()));
260   if (fSumPrimary)   fSumPrimary->Add(primary);
261   if (fSumCentral)   fSumCentral->Add(central);
262
263   PostData(1, fSums);
264 }
265
266 //________________________________________________________________________
267 void 
268 AliForwarddNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
269                                           const char* title, const char* ytitle)
270 {
271   h->SetTitle(title);
272   h->SetMarkerColor(colour);
273   h->SetMarkerStyle(marker);
274   h->SetMarkerSize(1);
275   h->SetFillStyle(0);
276   h->SetYTitle(ytitle);
277   h->GetXaxis()->SetTitleFont(132);
278   h->GetXaxis()->SetLabelFont(132);
279   h->GetXaxis()->SetNdivisions(10);
280   h->GetYaxis()->SetTitleFont(132);
281   h->GetYaxis()->SetLabelFont(132);
282   h->GetYaxis()->SetNdivisions(10);
283   h->GetYaxis()->SetDecimals();
284   h->SetStats(0);
285 }
286   
287 //________________________________________________________________________
288 void 
289 AliForwarddNdetaTask::Terminate(Option_t *) 
290 {
291   // Draw result to screen, or perform fitting, normalizations Called
292   // once at the end of the query
293         
294   fSums = dynamic_cast<TList*> (GetOutputData(1));
295   if(!fOutput) {
296     AliError("Could not retrieve TList fSums"); 
297     return; 
298   }
299   
300   if (!fOutput) { 
301     fOutput = new TList;
302     fOutput->SetName(Form("%s_result", GetName()));
303     fOutput->SetOwner();
304   }
305
306   fSumForward     = static_cast<TH2D*>(fSums->FindObject("forward"));
307   fSumForwardMC   = static_cast<TH2D*>(fSums->FindObject("forwardMC"));
308   fSumPrimary     = static_cast<TH2D*>(fSums->FindObject("primary"));
309   fSumCentral     = static_cast<TH2D*>(fSums->FindObject("central"));
310   fTriggers       = static_cast<TH1D*>(fSums->FindObject("triggers"));
311   fSNNString      = static_cast<TNamed*>(fSums->FindObject("sNN"));
312   fSysString      = static_cast<TNamed*>(fSums->FindObject("sys"));
313
314   if (!fTriggers) { 
315     AliError("Couldn't find histogram 'triggers' in list");
316     return;
317   }
318   if (!fSumForward) { 
319     AliError("Couldn't find histogram 'forward' in list");
320     return;
321   }
322
323   Int_t nAll        = Int_t(fTriggers->GetBinContent(kAll));
324   Int_t nB          = Int_t(fTriggers->GetBinContent(kB));
325   Int_t nA          = Int_t(fTriggers->GetBinContent(kA));
326   Int_t nC          = Int_t(fTriggers->GetBinContent(kC));
327   Int_t nE          = Int_t(fTriggers->GetBinContent(kE));
328   Int_t nMB         = Int_t(fTriggers->GetBinContent(kMB));
329   Int_t nTriggered  = Int_t(fTriggers->GetBinContent(kWithTrigger));
330   Int_t nWithVertex = Int_t(fTriggers->GetBinContent(kWithVertex));
331   Int_t nAccepted   = Int_t(fTriggers->GetBinContent(kAccepted));
332   Int_t nGood       = nB - nA - nC + 2 * nE;
333   Double_t vtxEff   = Double_t(nMB) / nTriggered * Double_t(nAccepted) / nGood;
334   AliInfo(Form("Total of %9d events\n"
335                "                   of these %9d are minimum bias\n"
336                "                   of these %9d has a %s trigger\n" 
337                "                   of these %9d has a vertex\n" 
338                "                   of these %9d were in [%+4.1f,%+4.1f]cm\n"
339                "                   Triggers by type:\n"
340                "                     B   = %9d\n"
341                "                     A|C = %9d (%9d+%-9d)\n"
342                "                     E   = %9d\n"
343                "                   Implies %9d good triggers\n"
344                "                   Vertex efficiency: %f",
345                nAll, nMB, nTriggered, 
346                AliAODForwardMult::GetTriggerString(fTriggerMask),
347                nWithVertex, nAccepted,
348                fVtxMin, fVtxMax, 
349                nB, nA+nC, nA, nC, nE, nGood, vtxEff));
350   
351   // Get acceptance normalisation from underflow bins 
352   TH1D* normForward   = fSumForward->ProjectionX("normForward", 0, 1, "");
353   // Project onto eta axis - _ignoring_underflow_bins_!
354   TH1D* dndetaForward = fSumForward->ProjectionX("dndetaForward", 1, -1, "e");
355   // Normalize to the acceptance 
356   dndetaForward->Divide(normForward);
357   // Scale by the vertex efficiency 
358   dndetaForward->Scale(vtxEff, "width");
359
360   SetHistogramAttributes(dndetaForward, kRed+1, 20, "ALICE Forward");
361   SetHistogramAttributes(normForward, kRed+1, 20, "ALICE Forward", "Normalisation");
362
363   fOutput->Add(fTriggers->Clone());
364   fOutput->Add(fSNNString->Clone());
365   fOutput->Add(fSysString->Clone());
366   fOutput->Add(dndetaForward);
367   fOutput->Add(normForward);
368   fOutput->Add(Rebin(dndetaForward));
369
370   if (fSumForwardMC) { 
371     // Get acceptance normalisation from underflow bins 
372     TH1D* normForwardMC   = fSumForwardMC->ProjectionX("normForwardMC", 0, 1, "");
373     // Project onto eta axis - _ignoring_underflow_bins_!
374     TH1D* dndetaForwardMC = fSumForwardMC->ProjectionX("dndetaForwardMC", 1, -1, "e");
375     // Normalize to the acceptance 
376     dndetaForwardMC->Divide(normForwardMC);
377     // Scale by the vertex efficiency 
378     dndetaForwardMC->Scale(vtxEff, "width");
379
380     SetHistogramAttributes(dndetaForwardMC, kRed+3, 21, "ALICE Forward (MC)");
381     SetHistogramAttributes(normForwardMC, kRed+3, 21, "ALICE Forward (MC)", 
382                            "Normalisation");
383
384     fOutput->Add(dndetaForwardMC);
385     fOutput->Add(normForwardMC);
386     fOutput->Add(Rebin(dndetaForwardMC));
387   }
388
389   if (fSumPrimary) { 
390     TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",-1,-1,"e");
391     dndetaTruth->Scale(1./nAll, "width");
392
393     SetHistogramAttributes(dndetaTruth, kGray+3, 22, "Monte-Carlo truth");
394
395     fOutput->Add(dndetaTruth);
396     fOutput->Add(Rebin(dndetaTruth));
397   }
398   if (fSumCentral) { 
399     TH1D* dndetaCentral = fSumCentral->ProjectionX("dndetaCentral",-1,-1,"e");
400     dndetaCentral->Scale(1./nGood, "width");
401
402     SetHistogramAttributes(dndetaCentral, kGray+3, 22, "ALICE Central - track(let)s");
403
404     dndetaCentral->GetXaxis()->SetRangeUser(-1,1);
405     fOutput->Add(dndetaCentral);
406     fOutput->Add(Rebin(dndetaCentral));
407   }
408
409   TNamed* trigString = 
410     new TNamed("trigString", AliAODForwardMult::GetTriggerString(fTriggerMask));
411   trigString->SetUniqueID(fTriggerMask);
412   fOutput->Add(trigString);
413
414   TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
415   vtxAxis->SetName("vtxAxis");
416   vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
417   fOutput->Add(vtxAxis);
418
419   // If only there was a away to get sqrt{s_NN} and beam type 
420   
421
422   PostData(2, fOutput);
423 }
424
425 //________________________________________________________________________
426 TH1D*
427 AliForwarddNdetaTask::Rebin(const TH1D* h) const
428 {
429   if (fRebin <= 1) return 0;
430
431   Int_t nBins = h->GetNbinsX();
432   if(nBins % fRebin != 0) {
433     Warning("Rebin", "Rebin factor %d is not a devisor of current number "
434             "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
435     return 0;
436   }
437     
438   // Make a copy 
439   TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d", 
440                                                h->GetName(), fRebin)));
441   tmp->Rebin(fRebin);
442   tmp->SetDirectory(0);
443
444   // The new number of bins 
445   Int_t nBinsNew = nBins / fRebin;
446   for(Int_t i = 1;i<= nBinsNew; i++) {
447     Double_t content = 0;
448     Double_t sumw    = 0;
449     Double_t wsum    = 0;
450     Int_t    nbins   = 0;
451     for(Int_t j = 1; j<=fRebin;j++) {
452       Int_t    bin = (i-1)*fRebin + j;
453       Double_t c   =  h->GetBinContent(bin);
454
455       if (c <= 0) continue;
456
457       if (fCutEdges) {
458         if (h->GetBinContent(bin+1)<=0 || 
459             h->GetBinContent(bin-1)) {
460           Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)", 
461                   bin, c, h->GetName(), 
462                   bin+1, h->GetBinContent(bin+1), 
463                   bin-1, h->GetBinContent(bin-1));
464           continue;
465         }       
466       }
467       Double_t e =  h->GetBinError(bin);
468       Double_t w =  1 / (e*e); // 1/c/c
469       content    += c;
470       sumw       += w;
471       wsum       += w * c;
472       nbins++;
473     }
474       
475     if(content > 0 && nbins > 0) {
476       tmp->SetBinContent(i, wsum / sumw);
477       tmp->SetBinError(i,1./TMath::Sqrt(sumw));
478     }
479   }
480   
481   return tmp;
482 }