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