AliPhysics  5eaf189 (5eaf189)
AliAodSkimTask.cxx
Go to the documentation of this file.
1 #include <Riostream.h>
2 #include <TChain.h>
3 #include <TFile.h>
4 #include <TGraph.h>
5 #include <TH1F.h>
6 #include <TKey.h>
7 #include <TList.h>
8 #include <TMath.h>
9 #include <TProfile.h>
10 #include <TSystem.h>
11 #include <AliAODEvent.h>
12 #include <AliAODHandler.h>
13 #include <AliAODInputHandler.h>
14 #include <AliAODMCHeader.h>
15 #include <AliAODMCParticle.h>
16 #include <AliAnalysisManager.h>
17 #include <AliLog.h>
18 #include "AliAodSkimTask.h"
19 
20 using namespace std;
21 ClassImp(AliAodSkimTask)
22 
24  AliAnalysisTaskSE(), fClusMinE(-1), fCutMC(1), fYCutMC(0.7),
25  fDoCopyHeader(1), fDoCopyVZERO(1), fDoCopyTZERO(1), fDoCopyVertices(1), fDoCopyTOF(1), fDoCopyTracks(1), fDoCopyTrigger(1), fDoCopyPTrigger(0),
26  fDoCopyCells(1), fDoCopyPCells(0), fDoCopyClusters(1), fDoCopyDiMuons(0), fDoCopyZDC(1), fDoCopyMC(1), fDoCopyMCHeader(1), fTrials(0), fPyxsec(0),
27  fPytrials(0), fPypthardbin(0), fAOD(0), fAODMcHeader(0), fOutputList(0), fHevs(0), fHclus(0)
28 {
29 }
30 
31 AliAodSkimTask::AliAodSkimTask(const char* name) :
32  AliAnalysisTaskSE(name), fClusMinE(-1), fCutMC(1), fYCutMC(0.7),
33  fDoCopyHeader(1), fDoCopyVZERO(1), fDoCopyTZERO(1), fDoCopyVertices(1), fDoCopyTOF(1), fDoCopyTracks(1), fDoCopyTrigger(1), fDoCopyPTrigger(0),
34  fDoCopyCells(1), fDoCopyPCells(0), fDoCopyClusters(1), fDoCopyDiMuons(0), fDoCopyZDC(1), fDoCopyMC(1), fDoCopyMCHeader(1), fTrials(0), fPyxsec(0),
35  fPytrials(0), fPypthardbin(0), fAOD(0), fAODMcHeader(0), fOutputList(0), fHevs(0), fHclus(0)
36 {
37  DefineInput(0, TChain::Class());
38  DefineOutput(1, TList::Class());
39 }
40 
42 {
43  if (fOutputList) {
44  delete fOutputList;
45  }
46  delete fHevs;
47  delete fHclus;
48 }
49 
51 {
52 
53  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
54  AliAODHandler *oh = (AliAODHandler*)man->GetOutputEventHandler();
55  if (oh) {
56  TFile *fout = oh->GetTree()->GetCurrentFile();
57  fout->SetCompressionLevel(2);
58  }
59 
60  fOutputList = new TList;
61  fOutputList->SetOwner();
62  fHevs = new TH1F("hEvs","",2,-0.5,1.5);
63  fOutputList->Add(fHevs);
64  fHclus = new TH1F("hClus",";E (GeV)",200,0,100);
65  fOutputList->Add(fHclus);
66  PostData(1, fOutputList);
67 }
68 
70 {
71  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
72  if (!fAOD)
73  return;
74 
75  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
76  AliAODHandler *oh = (AliAODHandler*)man->GetOutputEventHandler();
77  if (!oh) {
78  AliFatal(Form("%s: No output handler found", GetName()));
79  return;
80  }
81  oh->SetFillAOD(kFALSE);
82 
83  Bool_t store = kFALSE;
84  if (fClusMinE>0) {
85  TClonesArray *cls = fAOD->GetCaloClusters();
86  for (Int_t i=0; i<cls->GetEntriesFast(); ++i) {
87  AliAODCaloCluster *clus = static_cast<AliAODCaloCluster*>(cls->At(i));
88  if (!clus->IsEMCAL())
89  continue;
90  Double_t e = clus->E();
91  fHclus->Fill(e);
92  if (e>fClusMinE) {
93  store = kTRUE;
94  }
95  }
96  } else {
97  store = kTRUE;
98  }
99 
100  if (!store) {
101  ++fTrials;
102  fHevs->Fill(0);
103  return;
104  }
105 
106  fHevs->Fill(1);
107 
108  oh->SetFillAOD(kTRUE);
109  AliAODEvent *eout = dynamic_cast<AliAODEvent*>(oh->GetAOD());
110  AliAODEvent *evin = dynamic_cast<AliAODEvent*>(InputEvent());
111  TTree *tout = oh->GetTree();
112  if (tout) {
113  TList *lout = tout->GetUserInfo();
114  if (lout->FindObject("alirootVersion")==0) {
115  TList *lin = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetUserInfo();
116  TString apver(gSystem->BaseName(gSystem->Getenv("ALICE_PHYSICS")));
117  lout->Add(new TObjString(Form("AodSkim: ver %s, tag %s with settings %s",GetVersion(),apver.Data(),Str())));
118  AliInfo(Form("%s: Set user info %s", GetName(), lout->At(0)->GetName()));
119  for (Int_t jj=0;jj<lin->GetEntries()-1;++jj) {
120  lout->Add(lin->At(jj)->Clone(lin->At(jj)->GetName()));
121  }
122  }
123  }
124 
125  if (fDoCopyHeader) {
126  AliAODHeader *out = (AliAODHeader*)eout->GetHeader();
127  AliAODHeader *in = (AliAODHeader*)evin->GetHeader();
128  *out = *in;
129  out->SetUniqueID(fTrials);
130  }
131  if (fDoCopyVZERO) {
132  AliAODVZERO *out = eout->GetVZEROData();
133  AliAODVZERO *in = evin->GetVZEROData();
134  *out = *in;
135  }
136  if (fDoCopyTZERO) {
137  AliAODTZERO *out = eout->GetTZEROData();
138  AliAODTZERO *in = evin->GetTZEROData();
139  *out = *in;
140  }
141  if (fDoCopyVertices) {
142  TClonesArray *out = eout->GetVertices();
143  TClonesArray *in = evin->GetVertices();
144  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
145  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
146  }
147  out->AbsorbObjects(in);
148  }
149  if (fDoCopyTOF) {
150  AliTOFHeader *out = const_cast<AliTOFHeader*>(eout->GetTOFHeader());
151  const AliTOFHeader *in = evin->GetTOFHeader();
152  *out = *in;
153  }
154  if (fDoCopyTracks) {
155  TClonesArray *out = eout->GetTracks();
156  TClonesArray *in = evin->GetTracks();
157  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
158  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
159  }
160  out->AbsorbObjects(in);
161  }
162  if (fDoCopyTrigger) {
163  AliAODCaloTrigger *out = eout->GetCaloTrigger("EMCAL");
164  AliAODCaloTrigger *in = evin->GetCaloTrigger("EMCAL");
165  *out = *in;
166  }
167  if (fDoCopyPTrigger) {
168  AliAODCaloTrigger *out = eout->GetCaloTrigger("PHOS");
169  AliAODCaloTrigger *in = evin->GetCaloTrigger("PHOS");
170  *out = *in;
171  }
172  if (fDoCopyCells) {
173  AliAODCaloCells *out = eout->GetEMCALCells();
174  AliAODCaloCells *in = evin->GetEMCALCells();
175  *out = *in;
176  }
177  if (fDoCopyPCells) {
178  AliAODCaloCells *out = eout->GetPHOSCells();
179  AliAODCaloCells *in = evin->GetPHOSCells();
180  *out = *in;
181  }
182  if (fDoCopyClusters) {
183  TClonesArray *out = eout->GetCaloClusters();
184  TClonesArray *in = evin->GetCaloClusters();
185  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
186  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
187  }
188  out->AbsorbObjects(in);
189  }
190 
191  if (fDoCopyDiMuons) {
192  TClonesArray *out = eout->GetDimuons();
193  TClonesArray *in = evin->GetDimuons();
194  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
195  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
196  }
197  out->AbsorbObjects(in);
198  }
199 
200  if (fDoCopyZDC) {
201  AliAODZDC *out = eout->GetZDCData();
202  AliAODZDC *in = evin->GetZDCData();
203  *out = *in;
204  }
205 
206  if (fDoCopyMC) {
207  TClonesArray *out = static_cast<TClonesArray*>(eout->FindListObject(AliAODMCParticle::StdBranchName()));
208  TClonesArray *in = static_cast<TClonesArray*>(evin->FindListObject(AliAODMCParticle::StdBranchName()));
209  if (in && !out) {
210  fgAODMCParticles = new TClonesArray("AliAODMCParticle",2*in->GetEntries());
211  fgAODMCParticles->SetName(AliAODMCParticle::StdBranchName());
212  oh->AddBranch("TClonesArray", &fgAODMCParticles);
213  out = static_cast<TClonesArray*>(eout->FindListObject(AliAODMCParticle::StdBranchName()));
214  }
215  if (in && out) {
216  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
217  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
218  }
219  out->AbsorbObjects(in);
220  if (fCutMC) {
221  for (Int_t i=0;i<out->GetEntriesFast();++i) {
222  AliAODMCParticle *mc = static_cast<AliAODMCParticle*>(in->At(i));
223  if ((mc==0)&&(i==0)) {
224  AliError(Form("%s: No MC info, skipping this event!",GetName()));
225  oh->SetFillAOD(kFALSE);
226  return;
227  }
228  if ((mc==0)||(TMath::Abs(mc->Y())>fYCutMC))
229  new ((*out)[i]) AliAODMCParticle;
230  }
231  }
232  }
233  }
234 
235  if (fDoCopyMCHeader) {
236  AliAODMCHeader *out = static_cast<AliAODMCHeader*>(eout->FindListObject(AliAODMCHeader::StdBranchName()));
237  AliAODMCHeader *in = static_cast<AliAODMCHeader*>(evin->FindListObject(AliAODMCHeader::StdBranchName()));
238  if (in && !out) {
239  fAODMcHeader = new AliAODMCHeader();
240  fAODMcHeader->SetName(AliAODMCHeader::StdBranchName());
241  oh->AddBranch("AliAODMCHeader",&fAODMcHeader);
242  out = static_cast<AliAODMCHeader*>(eout->FindListObject(AliAODMCHeader::StdBranchName()));
243  }
244  if (in && out) {
245  *out = *in;
246  if ((in->GetCrossSection()==0) && (fPyxsec>0)) {
247  out->SetCrossSection(fPyxsec);
248  out->SetTrials(fPytrials);
249  out->SetPtHard(fPypthardbin);
250  }
251  }
252  }
253 
254  if (gDebug>10) {
255  Int_t run = eout->GetRunNumber();
256  AliAODVertex *v=(AliAODVertex*)eout->GetVertices()->At(0);
257  Int_t vzn = v->GetNContributors();
258  Double_t vz = v->GetZ();
259  cout << "debug run " << run << " " << vzn << " " << vz << endl;
260  }
261 
262  fTrials = 0;
263  PostData(1, fOutputList);
264 }
265 
267 {
268  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
269  if (!tree) {
270  AliError(Form("%s: No current tree!",GetName()));
271  return kFALSE;
272  }
273 
274  Float_t xsection = 0;
275  Float_t trials = 0;
276  Int_t pthardbin = 0;
277 
278  TFile *curfile = tree->GetCurrentFile();
279  if (!curfile) {
280  AliError(Form("%s: No current file!",GetName()));
281  return kFALSE;
282  }
283 
284  TChain *chain = dynamic_cast<TChain*>(tree);
285  if (chain) tree = chain->GetTree();
286  Int_t nevents = tree->GetEntriesFast();
287 
288  Int_t slevel = gErrorIgnoreLevel;
289  gErrorIgnoreLevel = kFatal;
290  Bool_t res = PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
291  gErrorIgnoreLevel=slevel;
292 
293  if (res) {
294  cout << "AliAodSkimTask " << GetName() << " found xsec info: " << xsection << " " << trials << " " << pthardbin << " " << nevents << endl;
295  fPyxsec = xsection;
296  fPytrials = trials;
297  fPypthardbin = pthardbin;
298  }
299 
300  return res;
301 }
302 
304 {
305  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
306  if (man->GetAnalysisType()!=0)
307  return;
308  if (fHevs==0)
309  return;
310  cout << "AliAodSkimTask " << GetName() << " terminated with accepted fraction of events: " << fHevs->GetBinContent(2)/fHevs->GetEntries()
311  << " (" << fHevs->GetBinContent(2) << "/" << fHevs->GetEntries() << ")" << endl;
312 }
313 
314 Bool_t AliAodSkimTask::PythiaInfoFromFile(const char* currFile, Float_t &xsec, Float_t &trials, Int_t &pthard)
315 {
316  TString file(currFile);
317  xsec = 0;
318  trials = 1;
319 
320  if (file.Contains(".zip#")) {
321  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
322  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
323  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
324  file.Replace(pos+1,pos2-pos1,"");
325  } else {
326  // not an archive take the basename....
327  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
328  }
329  AliDebug(1,Form("File name: %s",file.Data()));
330 
331  // Get the pt hard bin
332  TString strPthard(file);
333 
334  strPthard.Remove(strPthard.Last('/'));
335  strPthard.Remove(strPthard.Last('/'));
336  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
337  strPthard.Remove(0,strPthard.Last('/')+1);
338  if (strPthard.IsDec()) {
339  pthard = strPthard.Atoi();
340  }
341  else {
342  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
343  pthard = -1;
344  }
345 
346  // problem that we cannot really test the existance of a file in a archive so we have to live with open error message from root
347  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
348 
349  if (!fxsec) {
350  // next trial fetch the histgram file
351  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
352  if (!fxsec) {
353  // not a severe condition but inciate that we have no information
354  return kFALSE;
355  } else {
356  // find the tlist we want to be independtent of the name so use the Tkey
357  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
358  if (!key) {
359  fxsec->Close();
360  return kFALSE;
361  }
362  TList *list = dynamic_cast<TList*>(key->ReadObj());
363  if (!list) {
364  fxsec->Close();
365  return kFALSE;
366  }
367  xsec = static_cast<TProfile*>(list->FindObject("h1Xsec"))->GetBinContent(1);
368  trials = static_cast<TH1F*>(list->FindObject("h1Trials"))->GetBinContent(1);
369  fxsec->Close();
370  }
371  } else { // no tree pyxsec.root
372  TTree *xtree = static_cast<TTree*>(fxsec->Get("Xsection"));
373  if (!xtree) {
374  fxsec->Close();
375  return kFALSE;
376  }
377  UInt_t ntrials = 0;
378  Double_t xsection = 0;
379  xtree->SetBranchAddress("xsection",&xsection);
380  xtree->SetBranchAddress("ntrials",&ntrials);
381  xtree->GetEntry(0);
382  trials = ntrials;
383  xsec = xsection;
384  fxsec->Close();
385  }
386  return kTRUE;
387 }
388 
389 const char *AliAodSkimTask::Str() const
390 {
391  return Form("mine%.2f_%dycut%.2f_%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d",
392  fClusMinE,
393  fCutMC,
394  fYCutMC,
396  fDoCopyVZERO,
397  fDoCopyTZERO,
399  fDoCopyTOF,
403  fDoCopyCells,
407  fDoCopyZDC,
408  fDoCopyMC,
410 }
Double_t fYCutMC
Bool_t fDoCopyTrigger
Int_t fPypthardbin
pythia trials
Bool_t fDoCopyPTrigger
double Double_t
Definition: External.C:58
TSystem * gSystem
TH1F * fHevs
output list
Float_t fPyxsec
events seen since last acceptance
void Terminate(Option_t *option)
AliAODEvent * fAOD
pythia pthard bin
Float_t fPytrials
pythia xsection
Bool_t fDoCopyTracks
Bool_t fDoCopyHeader
Double_t fClusMinE
int Int_t
Definition: External.C:63
const char * Str() const
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
void UserExec(Option_t *option)
virtual ~AliAodSkimTask()
const char * GetVersion() const
cluster distribution
Bool_t fDoCopyClusters
Use to skim AOD files.
void UserCreateOutputObjects()
TList * fOutputList
MC header.
Bool_t fDoCopyPCells
TH1F * fHclus
events processed/accepted
Bool_t fDoCopyVertices
AliAODMCHeader * fAODMcHeader
input event
Bool_t PythiaInfoFromFile(const char *currFile, Float_t &xsec, Float_t &trials, Int_t &pthard)
Bool_t fDoCopyDiMuons
TFile * file
TList with histograms for a given trigger.
Int_t nevents[nsamples]
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
TFile * fout
input train file
Bool_t fDoCopyMCHeader