AliPhysics  ef3b16e (ef3b16e)
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), fDoCopyTracklets(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), fDoCopyTracklets(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 (fDoCopyTracklets) {
163  AliAODTracklets *out = eout->GetTracklets();
164  AliAODTracklets *in = evin->GetTracklets();
165  *out = *in;
166  }
167  if (fDoCopyTrigger) {
168  AliAODCaloTrigger *out = eout->GetCaloTrigger("EMCAL");
169  AliAODCaloTrigger *in = evin->GetCaloTrigger("EMCAL");
170  *out = *in;
171  }
172  if (fDoCopyPTrigger) {
173  AliAODCaloTrigger *out = eout->GetCaloTrigger("PHOS");
174  AliAODCaloTrigger *in = evin->GetCaloTrigger("PHOS");
175  *out = *in;
176  }
177  if (fDoCopyCells) {
178  AliAODCaloCells *out = eout->GetEMCALCells();
179  AliAODCaloCells *in = evin->GetEMCALCells();
180  *out = *in;
181  }
182  if (fDoCopyPCells) {
183  AliAODCaloCells *out = eout->GetPHOSCells();
184  AliAODCaloCells *in = evin->GetPHOSCells();
185  *out = *in;
186  }
187  if (fDoCopyClusters) {
188  TClonesArray *out = eout->GetCaloClusters();
189  TClonesArray *in = evin->GetCaloClusters();
190  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
191  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
192  }
193  out->AbsorbObjects(in);
194  }
195 
196  if (fDoCopyDiMuons) {
197  TClonesArray *out = eout->GetDimuons();
198  TClonesArray *in = evin->GetDimuons();
199  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
200  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
201  }
202  out->AbsorbObjects(in);
203  }
204 
205  if (fDoCopyZDC) {
206  AliAODZDC *out = eout->GetZDCData();
207  AliAODZDC *in = evin->GetZDCData();
208  *out = *in;
209  }
210 
211  if (fDoCopyMC) {
212  TClonesArray *out = static_cast<TClonesArray*>(eout->FindListObject(AliAODMCParticle::StdBranchName()));
213  TClonesArray *in = static_cast<TClonesArray*>(evin->FindListObject(AliAODMCParticle::StdBranchName()));
214  if (in && !out) {
215  fgAODMCParticles = new TClonesArray("AliAODMCParticle",2*in->GetEntries());
216  fgAODMCParticles->SetName(AliAODMCParticle::StdBranchName());
217  oh->AddBranch("TClonesArray", &fgAODMCParticles);
218  out = static_cast<TClonesArray*>(eout->FindListObject(AliAODMCParticle::StdBranchName()));
219  }
220  if (in && out) {
221  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
222  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
223  }
224  out->AbsorbObjects(in);
225  if (fCutMC) {
226  for (Int_t i=0;i<out->GetEntriesFast();++i) {
227  AliAODMCParticle *mc = static_cast<AliAODMCParticle*>(in->At(i));
228  if ((mc==0)&&(i==0)) {
229  AliError(Form("%s: No MC info, skipping this event!",GetName()));
230  oh->SetFillAOD(kFALSE);
231  return;
232  }
233  if ((mc==0)||(TMath::Abs(mc->Y())>fYCutMC))
234  new ((*out)[i]) AliAODMCParticle;
235  }
236  }
237  }
238  }
239 
240  if (fDoCopyMCHeader) {
241  AliAODMCHeader *out = static_cast<AliAODMCHeader*>(eout->FindListObject(AliAODMCHeader::StdBranchName()));
242  AliAODMCHeader *in = static_cast<AliAODMCHeader*>(evin->FindListObject(AliAODMCHeader::StdBranchName()));
243  if (in && !out) {
244  fAODMcHeader = new AliAODMCHeader();
245  fAODMcHeader->SetName(AliAODMCHeader::StdBranchName());
246  oh->AddBranch("AliAODMCHeader",&fAODMcHeader);
247  out = static_cast<AliAODMCHeader*>(eout->FindListObject(AliAODMCHeader::StdBranchName()));
248  }
249  if (in && out) {
250  *out = *in;
251  if ((in->GetCrossSection()==0) && (fPyxsec>0)) {
252  out->SetCrossSection(fPyxsec);
253  out->SetTrials(fPytrials);
254  out->SetPtHard(fPypthardbin);
255  }
256  }
257  }
258 
259  if (gDebug>10) {
260  Int_t run = eout->GetRunNumber();
261  AliAODVertex *v=(AliAODVertex*)eout->GetVertices()->At(0);
262  Int_t vzn = v->GetNContributors();
263  Double_t vz = v->GetZ();
264  cout << "debug run " << run << " " << vzn << " " << vz << endl;
265  }
266 
267  fTrials = 0;
268  PostData(1, fOutputList);
269 }
270 
272 {
273  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
274  if (!tree) {
275  AliError(Form("%s: No current tree!",GetName()));
276  return kFALSE;
277  }
278 
279  Float_t xsection = 0;
280  Float_t trials = 0;
281  Int_t pthardbin = 0;
282 
283  TFile *curfile = tree->GetCurrentFile();
284  if (!curfile) {
285  AliError(Form("%s: No current file!",GetName()));
286  return kFALSE;
287  }
288 
289  TChain *chain = dynamic_cast<TChain*>(tree);
290  if (chain) tree = chain->GetTree();
291  Int_t nevents = tree->GetEntriesFast();
292 
293  Int_t slevel = gErrorIgnoreLevel;
294  gErrorIgnoreLevel = kFatal;
295  Bool_t res = PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
296  gErrorIgnoreLevel=slevel;
297 
298  if (res) {
299  cout << "AliAodSkimTask " << GetName() << " found xsec info: " << xsection << " " << trials << " " << pthardbin << " " << nevents << endl;
300  fPyxsec = xsection;
301  fPytrials = trials;
302  fPypthardbin = pthardbin;
303  }
304 
305  return res;
306 }
307 
309 {
310  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
311  if (man->GetAnalysisType()!=0)
312  return;
313  if (fHevs==0)
314  return;
315  cout << "AliAodSkimTask " << GetName() << " terminated with accepted fraction of events: " << fHevs->GetBinContent(2)/fHevs->GetEntries()
316  << " (" << fHevs->GetBinContent(2) << "/" << fHevs->GetEntries() << ")" << endl;
317 }
318 
319 Bool_t AliAodSkimTask::PythiaInfoFromFile(const char* currFile, Float_t &xsec, Float_t &trials, Int_t &pthard)
320 {
321  TString file(currFile);
322  xsec = 0;
323  trials = 1;
324 
325  if (file.Contains(".zip#")) {
326  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
327  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
328  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
329  file.Replace(pos+1,pos2-pos1,"");
330  } else {
331  // not an archive take the basename....
332  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
333  }
334  AliDebug(1,Form("File name: %s",file.Data()));
335 
336  // Get the pt hard bin
337  TString strPthard(file);
338 
339  strPthard.Remove(strPthard.Last('/'));
340  strPthard.Remove(strPthard.Last('/'));
341  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
342  strPthard.Remove(0,strPthard.Last('/')+1);
343  if (strPthard.IsDec()) {
344  pthard = strPthard.Atoi();
345  }
346  else {
347  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
348  pthard = -1;
349  }
350 
351  // 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
352  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
353 
354  if (!fxsec) {
355  // next trial fetch the histgram file
356  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
357  if (!fxsec) {
358  // not a severe condition but inciate that we have no information
359  return kFALSE;
360  } else {
361  // find the tlist we want to be independtent of the name so use the Tkey
362  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
363  if (!key) {
364  fxsec->Close();
365  return kFALSE;
366  }
367  TList *list = dynamic_cast<TList*>(key->ReadObj());
368  if (!list) {
369  fxsec->Close();
370  return kFALSE;
371  }
372  xsec = static_cast<TProfile*>(list->FindObject("h1Xsec"))->GetBinContent(1);
373  trials = static_cast<TH1F*>(list->FindObject("h1Trials"))->GetBinContent(1);
374  fxsec->Close();
375  }
376  } else { // no tree pyxsec.root
377  TTree *xtree = static_cast<TTree*>(fxsec->Get("Xsection"));
378  if (!xtree) {
379  fxsec->Close();
380  return kFALSE;
381  }
382  UInt_t ntrials = 0;
383  Double_t xsection = 0;
384  xtree->SetBranchAddress("xsection",&xsection);
385  xtree->SetBranchAddress("ntrials",&ntrials);
386  xtree->GetEntry(0);
387  trials = ntrials;
388  xsec = xsection;
389  fxsec->Close();
390  }
391  return kTRUE;
392 }
393 
394 const char *AliAodSkimTask::Str() const
395 {
396  return Form("mine%.2f_%dycut%.2f_%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d",
397  fClusMinE,
398  fCutMC,
399  fYCutMC,
401  fDoCopyVZERO,
402  fDoCopyTZERO,
404  fDoCopyTOF,
409  fDoCopyCells,
413  fDoCopyZDC,
414  fDoCopyMC,
416 }
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()
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