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