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