Adding analysis task for the centrality trigger calibration.
[u/mrichter/AliRoot.git] / PWG1 / VZERO / AliAnaVZEROTrigger.cxx
1 #include <stdio.h>
2 #include <stdlib.h>
3
4 #include "TH1F.h"
5 #include "TH2F.h"
6 #include "TChain.h"
7 #include "TFile.h"
8 #include "TParameter.h"
9
10 #include "AliAnalysisTask.h"
11 #include "AliAnalysisManager.h"
12
13 #include "AliESDEvent.h"
14 #include "AliESDfriend.h"
15 #include "AliESDInputHandler.h"
16 #include "AliAnaVZEROTrigger.h"
17 #include "AliCentrality.h"
18
19 // VZERO includes
20 #include "AliESDVZERO.h"
21
22 ClassImp(AliAnaVZEROTrigger)
23
24 AliAnaVZEROTrigger::AliAnaVZEROTrigger() 
25   : AliAnalysisTaskSE("AliAnaVZEROTrigger"), fESD(0), fOutputList(0),
26   fMinThr(812.),
27   fMaxThr(6988.),
28   fRatio(2.08),
29   fNThr(20),
30   fMBTrigName("CPBI"),
31   fV0Percent(0),
32   fV0PercentAll(0),
33   fZvtx(0),
34   fV0Mult1d(0),
35   fV0Charge2d(0),
36   fV0Charge2dPercent(0),
37   fV0Charge2dAll(0),
38   fV0PercentBins(0),
39   fV0PercentBinsAll(0),
40   fV0Cent(0),
41   fV0CentAll(0),
42   fV0SemiCent(0),
43   fV0SemiCentAll(0),
44   fV0CentHw(0),
45   fV0CentHwAll(0),
46   fV0SemiCentHw(0),
47   fV0SemiCentHwAll(0),
48   fV0CentTr(0),
49   fV0CentTrAll(0),
50   fV0SemiCentTr(0),
51   fV0SemiCentTrAll(0),
52   fV0Percent63(0),
53   fV0Percent63All(0),
54   fV0MultAll(0),
55   fV0Mult63(0)
56 {
57   // Constructor
58   // Init default thr
59   fCentCuts[0] = 4152;
60   fCentCuts[1] = 8637;
61   fSemiCentCuts[0] = 617;
62   fSemiCentCuts[1] = 1283;
63
64   // Define input and output slots here
65   // Input slot #0 works with a TChain
66   DefineInput(0, TChain::Class());
67   // Output slot #0 id reserved by the base class for AOD
68   // Output slot #1 writes into a TH1 container
69   DefineOutput(1, TList::Class());
70 }
71
72 //________________________________________________________________________
73 AliAnaVZEROTrigger::AliAnaVZEROTrigger(const char *name) 
74   : AliAnalysisTaskSE(name), fESD(0), fOutputList(0),
75   fMinThr(812.),
76   fMaxThr(6988.),
77   fRatio(2.08),
78   fNThr(20),
79   fMBTrigName("CPBI"),
80   fV0Percent(0),
81   fV0PercentAll(0),
82   fZvtx(0),
83   fV0Mult1d(0),
84   fV0Charge2d(0),
85   fV0Charge2dPercent(0),
86   fV0Charge2dAll(0),
87   fV0PercentBins(0),
88   fV0PercentBinsAll(0),
89   fV0Cent(0),
90   fV0CentAll(0),
91   fV0SemiCent(0),
92   fV0SemiCentAll(0),
93   fV0CentHw(0),
94   fV0CentHwAll(0),
95   fV0SemiCentHw(0),
96   fV0SemiCentHwAll(0),
97   fV0CentTr(0),
98   fV0CentTrAll(0),
99   fV0SemiCentTr(0),
100   fV0SemiCentTrAll(0),
101   fV0Percent63(0),
102   fV0Percent63All(0),
103   fV0MultAll(0),
104   fV0Mult63(0)
105 {
106   // Constructor
107   // Init default thr
108   fCentCuts[0] = 4152;
109   fCentCuts[1] = 8637;
110   fSemiCentCuts[0] = 617;
111   fSemiCentCuts[1] = 1283;
112
113   // Define input and output slots here
114   // Input slot #0 works with a TChain
115   DefineInput(0, TChain::Class());
116   // Output slot #0 id reserved by the base class for AOD
117   // Output slot #1 writes into a TH1 container
118   DefineOutput(1, TList::Class());
119 }
120
121 //________________________________________________________________________
122 void AliAnaVZEROTrigger::Setup(const char *filename)
123 {
124   // Open the config file and
125   // read the parameters
126   FILE *fin = fopen(filename,"r");
127   if (!fin) {
128     AliInfo(Form("File %s is not found running with the default parameters.",filename));
129   }
130   else {
131     Int_t res = fscanf(fin,"%f %f %f %d %f %f %f %f",
132                        &fMinThr,&fMaxThr,&fRatio,&fNThr,
133                        &fSemiCentCuts[0],&fSemiCentCuts[1],&fCentCuts[0],&fCentCuts[1]);
134     if(res!=8) {
135       AliFatal("Failed to get values from the config file.\n");
136     }
137     fclose(fin);
138   }
139
140   AliInfo(Form("MinThr=%.1f MaxThr=%.1f Ratio=%.4f NThr=%d",fMinThr,fMaxThr,fRatio,fNThr));
141   AliInfo(Form("CTA1=%.1f CTC1=%.1f CTA2=%.1f CTC2=%.1f",fSemiCentCuts[0],fSemiCentCuts[1],fCentCuts[0],fCentCuts[1]));
142 }
143
144 //________________________________________________________________________
145 void AliAnaVZEROTrigger::UserCreateOutputObjects()
146 {
147   // Create histograms
148   // Called once
149
150   fOutputList = new TList();
151   fOutputList->SetOwner(kTRUE);
152
153   fV0Percent    = new TH1F("fV0Percent","Centrality percentile based on v0 mult",400,0,100);
154   fOutputList->Add(fV0Percent);
155   fV0PercentAll    = new TH1F("fV0PercentAll","Centrality percentile based on v0 mult (no event selection)",400,0,100);
156   fOutputList->Add(fV0PercentAll);
157   fZvtx    = new TH1F("fZvtx","",500,-20,20);
158   fOutputList->Add(fZvtx);
159   fV0Mult1d     = new TH1F("fV0Mult1d","Total v0 mult",500,0,25000);
160   fOutputList->Add(fV0Mult1d);
161
162   fV0Charge2d     = new TH2F("fV0Charge2d","V0C vs V0A charge",125,0,25000,125,0,25000);
163   fOutputList->Add(fV0Charge2d);
164   fV0Charge2dPercent     = new TH2F("fV0Charge2dPercent","V0C vs V0A charge vs centrality percentile",125,0,25000,125,0,25000);
165   fOutputList->Add(fV0Charge2dPercent);
166   fV0Charge2dAll  = new TH2F("fV0Charge2dAll","V0C vs V0A charge (all events)",125,0,25000,125,0,25000);
167   fOutputList->Add(fV0Charge2dAll);
168
169   fV0PercentBins = new TH1F*[fNThr];
170   fV0PercentBinsAll = new TH1F*[fNThr];
171   for(Int_t j = 0; j < fNThr; ++j) {
172     fV0PercentBins[j] = new TH1F(Form("fV0PercentBins_%d",j),"Centrality percentile with V0 charge thresholds",400,0,100);
173     fOutputList->Add(fV0PercentBins[j]);
174     fV0PercentBinsAll[j] = new TH1F(Form("fV0PercentBinsAll_%d",j),"Centrality percentile with V0 charge thresholds (no event selection)",100,0,100);
175     fOutputList->Add(fV0PercentBinsAll[j]);
176   }
177
178   fV0Cent = new TH1F("fV0Cent","Centrality percentile with custom V0 charge thresholds",400,0,100);
179   fOutputList->Add(fV0Cent);
180   fV0CentAll = new TH1F("fV0CentAll","Centrality percentile with custom V0 charge thresholds (no event selection)",100,0,100);
181   fOutputList->Add(fV0CentAll);
182   fV0SemiCent = new TH1F("fV0SemiCent","Centrality percentile with custom V0 charge thresholds",400,0,100);
183   fOutputList->Add(fV0SemiCent);
184   fV0SemiCentAll = new TH1F("fV0SemiCentAll","Centrality percentile with custom V0 charge thresholds (no event selection)",100,0,100);
185   fOutputList->Add(fV0SemiCentAll);
186
187   fV0CentHw = new TH1F("fV0CentHw","Centrality percentile with hardware V0 charge thresholds",400,0,100);
188   fOutputList->Add(fV0CentHw);
189   fV0CentHwAll = new TH1F("fV0CentHwAll","Centrality percentile with hardware V0 charge thresholds (no event selection)",100,0,100);
190   fOutputList->Add(fV0CentHwAll);
191   fV0SemiCentHw = new TH1F("fV0SemiCentHw","Centrality percentile with hardware V0 charge thresholds",400,0,100);
192   fOutputList->Add(fV0SemiCentHw);
193   fV0SemiCentHwAll = new TH1F("fV0SemiCentHwAll","Centrality percentile with hardware V0 charge thresholds (no event selection)",100,0,100);
194   fOutputList->Add(fV0SemiCentHwAll);
195
196   fV0CentTr = new TH1F("fV0CentTr","Centrality percentile with hardware V0 charge thresholds",400,0,100);
197   fOutputList->Add(fV0CentTr);
198   fV0CentTrAll = new TH1F("fV0CentTrAll","Centrality percentile with hardware V0 charge thresholds (no event selection)",100,0,100);
199   fOutputList->Add(fV0CentTrAll);
200   fV0SemiCentTr = new TH1F("fV0SemiCentTr","Centrality percentile with hardware V0 charge thresholds",400,0,100);
201   fOutputList->Add(fV0SemiCentTr);
202   fV0SemiCentTrAll = new TH1F("fV0SemiCentTrAll","Centrality percentile with hardware V0 charge thresholds (no event selection)",100,0,100);
203   fOutputList->Add(fV0SemiCentTrAll);
204
205   fV0Percent63    = new TH1F("fV0Percent63","Centrality percentile based on v0 mult",400,0,100);
206   fOutputList->Add(fV0Percent63);
207   fV0Percent63All    = new TH1F("fV0Percent63All","Centrality percentile based on v0 mult (no event selection)",400,0,100);
208   fOutputList->Add(fV0Percent63All);
209
210   fV0MultAll = new TH2F("fV0MultAll","",250,0,25000,250,0,25000);
211   fOutputList->Add(fV0MultAll);
212   fV0Mult63 = new TH2F("fV0Mult63","",250,0,25000,250,0,25000);
213   fOutputList->Add(fV0Mult63);
214
215   PostData(1, fOutputList);
216  }
217 //________________________________________________________________________
218 void AliAnaVZEROTrigger::Init()
219 {
220   // Nothing here
221   // ....
222 }
223
224 //________________________________________________________________________
225 void AliAnaVZEROTrigger::UserExec(Option_t *) 
226 {
227   // Main loop
228   // Called for each event
229
230
231   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
232   if (!fESD) {
233     printf("ERROR: fESD not available\n");
234     return;
235   }
236
237   AliESDVZERO* esdV0 = fESD->GetVZEROData();
238   if (!esdV0) {
239     Printf("ERROR: esd V0  not available");
240     return;
241   }
242
243   // Trigger
244   TString trigStr(fESD->GetFiredTriggerClasses());
245   if (!trigStr.Contains("-B-")) return;
246   if (!trigStr.Contains(fMBTrigName.Data())) return;
247
248   // Count V0 flags
249   Int_t nV0A = 0;
250   Int_t nV0C = 0;
251   for(Int_t i = 0; i < 32; ++i) {
252     if (esdV0->GetBBFlag(i)) nV0C++;
253     if (esdV0->GetBBFlag(i+32)) nV0A++;
254   }
255
256   // Phys sel
257   Bool_t goodEvent = kTRUE;
258   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);
259   if (!isSelected) goodEvent = kFALSE;
260
261   const AliESDVertex *primaryVtx = fESD->GetPrimaryVertex();
262   if (!primaryVtx) goodEvent = kFALSE;
263   if (!primaryVtx->GetStatus()) goodEvent = kFALSE;
264   Double_t tPrimaryVtxPosition[3];
265   primaryVtx->GetXYZ(tPrimaryVtxPosition);
266   if (goodEvent) fZvtx->Fill(tPrimaryVtxPosition[2]);
267   if (TMath::Abs(tPrimaryVtxPosition[2]) > 10.0) goodEvent = kFALSE;
268
269   AliCentrality *centrality = fESD->GetCentrality();
270   if (centrality->GetQuality()) goodEvent = kFALSE;
271   Float_t percentile = centrality->GetCentralityPercentile("V0M");
272   if (percentile < 0) percentile = 0;
273
274   fV0MultAll->Fill(esdV0->GetMTotV0A(),esdV0->GetMTotV0C());
275   if ((nV0A+nV0C)>=63)
276     fV0Mult63->Fill(esdV0->GetMTotV0A(),esdV0->GetMTotV0C());
277
278   fV0PercentAll->Fill(percentile);
279   if (goodEvent) fV0Percent->Fill(percentile);
280
281   if ((nV0A+nV0C)>=63) {
282     fV0Percent63All->Fill(percentile);
283     if (goodEvent) fV0Percent63->Fill(percentile);
284   }
285
286   Float_t multV0 = esdV0->GetMTotV0A()+esdV0->GetMTotV0C();
287   if (goodEvent) fV0Mult1d->Fill(multV0);
288
289   UShort_t chargeA = esdV0->GetTriggerChargeA();
290   UShort_t chargeC = esdV0->GetTriggerChargeC();
291   fV0Charge2dAll->Fill(chargeA,chargeC);
292   if (goodEvent) {
293     fV0Charge2d->Fill(chargeA,chargeC);
294     fV0Charge2dPercent->Fill(chargeA,chargeC,percentile);
295   }
296
297   for(Int_t j = 0; j < fNThr; ++j) {
298     Float_t thrA = GetThrA(j);
299     Float_t thrC = GetThrC(j);
300     if (chargeA >= thrA && chargeC >= thrC) {
301       fV0PercentBinsAll[j]->Fill(percentile);
302       if (goodEvent) fV0PercentBins[j]->Fill(percentile);
303     }
304   }
305
306   if (chargeA >= fCentCuts[0] && chargeC >= fCentCuts[1]) {
307     fV0CentAll->Fill(percentile);
308     if (goodEvent) fV0Cent->Fill(percentile);
309   }
310   if (chargeA >= fSemiCentCuts[0] && chargeC >= fSemiCentCuts[1]) {
311     fV0SemiCentAll->Fill(percentile);
312     if (goodEvent) fV0SemiCent->Fill(percentile);
313   }
314
315   if (esdV0->GetTriggerBits() & (1<<8)) {
316     fV0CentHwAll->Fill(percentile);
317     if (goodEvent) fV0CentHw->Fill(percentile);
318   }
319   if (esdV0->GetTriggerBits() & (1<<6)) {
320     fV0SemiCentHwAll->Fill(percentile);
321     if (goodEvent) fV0SemiCentHw->Fill(percentile);
322   }
323
324   if (trigStr.Contains("CVHN")) {
325     fV0CentHwAll->Fill(percentile);
326     if (goodEvent) fV0CentHw->Fill(percentile);
327   }
328   if (trigStr.Contains("CVLN")) {
329     fV0SemiCentHwAll->Fill(percentile);
330     if (goodEvent) fV0SemiCentHw->Fill(percentile);
331   }
332
333   PostData(1, fOutputList);
334 }      
335
336 //________________________________________________________________________
337 void AliAnaVZEROTrigger::Terminate(Option_t *) 
338 {
339   // Check the output list and store output config file
340   // Called once at the end of the query
341
342   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
343   if (!fOutputList) {
344     printf("ERROR: Output list not available\n");
345     return;
346   }
347
348   printf("\n\n");
349   FILE *fout=fopen("./trigger.txt","w");
350   if (!fout) {
351     printf("Failed to open local result file\n");
352     return;
353   }
354   printf("%.0f %.0f %f %d %.0f %.0f %.0f %.0f\n\n",
355          GetMinThr(),
356          GetMaxThr(),
357          GetRatio(),
358          GetNThr(),
359          GetSemiCentCuts(0),GetSemiCentCuts(1),
360          GetCentCuts(0),GetCentCuts(1));
361
362   fprintf(fout,"%.0f %.0f %f %d %.0f %.0f %.0f %.0f\n\n",
363           GetMinThr(),
364           GetMaxThr(),
365           GetRatio(),
366           GetNThr(),
367           GetSemiCentCuts(0),GetSemiCentCuts(1),
368           GetCentCuts(0),GetCentCuts(1));
369   fclose(fout);
370 }
371
372 //________________________________________________________________________
373 Float_t AliAnaVZEROTrigger::GetThrA(Int_t j) const
374 {
375   // Get the threshold on A side with
376   // index i
377   if (j<0 || j>= fNThr) return 0;
378
379   return (fMinThr + ((Float_t)j)*(fMaxThr-fMinThr)/((Float_t)fNThr-1.));
380 }
381
382 //________________________________________________________________________
383 Float_t AliAnaVZEROTrigger::GetThrC(Int_t j) const
384 {
385   // Get the threshold on C side with
386   // index i
387   if (j<0 || j>= fNThr) return 0;
388
389   Float_t thrA = (fMinThr + ((Float_t)j)*(fMaxThr-fMinThr)/((Float_t)fNThr-1.));
390   return (thrA*fRatio);
391 }
392