Next version of the analysis task used for extraction of centrality trigger thr.
[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(400.),
27   fMaxThr(14000.),
28   fRatio(2.0),
29   fNThr(30),
30   fMBTrigName("CPBI"),
31   fUsePhysSel(kFALSE),
32   fV0Percent(0),
33   fV0PercentAll(0),
34   fZvtx(0),
35   fXYvtx(0),
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),
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)
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 //________________________________________________________________________
85 AliAnaVZEROTrigger::AliAnaVZEROTrigger(const char *name) 
86   : AliAnalysisTaskSE(name), fESD(0), fOutputList(0),
87   fMinThr(400.),
88   fMaxThr(14000.),
89   fRatio(2.0),
90   fNThr(30),
91   fMBTrigName("CPBI"),
92   fUsePhysSel(kFALSE),
93   fV0Percent(0),
94   fV0PercentAll(0),
95   fZvtx(0),
96   fXYvtx(0),
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),
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)
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 //________________________________________________________________________
146 void 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 //________________________________________________________________________
169 void 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);
183   fXYvtx    = new TH2F("fXYvtx","",250,-0.5,0.5,250,-0.5,0.5);
184   fOutputList->Add(fXYvtx);
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
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
263   PostData(1, fOutputList);
264  }
265 //________________________________________________________________________
266 void AliAnaVZEROTrigger::Init()
267 {
268   // Nothing here
269   // ....
270 }
271
272 //________________________________________________________________________
273 void 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;
306   Bool_t isSelected;
307   if (fUsePhysSel)
308     isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);
309   else
310     isSelected = ((esdV0->GetV0ADecision()==1) && (esdV0->GetV0CDecision()==1));
311
312   if (!isSelected) goodEvent = kFALSE;
313
314   const AliESDVertex *primaryVtx = fESD->GetPrimaryVertex();
315   if (!primaryVtx) goodEvent = kFALSE;
316   if (!primaryVtx->GetStatus()) goodEvent = kFALSE;
317   Double_t tPrimaryVtxPosition[3];
318   primaryVtx->GetXYZ(tPrimaryVtxPosition);
319   if (goodEvent) {
320     fZvtx->Fill(tPrimaryVtxPosition[2]);
321     fXYvtx->Fill(tPrimaryVtxPosition[0],tPrimaryVtxPosition[1]);
322   }
323   if (TMath::Abs(tPrimaryVtxPosition[2]) > 10.0) goodEvent = kFALSE;
324
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;
351   AliCentrality *centrality = fESD->GetCentrality();
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   }
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
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")) {
411     fV0CentTrAll->Fill(percentile);
412     if (goodEvent) fV0CentTr->Fill(percentile);
413   }
414   if (trigStr.Contains("CVLN")) {
415     fV0SemiCentTrAll->Fill(percentile);
416     if (goodEvent) fV0SemiCentTr->Fill(percentile);
417   }
418
419   PostData(1, fOutputList);
420 }      
421
422 //________________________________________________________________________
423 void 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 //________________________________________________________________________
459 Float_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 //________________________________________________________________________
469 Float_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