AliPhysics  63e47e1 (63e47e1)
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), fGammaBr(""),
25  fDoCopyHeader(1), fDoCopyVZERO(1), fDoCopyTZERO(1), fDoCopyVertices(1), fDoCopyTOF(1), fDoCopyTracklets(1), fDoCopyTracks(1), fDoCopyTrigger(1), fDoCopyPTrigger(0),
26  fDoCopyCells(1), fDoCopyPCells(0), fDoCopyClusters(1), fDoCopyDiMuons(0), fDoCopyTrdTracks(0), fDoCopyV0s(0), fDoCopyCascades(0), fDoCopyZDC(1), fDoCopyConv(0),
27  fDoCopyMC(1), fDoCopyMCHeader(1), fTrials(0), fPyxsec(0), 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), fGammaBr(""),
33  fDoCopyHeader(1), fDoCopyVZERO(1), fDoCopyTZERO(1), fDoCopyVertices(1), fDoCopyTOF(1), fDoCopyTracklets(1), fDoCopyTracks(1), fDoCopyTrigger(1), fDoCopyPTrigger(0),
34  fDoCopyCells(1), fDoCopyPCells(0), fDoCopyClusters(1), fDoCopyDiMuons(0), fDoCopyTrdTracks(0), fDoCopyV0s(0), fDoCopyCascades(0), fDoCopyZDC(1), fDoCopyConv(0),
35  fDoCopyMC(1), fDoCopyMCHeader(1), fTrials(0), fPyxsec(0), 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 
132  if (fDoCopyVZERO) {
133  AliAODVZERO *out = eout->GetVZEROData();
134  AliAODVZERO *in = evin->GetVZEROData();
135  *out = *in;
136  }
137 
138  if (fDoCopyTZERO) {
139  AliAODTZERO *out = eout->GetTZEROData();
140  AliAODTZERO *in = evin->GetTZEROData();
141  *out = *in;
142  }
143 
144  if (fDoCopyVertices) {
145  TClonesArray *out = eout->GetVertices();
146  TClonesArray *in = evin->GetVertices();
147  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
148  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
149  }
150  out->AbsorbObjects(in);
151  }
152 
153  if (fDoCopyTOF) {
154  AliTOFHeader *out = const_cast<AliTOFHeader*>(eout->GetTOFHeader());
155  const AliTOFHeader *in = evin->GetTOFHeader();
156  *out = *in;
157  }
158 
159  if (fDoCopyTracks) {
160  TClonesArray *out = eout->GetTracks();
161  TClonesArray *in = evin->GetTracks();
162  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
163  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
164  }
165  out->AbsorbObjects(in);
166  }
167 
168  if (fDoCopyTracklets) {
169  AliAODTracklets *out = eout->GetTracklets();
170  AliAODTracklets *in = evin->GetTracklets();
171  *out = *in;
172  }
173 
174  if (fDoCopyTrigger) {
175  AliAODCaloTrigger *out = eout->GetCaloTrigger("EMCAL");
176  AliAODCaloTrigger *in = evin->GetCaloTrigger("EMCAL");
177  *out = *in;
178  }
179 
180  if (fDoCopyPTrigger) {
181  AliAODCaloTrigger *out = eout->GetCaloTrigger("PHOS");
182  AliAODCaloTrigger *in = evin->GetCaloTrigger("PHOS");
183  *out = *in;
184  }
185 
186  if (fDoCopyCells) {
187  AliAODCaloCells *out = eout->GetEMCALCells();
188  AliAODCaloCells *in = evin->GetEMCALCells();
189  *out = *in;
190  }
191 
192  if (fDoCopyPCells) {
193  AliAODCaloCells *out = eout->GetPHOSCells();
194  AliAODCaloCells *in = evin->GetPHOSCells();
195  *out = *in;
196  }
197 
198  if (fDoCopyClusters) {
199  TClonesArray *out = eout->GetCaloClusters();
200  TClonesArray *in = evin->GetCaloClusters();
201  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
202  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
203  }
204  out->AbsorbObjects(in);
205  }
206 
207  if (fDoCopyDiMuons) {
208  TClonesArray *out = eout->GetDimuons();
209  TClonesArray *in = evin->GetDimuons();
210  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
211  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
212  }
213  out->AbsorbObjects(in);
214  }
215 
216  if (fDoCopyTrdTracks) {
217  TClonesArray *out = static_cast<TClonesArray*>(eout->FindListObject("trdTracks"));
218  TClonesArray *in = static_cast<TClonesArray*>(eout->FindListObject("trdTracks"));
219  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
220  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
221  }
222  out->AbsorbObjects(in);
223  }
224 
225  if (fDoCopyV0s) {
226  TClonesArray *out = eout->GetV0s();
227  TClonesArray *in = evin->GetV0s();
228  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
229  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
230  }
231  out->AbsorbObjects(in);
232  }
233 
234  if (fDoCopyCascades) {
235  TClonesArray *out = eout->GetCascades();
236  TClonesArray *in = evin->GetCascades();
237  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
238  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
239  }
240  out->AbsorbObjects(in);
241  }
242 
243  if (fDoCopyZDC) {
244  AliAODZDC *out = eout->GetZDCData();
245  AliAODZDC *in = evin->GetZDCData();
246  *out = *in;
247  }
248 
249  if (fDoCopyDiMuons) {
250  TClonesArray *out = eout->GetDimuons();
251  TClonesArray *in = evin->GetDimuons();
252  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
253  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
254  }
255  out->AbsorbObjects(in);
256  }
257 
258  if (fDoCopyConv) {
259  TClonesArray *out = dynamic_cast<TClonesArray*>(eout->FindListObject(fGammaBr));
260  TClonesArray *in = dynamic_cast<TClonesArray*>(evin->FindListObject(fGammaBr));
261  if (!in) {
262  evin->GetList()->ls();
263  AliFatal(Form("%s: Could not find conversion branch with name %s!",GetName(), fGammaBr.Data()));
264  }
265  if (in && !out) {
266  out = new TClonesArray("AliAODConversionPhoton",2*in->GetEntries());
267  out->SetName(fGammaBr);
268  oh->AddBranch("TClonesArray", &out);
269  }
270  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
271  out->Delete();
272  }
273  out->AbsorbObjects(in);
274  }
275 
276  if (fDoCopyMC) {
277  TClonesArray *out = static_cast<TClonesArray*>(eout->FindListObject(AliAODMCParticle::StdBranchName()));
278  TClonesArray *in = static_cast<TClonesArray*>(evin->FindListObject(AliAODMCParticle::StdBranchName()));
279  if (in && !out) {
280  fgAODMCParticles = new TClonesArray("AliAODMCParticle",2*in->GetEntries());
281  fgAODMCParticles->SetName(AliAODMCParticle::StdBranchName());
282  oh->AddBranch("TClonesArray", &fgAODMCParticles);
283  out = static_cast<TClonesArray*>(eout->FindListObject(AliAODMCParticle::StdBranchName()));
284  }
285  if (in && out) {
286  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
287  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
288  }
289  out->AbsorbObjects(in);
290  if (fCutMC) {
291  for (Int_t i=0;i<out->GetEntriesFast();++i) {
292  AliAODMCParticle *mc = static_cast<AliAODMCParticle*>(in->At(i));
293  if ((mc==0)&&(i==0)) {
294  AliError(Form("%s: No MC info, skipping this event!",GetName()));
295  oh->SetFillAOD(kFALSE);
296  return;
297  }
298  if ((mc==0)||(TMath::Abs(mc->Y())>fYCutMC))
299  new ((*out)[i]) AliAODMCParticle;
300  }
301  }
302  }
303  }
304 
305  if (fDoCopyMCHeader) {
306  AliAODMCHeader *out = static_cast<AliAODMCHeader*>(eout->FindListObject(AliAODMCHeader::StdBranchName()));
307  AliAODMCHeader *in = static_cast<AliAODMCHeader*>(evin->FindListObject(AliAODMCHeader::StdBranchName()));
308  if (in && !out) {
309  fAODMcHeader = new AliAODMCHeader();
310  fAODMcHeader->SetName(AliAODMCHeader::StdBranchName());
311  oh->AddBranch("AliAODMCHeader",&fAODMcHeader);
312  out = static_cast<AliAODMCHeader*>(eout->FindListObject(AliAODMCHeader::StdBranchName()));
313  }
314  if (in && out) {
315  *out = *in;
316  if ((in->GetCrossSection()==0) && (fPyxsec>0)) {
317  out->SetCrossSection(fPyxsec);
318  out->SetTrials(fPytrials);
319  out->SetPtHard(fPypthardbin);
320  }
321  }
322  }
323 
324  if (gDebug>10) {
325  Int_t run = eout->GetRunNumber();
326  AliAODVertex *v=(AliAODVertex*)eout->GetVertices()->At(0);
327  Int_t vzn = v->GetNContributors();
328  Double_t vz = v->GetZ();
329  cout << "debug run " << run << " " << vzn << " " << vz << endl;
330  }
331 
332  fTrials = 0;
333  PostData(1, fOutputList);
334 }
335 
337 {
338  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
339  if (!tree) {
340  AliError(Form("%s: No current tree!",GetName()));
341  return kFALSE;
342  }
343 
344  Float_t xsection = 0;
345  Float_t trials = 0;
346  Int_t pthardbin = 0;
347 
348  TFile *curfile = tree->GetCurrentFile();
349  if (!curfile) {
350  AliError(Form("%s: No current file!",GetName()));
351  return kFALSE;
352  }
353 
354  TChain *chain = dynamic_cast<TChain*>(tree);
355  if (chain) tree = chain->GetTree();
356  Int_t nevents = tree->GetEntriesFast();
357 
358  Int_t slevel = gErrorIgnoreLevel;
359  gErrorIgnoreLevel = kFatal;
360  Bool_t res = PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
361  gErrorIgnoreLevel=slevel;
362 
363  if (res) {
364  cout << "AliAodSkimTask " << GetName() << " found xsec info: " << xsection << " " << trials << " " << pthardbin << " " << nevents << endl;
365  fPyxsec = xsection;
366  fPytrials = trials;
367  fPypthardbin = pthardbin;
368  }
369 
370  return res;
371 }
372 
374 {
375  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
376  if (man->GetAnalysisType()!=0)
377  return;
378  if (fHevs==0)
379  return;
380  cout << "AliAodSkimTask " << GetName() << " terminated with accepted fraction of events: " << fHevs->GetBinContent(2)/fHevs->GetEntries()
381  << " (" << fHevs->GetBinContent(2) << "/" << fHevs->GetEntries() << ")" << endl;
382 }
383 
384 Bool_t AliAodSkimTask::PythiaInfoFromFile(const char* currFile, Float_t &xsec, Float_t &trials, Int_t &pthard)
385 {
386  TString file(currFile);
387  xsec = 0;
388  trials = 1;
389 
390  if (file.Contains(".zip#")) {
391  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
392  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
393  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
394  file.Replace(pos+1,pos2-pos1,"");
395  } else {
396  // not an archive take the basename....
397  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
398  }
399  AliDebug(1,Form("File name: %s",file.Data()));
400 
401  // Get the pt hard bin
402  TString strPthard(file);
403 
404  strPthard.Remove(strPthard.Last('/'));
405  strPthard.Remove(strPthard.Last('/'));
406  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
407  strPthard.Remove(0,strPthard.Last('/')+1);
408  if (strPthard.IsDec()) {
409  pthard = strPthard.Atoi();
410  }
411  else {
412  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
413  pthard = -1;
414  }
415 
416  // 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
417  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
418 
419  if (!fxsec) {
420  // next trial fetch the histgram file
421  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
422  if (!fxsec) {
423  // not a severe condition but inciate that we have no information
424  return kFALSE;
425  } else {
426  // find the tlist we want to be independtent of the name so use the Tkey
427  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
428  if (!key) {
429  fxsec->Close();
430  return kFALSE;
431  }
432  TList *list = dynamic_cast<TList*>(key->ReadObj());
433  if (!list) {
434  fxsec->Close();
435  return kFALSE;
436  }
437  xsec = static_cast<TProfile*>(list->FindObject("h1Xsec"))->GetBinContent(1);
438  trials = static_cast<TH1F*>(list->FindObject("h1Trials"))->GetBinContent(1);
439  fxsec->Close();
440  }
441  } else { // no tree pyxsec.root
442  TTree *xtree = static_cast<TTree*>(fxsec->Get("Xsection"));
443  if (!xtree) {
444  fxsec->Close();
445  return kFALSE;
446  }
447  UInt_t ntrials = 0;
448  Double_t xsection = 0;
449  xtree->SetBranchAddress("xsection",&xsection);
450  xtree->SetBranchAddress("ntrials",&ntrials);
451  xtree->GetEntry(0);
452  trials = ntrials;
453  xsec = xsection;
454  fxsec->Close();
455  }
456  return kTRUE;
457 }
458 
459 const char *AliAodSkimTask::Str() const
460 {
461  return Form("mine%.2f_%dycut%.2f_%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d",
462  fClusMinE,
463  fCutMC,
464  fYCutMC,
466  fDoCopyVZERO,
467  fDoCopyTZERO,
469  fDoCopyTOF,
474  fDoCopyCells,
479  fDoCopyV0s,
481  fDoCopyZDC,
482  fDoCopyConv,
483  fDoCopyMC,
485 }
486 
Double_t fYCutMC
Bool_t fDoCopyTrigger
Int_t fPypthardbin
pythia trials
Bool_t fDoCopyPTrigger
Bool_t fDoCopyTracklets
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()
Bool_t fDoCopyCascades
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_t fDoCopyTrdTracks
bool Bool_t
Definition: External.C:53
TFile * fout
input train file
Bool_t fDoCopyMCHeader