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       for (Int_t i = 0; i < fPtHardBinScaling.GetSize(); i++)
 120         fPtHardBinScaling[i] /= sum;
 121     }
 122     else {
 123       for (Int_t i = 0; i < fPtHardBinScaling.GetSize(); i++)
 124         fPtHardBinScaling[i] /= sum;
 125     }
 126   }
 127
 128   fPtHardBinParam = static_cast<TParameter<int>*>(InputEvent()->FindListObject("PYTHIAPtHardBin"));
 129   if (!fPtHardBinParam) {
 130     fPtHardBinParam = new TParameter<int>("PYTHIAPtHardBin", 0);
 131     AliDebug(3,"Adding pt hard bin param object to the event list...");
 132     InputEvent()->AddObject(fPtHardBinParam);
 133   }
 134
 135   return AliJetEmbeddingFromAODTask::ExecOnce();
 136 }
 137
 138 //________________________________________________________________________
 139 Bool_t AliJetEmbeddingFromPYTHIATask::GetNextEntry()
 140 {
 141   if (fCurrentPtHardBin < 0 || (fMinEntriesPerPtHardBin >= 0 && fPtHardBinCount >= fMinEntriesPerPtHardBin)) {
 142     fPtHardBinCount = 0;
 143
 144     Int_t newPtHard = GetRandomPtHardBin();
 145
 146     if (newPtHard != fCurrentPtHardBin) {
 147       fPtHardBinParam->SetVal(newPtHard);
 148       fCurrentPtHardBin = newPtHard;
 149       if (!OpenNextFile()) return kFALSE;
 150     }
 151   }
 152
 153   fPtHardBinCount++;
 154   if (fHistPtHardBins) fHistPtHardBins->SetBinContent(fCurrentPtHardBin+1, fHistPtHardBins->GetBinContent(fCurrentPtHardBin+1)+1);
 155
 156   return AliJetEmbeddingFromAODTask::GetNextEntry();
 157 }
 158
 159 //________________________________________________________________________
 160 Int_t AliJetEmbeddingFromPYTHIATask::GetRandomPtHardBin()
 161 {
 162   static Int_t order[20]={-1};
 163
 164   if (order[0] == -1)
 165     TMath::Sort(fPtHardBinScaling.GetSize(), fPtHardBinScaling.GetArray(), order);
 166
 167   Double_t rnd = gRandom->Rndm();
 168   Double_t sum = 0;
 169   Int_t ptHard = -1;
 170   for (Int_t i = 0; i < fPtHardBinScaling.GetSize(); i++) {
 171     sum += fPtHardBinScaling[order[i]];
 172     if (sum >= rnd) {
 173       ptHard = order[i];
 174       break;
 175     }
 176   }
 177
 178   return ptHard;
 179 }
 180
 181 //________________________________________________________________________
 182 Bool_t AliJetEmbeddingFromPYTHIATask::UserNotify()
 183 {
 184   if (!fLHC11hAnchorRun)
 185     return kTRUE;
 186
 187   Int_t runNumber = InputEvent()->GetRunNumber();
 188
 189   Int_t semiGoodRunList[32] = {169975, 169981, 170038, 170040, 170083, 170084, 170085,
 190                                170088, 170089, 170091, 170152, 170155, 170159, 170163,
 191                                170193, 170195, 170203, 170204, 170228, 170230, 170268,
 192                                170269, 170270, 170306, 170308, 170309, 169238, 169160,
 193                                169156, 169148, 169145, 169144};
 194
 195   fAnchorRun = 169838; // Assume it is a good run
 196
 197   for (Int_t i = 0; i < 32; i++) {
 198     if (runNumber == semiGoodRunList[i]) { // If it is semi good, change the anchor run
 199       fAnchorRun = 170040;
 200       break;
 201     }
 202   }
 203
 204   return kTRUE;
 205 }
 206
 207 //________________________________________________________________________
 208 TFile* AliJetEmbeddingFromPYTHIATask::GetNextFile()
 209 {
 210   fCurrentAODFileID = TMath::Nint(gRandom->Rndm()*(fTotalFiles-1))+1;
 211
 212   if (fMinEntriesPerPtHardBin < 0) {
 213     fCurrentPtHardBin = GetRandomPtHardBin();
 214     fPtHardBinParam->SetVal(fCurrentPtHardBin);
 215   }
 216
 217   TString fileName;
 218
 219   if (fAnchorRun>0)
 220     fileName = Form(fPYTHIAPath.Data(), fAnchorRun, fCurrentPtHardBin, fCurrentAODFileID);
 221   else
 222     fileName = Form(fPYTHIAPath.Data(), fCurrentPtHardBin, fCurrentAODFileID);
 223
 224   if (fFileTable && fFileTable->GetEntries() > 0) {
 225     TObject* obj = fFileTable->FindObject(fileName);
 226     if (obj != 0 && fUseAsVetoTable) {
 227       AliWarning(Form("File %s found in the vetoed file table. Skipping...", fileName.Data()));
 228       return 0;
 229     }
 230     if (obj == 0 && !fUseAsVetoTable) {
 231       AliWarning(Form("File %s not found in the allowed file table. Skipping...", fileName.Data()));
 232       return 0;
 233     }
 234   }
 235
 236   if (fileName.BeginsWith("alien://") && !gGrid) {
 237     AliInfo("Trying to connect to AliEn ...");
 238     TGrid::Connect("alien://");
 239   }
 240
 241   TString baseFileName(fileName);
 242   if (baseFileName.Contains(".zip#")) {
 243     Ssiz_t pos = baseFileName.Last('#');
 244     baseFileName.Remove(pos);
 245   }
 246
 247   if (gSystem->AccessPathName(baseFileName)) {
 248     AliError(Form("File %s does not exist!", baseFileName.Data()));
 249     return 0;
 250   }
 251
 252   AliDebug(3,Form("Trying to open file %s...", fileName.Data()));
 253   TFile *file = TFile::Open(fileName);
 254
 255   if (!file || file->IsZombie()) {
 256     AliError(Form("Unable to open file: %s!", fileName.Data()));
 257     return 0;
 258   }
 259
 260   return file;
 261 }