PWGPP/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::kMB);
 309   else
 310     isSelected = ((esdV0->GetV0ADecision()==1) && (esdV0->GetV0CDecision()==1));
 311
 312   if (!isSelected) goodEvent = kFALSE;
 313
 314   const AliESDVertex *primaryVtx = fESD->GetPrimaryVertexSPD();
 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