setter to assume pion mass for clusters
[u/mrichter/AliRoot.git] / 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 || fCurrentPtHardBin < 0) {
139     fPtHardBinCount = 0;
140
141     Int_t newPtHard = GetRandomPtHardBin();
142     
143     if (newPtHard != fCurrentPtHardBin) {
144       new (fPtHardBinParam) TParameter<int>("PYTHIAPtHardBin", newPtHard);
145       fCurrentPtHardBin = newPtHard;
146       if (!OpenNextFile()) return kFALSE;
147     }
148   }
149
150   fPtHardBinCount++;
151   fHistPtHardBins->SetBinContent(fCurrentPtHardBin+1, fHistPtHardBins->GetBinContent(fCurrentPtHardBin+1)+1);
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 }