PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalHighMultTrigger.cxx

   1 // $Id$
   2 //
   3 // High multiplicity pp trigger analysis task.
   4 //
   5 // Author: M. Verweij
   6
   7 #include <TClonesArray.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13
  14 #include "AliLog.h"
  15 #include "AliEmcalTriggerPatchInfo.h"
  16 #include "AliParticleContainer.h"
  17 #include "AliVVZERO.h"
  18
  19 #include "AliAnalysisTaskEmcalHighMultTrigger.h"
  20
  21 ClassImp(AliAnalysisTaskEmcalHighMultTrigger)
  22
  23 //________________________________________________________________________
  24 AliAnalysisTaskEmcalHighMultTrigger::AliAnalysisTaskEmcalHighMultTrigger() :
  25   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalHighMultTrigger", kTRUE),
  26   fNExLP(1),
  27   fNAccPatches(-1),
  28   fMedianEnergy(0),
  29   fMedianEnergyExLP(0),
  30   fSumEnergy(0),
  31   fSumEnergyExLP(0),
  32   fTruncatedMean(0),
  33   fTruncateThreshold(6.),
  34   fHistPatchEtaPhiE(0),
  35   fHistEnergyMedian(0),
  36   fHistEnergyMedianExLP(0),
  37   fHistEnergySum(0),
  38   fHistEnergySumExLP(0),
  39   fHistTruncatedMean(0),
  40   fHistTracks(0),
  41   fHistTracklets(0),
  42   fHistV0MultSum(0),
  43   fHistTracksTracklets(0),
  44   fHistTracksV0MultSum(0),
  45   fHistSPDTrkClsSum(0),
  46   fHistSPDTrkClsSumExLP(0),
  47   fHistSPDTrkClsMedian(0),
  48   fHistSPDTrkClsMedianExLP(0),
  49   fHistSPDTrkClsTruncMean(0)
  50 {
  51   // Default constructor.
  52
  53   const Int_t nMultEst = 3;
  54   for(Int_t i = 0; i<nMultEst; i++) {
  55     fHistEnergyMedianEst[i] = 0;
  56     fHistEnergyMedianExLPEst[i] = 0;
  57     fHistEnergySumEst[i] = 0;
  58     fHistEnergySumExLPEst[i] = 0;
  59     fHistEnergySumAvgEst[i] = 0;
  60     fHistEnergySumAvgExLPEst[i] = 0;
  61     fHistTruncatedMeanEst[i] = 0;
  62   }
  63
  64   SetMakeGeneralHistograms(kTRUE);
  65 }
  66
  67 //________________________________________________________________________
  68 AliAnalysisTaskEmcalHighMultTrigger::AliAnalysisTaskEmcalHighMultTrigger(const char *name) :
  69   AliAnalysisTaskEmcalJet(name, kTRUE),
  70   fNExLP(1),
  71   fNAccPatches(-1),
  72   fMedianEnergy(0),
  73   fMedianEnergyExLP(0),
  74   fSumEnergy(0),
  75   fSumEnergyExLP(0),
  76   fTruncatedMean(0),
  77   fTruncateThreshold(6.),
  78   fHistPatchEtaPhiE(0),
  79   fHistEnergyMedian(0),
  80   fHistEnergyMedianExLP(0),
  81   fHistEnergySum(0),
  82   fHistEnergySumExLP(0),
  83   fHistTruncatedMean(0),
  84   fHistTracks(0),
  85   fHistTracklets(0),
  86   fHistV0MultSum(0),
  87   fHistTracksTracklets(0),
  88   fHistTracksV0MultSum(0),
  89   fHistSPDTrkClsSum(0),
  90   fHistSPDTrkClsSumExLP(0),
  91   fHistSPDTrkClsMedian(0),
  92   fHistSPDTrkClsMedianExLP(0),
  93   fHistSPDTrkClsTruncMean(0)
  94 {
  95   // Standard constructor.
  96
  97   const Int_t nMultEst = 3;
  98   for(Int_t i = 0; i<nMultEst; i++) {
  99     fHistEnergyMedianEst[i] = 0;
 100     fHistEnergyMedianExLPEst[i] = 0;
 101     fHistEnergySumEst[i] = 0;
 102     fHistEnergySumExLPEst[i] = 0;
 103     fHistEnergySumAvgEst[i] = 0;
 104     fHistEnergySumAvgExLPEst[i] = 0;
 105     fHistTruncatedMeanEst[i] = 0;
 106   }
 107
 108   SetMakeGeneralHistograms(kTRUE);
 109 }
 110
 111 //________________________________________________________________________
 112 AliAnalysisTaskEmcalHighMultTrigger::~AliAnalysisTaskEmcalHighMultTrigger()
 113 {
 114   // Destructor.
 115 }
 116
 117 //________________________________________________________________________
 118 void AliAnalysisTaskEmcalHighMultTrigger::UserCreateOutputObjects()
 119 {
 120   // Create user output.
 121
 122   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 123
 124   Int_t nE = 500;
 125   Double_t minE = 0.;
 126   Double_t maxE = 50.;
 127
 128   TString histName = "";
 129   TString histTitle = "";
 130
 131   histName = Form("fHistPatchEtaPhiE");
 132   histTitle = Form("%s;#eta;#phi;E",histName.Data());
 133   fHistPatchEtaPhiE = new TH3F(histName.Data(),histTitle.Data(),100,-1.,1.,18*8,0.,TMath::TwoPi(),100,0.,100.);
 134   fOutput->Add(fHistPatchEtaPhiE);
 135
 136   histName = Form("fHistEnergyMedian");
 137   histTitle = Form("%s;med[#it{E}]",histName.Data());
 138   fHistEnergyMedian = new TH1F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.);
 139   fOutput->Add(fHistEnergyMedian);
 140
 141   histName = Form("fHistEnergyMedianExLP");
 142   histTitle = Form("%s;med[#it{E}]",histName.Data());
 143   fHistEnergyMedianExLP = new TH1F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.);
 144   fOutput->Add(fHistEnergyMedianExLP);
 145
 146   histName = Form("fHistEnergySum");
 147   histTitle = Form("%s;#sum[#it{E}]",histName.Data());
 148   fHistEnergySum = new TH1F(histName.Data(),histTitle.Data(),nE,minE,maxE);
 149   fOutput->Add(fHistEnergySum);
 150
 151   histName = Form("fHistEnergySumExLP");
 152   histTitle = Form("%s;#sum[#it{E}]",histName.Data());
 153   fHistEnergySumExLP = new TH1F(histName.Data(),histTitle.Data(),nE,minE,maxE);
 154   fOutput->Add(fHistEnergySumExLP);
 155
 156   histName = Form("fHistTruncatedMean");
 157   histTitle = Form("%s;#sum[#it{E}]",histName.Data());
 158   fHistTruncatedMean = new TH1F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.);
 159   fOutput->Add(fHistTruncatedMean);
 160
 161   histName = Form("fHistTracks");
 162   histTitle = Form("%s;#it{N}_{tracks}",histName.Data());
 163   fHistTracks = new TH1F(histName.Data(),histTitle.Data(),300,0.,300.);
 164   fOutput->Add(fHistTracks);
 165
 166   histName = Form("fHistTracklets");
 167   histTitle = Form("%s;#it{N}_{tracklets}",histName.Data());
 168   fHistTracklets = new TH1F(histName.Data(),histTitle.Data(),300,0.,300.);
 169   fOutput->Add(fHistTracklets);
 170
 171   histName = Form("fHistV0MultSum");
 172   histTitle = Form("%s;mult[V0A+V0C]",histName.Data());
 173   fHistV0MultSum = new TH1F(histName.Data(),histTitle.Data(),500,0.,500.);
 174   fOutput->Add(fHistV0MultSum);
 175
 176   const Int_t nMultEst = 3;
 177   Int_t nBinsMultEst[nMultEst] = {300,200,500};
 178   Double_t multEstMax[nMultEst] = {300.,200.,500.};
 179   TString strMultEst[nMultEst] = {"Tracks","Tracklets","V0MultSum"};
 180   for(Int_t i = 0; i<nMultEst; i++) {
 181     histName = Form("fHistEnergyMedianEst%s",strMultEst[i].Data());
 182     histTitle = Form("%s;med[#it{E}];%s",histName.Data(),strMultEst[i].Data());
 183     fHistEnergyMedianEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.,nBinsMultEst[i],0.,multEstMax[i]);
 184     fOutput->Add(fHistEnergyMedianEst[i]);
 185
 186     histName = Form("fHistEnergyMedianExLPEst%s",strMultEst[i].Data());
 187     histTitle = Form("%s;med[#it{E}];%s",histName.Data(),strMultEst[i].Data());
 188     fHistEnergyMedianExLPEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.,nBinsMultEst[i],0.,multEstMax[i]);
 189     fOutput->Add(fHistEnergyMedianExLPEst[i]);
 190
 191     histName = Form("fHistEnergySumEst%s",strMultEst[i].Data());
 192     histTitle = Form("%s;#sum[#it{E}];%s",histName.Data(),strMultEst[i].Data());
 193     fHistEnergySumEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE,nBinsMultEst[i],0.,multEstMax[i]);
 194     fOutput->Add(fHistEnergySumEst[i]);
 195
 196     histName = Form("fHistEnergySumExLPEst%s",strMultEst[i].Data());
 197     histTitle = Form("%s;#sum[#it{E}];%s",histName.Data(),strMultEst[i].Data());
 198     fHistEnergySumExLPEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE,nBinsMultEst[i],0.,multEstMax[i]);
 199     fOutput->Add(fHistEnergySumExLPEst[i]);
 200
 201     histName = Form("fHistEnergySumAvgEst%s",strMultEst[i].Data());
 202     histTitle = Form("%s;#sum[#it{E}];%s",histName.Data(),strMultEst[i].Data());
 203     fHistEnergySumAvgEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.,nBinsMultEst[i],0.,multEstMax[i]);
 204     fOutput->Add(fHistEnergySumAvgEst[i]);
 205
 206     histName = Form("fHistEnergySumAvgExLPEst%s",strMultEst[i].Data());
 207     histTitle = Form("%s;#sum[#it{E}];%s",histName.Data(),strMultEst[i].Data());
 208     fHistEnergySumAvgExLPEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.,nBinsMultEst[i],0.,multEstMax[i]);
 209     fOutput->Add(fHistEnergySumAvgExLPEst[i]);
 210
 211     histName = Form("fHistTruncatedMeanEst%s",strMultEst[i].Data());
 212     histTitle = Form("%s;#LT#it{E}#GT_{trunc};%s",histName.Data(),strMultEst[i].Data());
 213     fHistTruncatedMeanEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.,nBinsMultEst[i],0.,multEstMax[i]);
 214     fOutput->Add(fHistTruncatedMeanEst[i]);
 215   }
 216
 217   histName = Form("fHistTracksTracklets");
 218   histTitle = Form("%s;%s;%s",histName.Data(),strMultEst[0].Data(),strMultEst[1].Data());
 219   fHistTracksTracklets = new TH2F(histName.Data(),histTitle.Data(),nBinsMultEst[0],0.,multEstMax[0],nBinsMultEst[1],0.,multEstMax[1]);
 220   fOutput->Add(fHistTracksTracklets);
 221
 222   histName = Form("fHistTracksV0MultSum");
 223   histTitle = Form("%s;%s;%s",histName.Data(),strMultEst[0].Data(),strMultEst[2].Data());
 224   fHistTracksV0MultSum = new TH2F(histName.Data(),histTitle.Data(),nBinsMultEst[0],0.,multEstMax[0],nBinsMultEst[2],0.,multEstMax[2]);
 225   fOutput->Add(fHistTracksV0MultSum);
 226
 227   histName = Form("fHistSPDTrkClsSum");
 228   histTitle = Form("%s;#it{N}_{tracklets,SPD};#it{N}_{clusters,SPD};#sum[#it{E}]",histName.Data());
 229   fHistSPDTrkClsSum = new TH3F(histName.Data(),histTitle.Data(),nBinsMultEst[1],0.,multEstMax[1],200,0.,1000.,nE,minE,maxE);
 230   fOutput->Add(fHistSPDTrkClsSum);
 231
 232   histName = Form("fHistSPDTrkClsSumExLP");
 233   histTitle = Form("%s;#it{N}_{tracklets,SPD};#it{N}_{clusters,SPD};#sum[#it{E}]",histName.Data());
 234   fHistSPDTrkClsSumExLP = new TH3F(histName.Data(),histTitle.Data(),nBinsMultEst[1],0.,multEstMax[1],200,0.,1000.,nE,minE,maxE);
 235   fOutput->Add(fHistSPDTrkClsSumExLP);
 236
 237   histName = Form("fHistSPDTrkClsMedian");
 238   histTitle = Form("%s;#it{N}_{tracklets,SPD};#it{N}_{clusters,SPD};med[#it{E}]",histName.Data());
 239   fHistSPDTrkClsMedian = new TH3F(histName.Data(),histTitle.Data(),nBinsMultEst[1],0.,multEstMax[1],200,0.,1000.,nE,minE,maxE/4.);
 240   fOutput->Add(fHistSPDTrkClsMedian);
 241
 242   histName = Form("fHistSPDTrkClsMedianExLP");
 243   histTitle = Form("%s;#it{N}_{tracklets,SPD};#it{N}_{clusters,SPD};med[#it{E}]",histName.Data());
 244   fHistSPDTrkClsMedianExLP = new TH3F(histName.Data(),histTitle.Data(),nBinsMultEst[1],0.,multEstMax[1],200,0.,1000.,nE,minE,maxE/4.);
 245   fOutput->Add(fHistSPDTrkClsMedianExLP);
 246
 247   histName = Form("fHistSPDTrkClsTruncMean");
 248   histTitle = Form("%s;#it{N}_{tracklets,SPD};#it{N}_{clusters,SPD};#LT#it{E}#GT_{trunc}",histName.Data());
 249   fHistSPDTrkClsTruncMean = new TH3F(histName.Data(),histTitle.Data(),nBinsMultEst[1],0.,multEstMax[1],200,0.,1000.,nE,minE,maxE/4.);
 250   fOutput->Add(fHistSPDTrkClsTruncMean);
 251
 252   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 253 }
 254
 255 //________________________________________________________________________
 256 Bool_t AliAnalysisTaskEmcalHighMultTrigger::FillHistograms()
 257 {
 258   // Fill histograms.
 259
 260   fHistEnergyMedian->Fill(fMedianEnergy);
 261   fHistEnergyMedianExLP->Fill(fMedianEnergyExLP);
 262   fHistEnergySum->Fill(fSumEnergy);
 263   fHistEnergySumExLP->Fill(fSumEnergyExLP);
 264   fHistTruncatedMean->Fill(fTruncatedMean);
 265
 266   //Multiplicity estimators
 267   Int_t nTracks   = -1;
 268   if (GetParticleContainer(0)) nTracks = GetParticleContainer(0)->GetNAcceptedParticles();
 269   Int_t nTracklets = InputEvent()->GetMultiplicity()->GetNumberOfTracklets();
 270
 271   AliVVZERO* vV0 = InputEvent()->GetVZEROData();
 272   Float_t multV0A=vV0->GetMTotV0A();
 273   Float_t multV0C=vV0->GetMTotV0C();
 274
 275   fHistTracks->Fill(nTracks);
 276   fHistTracklets->Fill(nTracklets);
 277   fHistV0MultSum->Fill(multV0A+multV0C);
 278
 279   const Int_t nMultEst = 3;
 280   Float_t multEst[nMultEst] = {(Float_t)(nTracks),(Float_t)(nTracklets),(Float_t)(multV0A+multV0C)};
 281   for(Int_t i = 0; i<nMultEst; i++) {
 282     fHistEnergyMedianEst[i]->Fill(fMedianEnergy,multEst[i]);
 283     fHistEnergyMedianExLPEst[i]->Fill(fMedianEnergyExLP,multEst[i]);
 284     fHistEnergySumEst[i]->Fill(fSumEnergy,multEst[i]);
 285     fHistEnergySumExLPEst[i]->Fill(fSumEnergyExLP,multEst[i]);
 286     if(fNAccPatches>0)
 287       fHistEnergySumAvgEst[i]->Fill(fSumEnergy/((Double_t)fNAccPatches),multEst[i]);
 288     if(fNAccPatches>1)
 289       fHistEnergySumAvgExLPEst[i]->Fill(fSumEnergyExLP/((Double_t)fNAccPatches-1.),multEst[i]);
 290     fHistTruncatedMeanEst[i]->Fill(fTruncatedMean,multEst[i]);
 291   }
 292
 293   fHistTracksTracklets->Fill(multEst[0],multEst[1]);
 294   fHistTracksV0MultSum->Fill(multEst[0],multEst[2]);
 295
 296   Int_t nClustersLayer0 = InputEvent()->GetNumberOfITSClusters(0);
 297   Int_t nClustersLayer1 = InputEvent()->GetNumberOfITSClusters(1);
 298   fHistSPDTrkClsSum->Fill(multEst[1],(Float_t)(nClustersLayer0+nClustersLayer1),fSumEnergy);
 299   fHistSPDTrkClsSumExLP->Fill(multEst[1],(Float_t)(nClustersLayer0+nClustersLayer1),fSumEnergyExLP);
 300   fHistSPDTrkClsMedian->Fill(multEst[1],(Float_t)(nClustersLayer0+nClustersLayer1),fMedianEnergy);
 301   fHistSPDTrkClsMedianExLP->Fill(multEst[1],(Float_t)(nClustersLayer0+nClustersLayer1),fMedianEnergyExLP);
 302   fHistSPDTrkClsTruncMean->Fill(multEst[1],(Float_t)(nClustersLayer0+nClustersLayer1),fTruncatedMean);
 303
 304   return kTRUE;
 305 }
 306
 307 //________________________________________________________________________
 308 void AliAnalysisTaskEmcalHighMultTrigger::ExecOnce() {
 309
 310   AliAnalysisTaskEmcalJet::ExecOnce();
 311
 312 }
 313
 314 //________________________________________________________________________
 315 Bool_t AliAnalysisTaskEmcalHighMultTrigger::Run()
 316 {
 317   // Run analysis code here, if needed. It will be executed before FillHistograms().
 318
 319   if(!fTriggerPatchInfo) {
 320     AliDebug(11,Form("Couldn't find patch info object %s",fCaloTriggerPatchInfoName.Data()));
 321     return kFALSE;
 322   }
 323
 324   Int_t nPatch = fTriggerPatchInfo->GetEntriesFast();
 325
 326   //sort patches
 327   Double_t ptarr[999] = {0};
 328   Int_t indexes[999] = {0};
 329   Int_t iacc = 0;
 330   for(Int_t i = 0; i<nPatch; i++) {
 331     AliEmcalTriggerPatchInfo *patch = dynamic_cast<AliEmcalTriggerPatchInfo*>(fTriggerPatchInfo->At(i));
 332     if(!patch) continue;
 333     if(patch->GetPatchE()>0.) {
 334       ptarr[iacc] = patch->GetPatchE();
 335       iacc++;
 336     }
 337     fHistPatchEtaPhiE->Fill(patch->GetEtaMin(),patch->GetPhiMin(),patch->GetPatchE());
 338   }
 339
 340   TMath::Sort(nPatch,ptarr,indexes);
 341   Double_t ptarrSort[999];
 342   // Double_t indexesSort[999];
 343   Double_t ptarrSortExLP[999];
 344   //  Double_t indexesSortExLP[999];
 345   for(Int_t i = 0; i<iacc; i++) {
 346     ptarrSort[i] = ptarr[indexes[i]];
 347     //  indexesSort[i] = indexes[i];
 348     if(i>=fNExLP) {
 349       ptarrSortExLP[i] = ptarr[indexes[i]];
 350       //  indexesSortExLP[i] = indexes[i];
 351     }
 352   }
 353
 354   //Calculate sum energy
 355   fSumEnergy = 0.;
 356   fSumEnergyExLP = 0.;
 357   for(Int_t i = 0; i<iacc; i++) {
 358     fSumEnergy+= ptarr[indexes[i]];
 359     if(i>=fNExLP)
 360       fSumEnergyExLP+= ptarr[indexes[i]];
 361   }
 362
 363   //calculate median of all and excluding leading patch(es)
 364   fMedianEnergy = TMath::Median(iacc,ptarrSort);
 365   fMedianEnergyExLP = TMath::Median(iacc-fNExLP,ptarrSortExLP);
 366   fNAccPatches = iacc;
 367
 368   //calculate truncated mean
 369   Double_t it = 0.;
 370   Double_t sum = 0.;
 371   for(Int_t i = 0; i<iacc; i++) {
 372     if(ptarr[indexes[i]]<fTruncateThreshold) {
 373       sum+= ptarr[indexes[i]];
 374       it+=1.;
 375     }
 376   }
 377   if(it>0)
 378     fTruncatedMean = sum/it;
 379
 380    return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 381 }
 382
 383 //________________________________________________________________________
 384 void AliAnalysisTaskEmcalHighMultTrigger::Terminate(Option_t *)
 385 {
 386   // Called once at the end of the analysis.
 387 }
 388