PWGJE/EMCALJetTasks/AliJetEmbeddingFromPYTHIATask.cxx

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