PWGJE/EMCALJetTasks/AliJetEmbeddingFromPYTHIATask.cxx

   1 // $Id: AliJetEmbeddingFromPYTHIATask.cxx $
   2 //
   3 // Jet embedding from PYTHIA task.
   4 //
   5 // Author: S.Aiola, C.Loizides
   6
   7 #include "AliJetEmbeddingFromPYTHIATask.h"
   8
   9 #include <TFile.h>
  10 #include <TMath.h>
  11 #include <TString.h>
  12 #include <TRandom.h>
  13 #include <TParameter.h>
  14 #include <TSystem.h>
  15 #include <TH1F.h>
  16
  17 #include "AliVEvent.h"
  18 #include "AliLog.h"
  19
  20 ClassImp(AliJetEmbeddingFromPYTHIATask)
  21
  22 //________________________________________________________________________
  23 AliJetEmbeddingFromPYTHIATask::AliJetEmbeddingFromPYTHIATask() :
  24   AliJetEmbeddingFromAODTask("AliJetEmbeddingFromPYTHIATask"),
  25   fPYTHIAPath(),
  26   fPtHardBinScaling(),
  27   fLHC11hAnchorRun(kTRUE),
  28   fAnchorRun(-1),
  29   fCurrentPtHardBin(-1),
  30   fPtHardBinParam(0),
  31   fHistPtHardBins(0)
  32 {
  33   // Default constructor.
  34   SetSuffix("PYTHIAEmbedding");
  35   fTotalFiles = 2000;
  36   fRandomAccess = kTRUE;
  37   SetMC(kTRUE);
  38 }
  39
  40 //________________________________________________________________________
  41 AliJetEmbeddingFromPYTHIATask::AliJetEmbeddingFromPYTHIATask(const char *name, Bool_t drawqa) :
  42   AliJetEmbeddingFromAODTask(name, drawqa),
  43   fPYTHIAPath("/alice/sim/2012/LHC12a15e"),
  44   fPtHardBinScaling(),
  45   fLHC11hAnchorRun(kTRUE),
  46   fAnchorRun(-1),
  47   fCurrentPtHardBin(-1),
  48   fPtHardBinParam(0),
  49   fHistPtHardBins(0)
  50 {
  51   // Standard constructor.
  52   SetSuffix("PYTHIAEmbedding");
  53   fTotalFiles = 2000;
  54   fRandomAccess = kTRUE;
  55   SetMC(kTRUE);
  56 }
  57
  58 //________________________________________________________________________
  59 AliJetEmbeddingFromPYTHIATask::~AliJetEmbeddingFromPYTHIATask()
  60 {
  61   // Destructor
  62 }
  63
  64 //________________________________________________________________________
  65 void AliJetEmbeddingFromPYTHIATask::UserCreateOutputObjects()
  66 {
  67   if (!fQAhistos)
  68     return;
  69
  70   AliJetModelBaseTask::UserCreateOutputObjects();
  71
  72   fHistPtHardBins = new TH1F("fHistPtHardBins", "fHistPtHardBins", 11, 0, 11);
  73   fHistPtHardBins->GetXaxis()->SetTitle("p_{T} hard bin");
  74   fHistPtHardBins->GetYaxis()->SetTitle("total events");
  75   fOutput->Add(fHistPtHardBins);
  76
  77   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
  78   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
  79
  80   for (Int_t i = 1; i < 12; i++)
  81     fHistPtHardBins->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
  82
  83   PostData(1, fOutput);
  84 }
  85
  86 //________________________________________________________________________
  87 Bool_t AliJetEmbeddingFromPYTHIATask::ExecOnce()
  88 {
  89   if (fPtHardBinScaling.GetSize() > 0) {
  90     Double_t sum = 0;
  91     for (Int_t i = 0; i < fPtHardBinScaling.GetSize(); i++)
  92       sum += fPtHardBinScaling[i];
  93
  94     if (sum == 0) {
  95       AliWarning("No hard pt bin scaling!");
  96       sum = fPtHardBinScaling.GetSize();
  97     }
  98
  99     for (Int_t i = 0; i < fPtHardBinScaling.GetSize(); i++)
 100       fPtHardBinScaling[i] /= sum;
 101   }
 102
 103   fPtHardBinParam = static_cast<TParameter<int>*>(InputEvent()->FindListObject("PYTHIAPtHardBin"));
 104   if (!fPtHardBinParam) {
 105     fPtHardBinParam = new TParameter<int>("PYTHIAPtHardBin", 0);
 106     AliDebug(2,"Adding pt hard bin param object to the event list...");
 107     InputEvent()->AddObject(fPtHardBinParam);
 108   }
 109
 110   return AliJetEmbeddingFromAODTask::ExecOnce();
 111 }
 112
 113 //________________________________________________________________________
 114 Bool_t AliJetEmbeddingFromPYTHIATask::GetNextEntry()
 115 {
 116   Int_t newPtHard = GetRandomPtHardBin();
 117
 118   new (fPtHardBinParam) TParameter<int>("PYTHIAPtHardBin", newPtHard);
 119
 120   if (fHistPtHardBins)
 121     fHistPtHardBins->SetBinContent(newPtHard+1, fHistPtHardBins->GetBinContent(newPtHard+1)+1);
 122
 123   if (newPtHard != fCurrentPtHardBin) {
 124     fCurrentPtHardBin = newPtHard;
 125     if (!OpenNextFile()) return kFALSE;
 126   }
 127
 128   return AliJetEmbeddingFromAODTask::GetNextEntry();
 129 }
 130
 131 //________________________________________________________________________
 132 Int_t AliJetEmbeddingFromPYTHIATask::GetRandomPtHardBin()
 133 {
 134   static Int_t order[20]={-1};
 135
 136   if (order[0] == -1)
 137     TMath::Sort(fPtHardBinScaling.GetSize(), fPtHardBinScaling.GetArray(), order);
 138
 139   Double_t rnd = gRandom->Rndm();
 140   Double_t sum = 0;
 141   Int_t ptHard = -1;
 142   for (Int_t i = 0; i < fPtHardBinScaling.GetSize(); i++) {
 143     sum += fPtHardBinScaling[order[i]];
 144     if (sum >= rnd) {
 145       ptHard = order[i];
 146       break;
 147     }
 148   }
 149
 150   return ptHard;
 151 }
 152
 153 //________________________________________________________________________
 154 Bool_t AliJetEmbeddingFromPYTHIATask::UserNotify()
 155 {
 156   if (!fLHC11hAnchorRun)
 157     return kTRUE;
 158
 159   Int_t runNumber = InputEvent()->GetRunNumber();
 160
 161   Int_t semiGoodRunList[28] = {169975, 169981, 170038, 170040, 170083, 170084, 170085, 170088,
 162                                170089, 170091, 170152, 170155, 170159, 170163, 170193, 170195,
 163                                170203, 170204, 170205, 170228, 170230, 170264, 170268, 170269,
 164                                170270, 170306, 170308, 170309};
 165
 166   fAnchorRun = 169838; // Assume it is a good run
 167
 168   for (Int_t i = 0; i < 28; i++) {
 169     if (runNumber == semiGoodRunList[i]) { // If it is semi good, change the anchor run
 170       fAnchorRun = 170040;
 171       break;
 172     }
 173   }
 174
 175   return kTRUE;
 176 }
 177
 178 //________________________________________________________________________
 179 TString AliJetEmbeddingFromPYTHIATask::GetNextFileName()
 180 {
 181   fCurrentAODFileID = TMath::Nint(gRandom->Rndm()*(fTotalFiles-1))+1;
 182
 183   TString fileName(fPYTHIAPath);
 184   if (!fileName.EndsWith("/"))
 185     fileName += "/";
 186
 187   fileName += fAnchorRun;
 188   fileName += "/";
 189   fileName += fCurrentPtHardBin;
 190   fileName += "/";
 191   if (fCurrentAODFileID < 10)
 192     fileName += "00";
 193   else if (fCurrentAODFileID < 100)
 194     fileName += "0";
 195   fileName += fCurrentAODFileID;
 196
 197   if (!gSystem->AccessPathName(Form("%s/AliAOD.root", fileName.Data())))
 198     fileName += "/AliAOD.root";
 199   else if (!gSystem->AccessPathName(Form("%s/aod_archive.zip", fileName.Data())))
 200     fileName += "/aod_archive.zip#AliAOD.root";
 201   else
 202     fileName += "/root_archive.zip#AliAOD.root";
 203
 204   return fileName;
 205 }