AliPhysics  a6017e1 (a6017e1)
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 
23 AliAodSkimTask::AliAodSkimTask(const char* name) :
24  AliAnalysisTaskSE(name), fClusMinE(-1), fCutMC(1), fYCutMC(0.7), fCutMinPt(0), fCutFilterBit(-1), fGammaBr(""),
25  fDoCopyHeader(1), fDoCopyVZERO(1), fDoCopyTZERO(1), fDoCopyVertices(1), fDoCopyTOF(1), fDoCopyTracklets(1), fDoCopyTracks(1), fDoRemoveTracks(0), fDoCleanTracks(0),
26  fDoRemCovMat(0), fDoRemPid(0), fDoCopyTrigger(1), fDoCopyPTrigger(0), fDoCopyCells(1), fDoCopyPCells(0), fDoCopyClusters(1), fDoCopyDiMuons(0), fDoCopyTrdTracks(0),
27  fDoCopyV0s(0), fDoCopyCascades(0), fDoCopyZDC(1), fDoCopyConv(0), fDoCopyMC(1), fDoCopyMCHeader(1),
28  fTrials(0), fPyxsec(0), fPytrials(0), fPypthardbin(0), fAOD(0), fAODMcHeader(0), fOutputList(0), fHevs(0), fHclus(0)
29 {
30  if (name) {
31  DefineInput(0, TChain::Class());
32  DefineOutput(1, TList::Class());
33  }
34 }
35 
37 {
38  if (fOutputList) {
39  delete fOutputList;
40  }
41  delete fHevs;
42  delete fHclus;
43 }
44 
46 {
47  if (t->IsMuonGlobalTrack())
48  return kFALSE;
49  if (t->Pt()<fCutMinPt)
50  return kFALSE;
51  if (t->TestFilterBit(fCutFilterBit)==0)
52  return kFALSE;
53 
54  return kTRUE;
55 }
56 
57 void AliAodSkimTask::CleanTrack(AliAODTrack *t)
58 {
59  if (!fDoCleanTracks)
60  return;
61  //cout << "Clean tracks " << endl;
62  t->SetRAtAbsorberEnd(0);
63  t->SetChi2MatchTrigger(0);
64  t->SetMuonClusterMap(0);
65  t->SetITSMuonClusterMap(0);
66  t->SetMUONtrigHitsMapTrg(0);
67  t->SetMUONtrigHitsMapTrk(0);
68  t->SetMFTClusterPattern(0);
69  t->SetMatchTrigger(0);
70  t->SetIsMuonGlobalTrack(0);
71  t->SetXYAtDCA(-999., -999.);
72  t->SetPxPyPzAtDCA(-999., -999., -999.);
73  if (fDoRemCovMat)
74  t->RemoveCovMatrix();
75  if (fDoRemPid) {
76  AliAODPid *pid = t->GetDetPid();
77  delete pid;
78  t->SetDetPID(0);
79  }
80 }
81 
83 {
84 
85  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
86  AliAODHandler *oh = (AliAODHandler*)man->GetOutputEventHandler();
87  if (oh) {
88  TFile *fout = oh->GetTree()->GetCurrentFile();
89  fout->SetCompressionLevel(2);
90  }
91 
92  fOutputList = new TList;
93  fOutputList->SetOwner();
94  fHevs = new TH1F("hEvs","",2,-0.5,1.5);
95  fOutputList->Add(fHevs);
96  fHclus = new TH1F("hClus",";E (GeV)",200,0,100);
97  fOutputList->Add(fHclus);
98  PostData(1, fOutputList);
99 }
100 
102 {
103  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
104  if (!fAOD)
105  return;
106 
107  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
108  AliAODHandler *oh = (AliAODHandler*)man->GetOutputEventHandler();
109  if (!oh) {
110  AliFatal(Form("%s: No output handler found", GetName()));
111  return;
112  }
113  oh->SetFillAOD(kFALSE);
114 
115  Bool_t store = kFALSE;
116  if (fClusMinE>0) {
117  TClonesArray *cls = fAOD->GetCaloClusters();
118  for (Int_t i=0; i<cls->GetEntriesFast(); ++i) {
119  AliAODCaloCluster *clus = static_cast<AliAODCaloCluster*>(cls->At(i));
120  if (!clus->IsEMCAL())
121  continue;
122  Double_t e = clus->E();
123  fHclus->Fill(e);
124  if (e>fClusMinE) {
125  store = kTRUE;
126  }
127  }
128  } else {
129  store = kTRUE;
130  }
131 
132  if (!store) {
133  ++fTrials;
134  fHevs->Fill(0);
135  return;
136  }
137 
138  fHevs->Fill(1);
139 
140  oh->SetFillAOD(kTRUE);
141  AliAODEvent *eout = dynamic_cast<AliAODEvent*>(oh->GetAOD());
142  AliAODEvent *evin = dynamic_cast<AliAODEvent*>(InputEvent());
143  TTree *tout = oh->GetTree();
144  if (tout) {
145  TList *lout = tout->GetUserInfo();
146  if (lout->FindObject("alirootVersion")==0) {
147  TList *lin = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetUserInfo();
148  TString apver(gSystem->BaseName(gSystem->Getenv("ALICE_PHYSICS")));
149  lout->Add(new TObjString(Form("AodSkim: ver %s, tag %s with settings %s",GetVersion(),apver.Data(),Str())));
150  AliInfo(Form("%s: Set user info %s", GetName(), lout->At(0)->GetName()));
151  for (Int_t jj=0;jj<lin->GetEntries()-1;++jj) {
152  lout->Add(lin->At(jj)->Clone(lin->At(jj)->GetName()));
153  }
154  }
155  }
156 
157  if (fDoCopyHeader) {
158  AliAODHeader *out = (AliAODHeader*)eout->GetHeader();
159  AliAODHeader *in = (AliAODHeader*)evin->GetHeader();
160  *out = *in;
161  out->SetUniqueID(fTrials);
162  }
163 
164  if (fDoCopyVZERO) {
165  AliAODVZERO *out = eout->GetVZEROData();
166  AliAODVZERO *in = evin->GetVZEROData();
167  *out = *in;
168  }
169 
170  if (fDoCopyTZERO) {
171  AliAODTZERO *out = eout->GetTZEROData();
172  AliAODTZERO *in = evin->GetTZEROData();
173  *out = *in;
174  }
175 
176  if (fDoCopyVertices) {
177  TClonesArray *out = eout->GetVertices();
178  TClonesArray *in = evin->GetVertices();
179  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
180  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
181  }
182  out->AbsorbObjects(in);
183  }
184 
185  if (fDoCopyTOF) {
186  AliTOFHeader *out = const_cast<AliTOFHeader*>(eout->GetTOFHeader());
187  const AliTOFHeader *in = evin->GetTOFHeader();
188  *out = *in;
189  }
190 
191  if (fDoCopyTracks) {
192  TClonesArray *out = eout->GetTracks();
193  TClonesArray *in = evin->GetTracks();
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  if (fDoRemoveTracks||fDoCleanTracks) {
199  for (Int_t i=0;i<out->GetEntriesFast();++i) {
200  AliAODTrack *t = static_cast<AliAODTrack*>(out->At(i));
201  if (KeepTrack(t))
202  CleanTrack(t);
203  else {
204  new ((*out)[i]) AliAODTrack;
205  }
206  }
207  }
208  }
209 
210  if (fDoCopyTracklets) {
211  AliAODTracklets *out = eout->GetTracklets();
212  AliAODTracklets *in = evin->GetTracklets();
213  *out = *in;
214  }
215 
216  if (fDoCopyTrigger) {
217  AliAODCaloTrigger *out = eout->GetCaloTrigger("EMCAL");
218  AliAODCaloTrigger *in = evin->GetCaloTrigger("EMCAL");
219  *out = *in;
220  }
221 
222  if (fDoCopyPTrigger) {
223  AliAODCaloTrigger *out = eout->GetCaloTrigger("PHOS");
224  AliAODCaloTrigger *in = evin->GetCaloTrigger("PHOS");
225  *out = *in;
226  }
227 
228  if (fDoCopyCells) {
229  AliAODCaloCells *out = eout->GetEMCALCells();
230  AliAODCaloCells *in = evin->GetEMCALCells();
231  *out = *in;
232  }
233 
234  if (fDoCopyPCells) {
235  AliAODCaloCells *out = eout->GetPHOSCells();
236  AliAODCaloCells *in = evin->GetPHOSCells();
237  *out = *in;
238  }
239 
240  if (fDoCopyClusters) {
241  TClonesArray *out = eout->GetCaloClusters();
242  TClonesArray *in = evin->GetCaloClusters();
243  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
244  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
245  }
246  out->AbsorbObjects(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 (fDoCopyTrdTracks) {
259  TClonesArray *out = static_cast<TClonesArray*>(eout->FindListObject("trdTracks"));
260  TClonesArray *in = static_cast<TClonesArray*>(eout->FindListObject("trdTracks"));
261  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
262  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
263  }
264  out->AbsorbObjects(in);
265  }
266 
267  if (fDoCopyV0s) {
268  TClonesArray *out = eout->GetV0s();
269  TClonesArray *in = evin->GetV0s();
270  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
271  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
272  }
273  out->AbsorbObjects(in);
274  }
275 
276  if (fDoCopyCascades) {
277  TClonesArray *out = eout->GetCascades();
278  TClonesArray *in = evin->GetCascades();
279  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
280  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
281  }
282  out->AbsorbObjects(in);
283  }
284 
285  if (fDoCopyZDC) {
286  AliAODZDC *out = eout->GetZDCData();
287  AliAODZDC *in = evin->GetZDCData();
288  *out = *in;
289  }
290 
291  if (fDoCopyDiMuons) {
292  TClonesArray *out = eout->GetDimuons();
293  TClonesArray *in = evin->GetDimuons();
294  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
295  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
296  }
297  out->AbsorbObjects(in);
298  }
299 
300  if (fDoCopyConv) {
301  TClonesArray *out = dynamic_cast<TClonesArray*>(eout->FindListObject(fGammaBr));
302  TClonesArray *in = dynamic_cast<TClonesArray*>(evin->FindListObject(fGammaBr));
303  if (!in) {
304  evin->GetList()->ls();
305  AliFatal(Form("%s: Could not find conversion branch with name %s!",GetName(), fGammaBr.Data()));
306  }
307  if (in && !out) {
308  out = new TClonesArray("AliAODConversionPhoton",2*in->GetEntries());
309  out->SetName(fGammaBr);
310  oh->AddBranch("TClonesArray", &out);
311  }
312  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
313  out->Delete();
314  }
315  out->AbsorbObjects(in);
316  }
317 
318  if (fDoCopyMC) {
319  TClonesArray *out = static_cast<TClonesArray*>(eout->FindListObject(AliAODMCParticle::StdBranchName()));
320  TClonesArray *in = static_cast<TClonesArray*>(evin->FindListObject(AliAODMCParticle::StdBranchName()));
321  if (in && !out) {
322  fgAODMCParticles = new TClonesArray("AliAODMCParticle",2*in->GetEntries());
323  fgAODMCParticles->SetName(AliAODMCParticle::StdBranchName());
324  oh->AddBranch("TClonesArray", &fgAODMCParticles);
325  out = static_cast<TClonesArray*>(eout->FindListObject(AliAODMCParticle::StdBranchName()));
326  }
327  if (in && out) {
328  if (out->GetEntries()>0) { // just checking if the deletion of previous event worked
329  AliFatal(Form("%s: Previous event not deleted. This should not happen!",GetName()));
330  }
331  out->AbsorbObjects(in);
332  if (fCutMC) {
333  for (Int_t i=0;i<out->GetEntriesFast();++i) {
334  AliAODMCParticle *mc = static_cast<AliAODMCParticle*>(in->At(i));
335  if ((mc==0)&&(i==0)) {
336  AliError(Form("%s: No MC info, skipping this event!",GetName()));
337  oh->SetFillAOD(kFALSE);
338  return;
339  }
340  if ((mc==0)||(TMath::Abs(mc->Y())>fYCutMC))
341  new ((*out)[i]) AliAODMCParticle;
342  }
343  }
344  }
345  }
346 
347  if (fDoCopyMCHeader) {
348  AliAODMCHeader *out = static_cast<AliAODMCHeader*>(eout->FindListObject(AliAODMCHeader::StdBranchName()));
349  AliAODMCHeader *in = static_cast<AliAODMCHeader*>(evin->FindListObject(AliAODMCHeader::StdBranchName()));
350  if (in && !out) {
351  fAODMcHeader = new AliAODMCHeader();
352  fAODMcHeader->SetName(AliAODMCHeader::StdBranchName());
353  oh->AddBranch("AliAODMCHeader",&fAODMcHeader);
354  out = static_cast<AliAODMCHeader*>(eout->FindListObject(AliAODMCHeader::StdBranchName()));
355  }
356  if (in && out) {
357  *out = *in;
358  if ((in->GetCrossSection()==0) && (fPyxsec>0)) {
359  out->SetCrossSection(fPyxsec);
360  out->SetTrials(fPytrials);
361  out->SetPtHard(fPypthardbin);
362  }
363  }
364  }
365 
366  if (gDebug>10) {
367  Int_t run = eout->GetRunNumber();
368  AliAODVertex *v=(AliAODVertex*)eout->GetVertices()->At(0);
369  Int_t vzn = v->GetNContributors();
370  Double_t vz = v->GetZ();
371  cout << "debug run " << run << " " << vzn << " " << vz << endl;
372  }
373 
374  fTrials = 0;
375  PostData(1, fOutputList);
376 }
377 
379 {
380  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
381  if (!tree) {
382  AliError(Form("%s: No current tree!",GetName()));
383  return kFALSE;
384  }
385 
386  Float_t xsection = 0;
387  Float_t trials = 0;
388  Int_t pthardbin = 0;
389 
390  TFile *curfile = tree->GetCurrentFile();
391  if (!curfile) {
392  AliError(Form("%s: No current file!",GetName()));
393  return kFALSE;
394  }
395 
396  TChain *chain = dynamic_cast<TChain*>(tree);
397  if (chain) tree = chain->GetTree();
398  Int_t nevents = tree->GetEntriesFast();
399 
400  Int_t slevel = gErrorIgnoreLevel;
401  gErrorIgnoreLevel = kFatal;
402  Bool_t res = PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
403  gErrorIgnoreLevel=slevel;
404 
405  if (res) {
406  cout << "AliAodSkimTask " << GetName() << " found xsec info: " << xsection << " " << trials << " " << pthardbin << " " << nevents << endl;
407  fPyxsec = xsection;
408  fPytrials = trials;
409  fPypthardbin = pthardbin;
410  }
411 
412  return res;
413 }
414 
416 {
417  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
418  if (man->GetAnalysisType()!=0)
419  return;
420  if (fHevs==0)
421  return;
422  cout << "AliAodSkimTask " << GetName() << " terminated with accepted fraction of events: " << fHevs->GetBinContent(2)/fHevs->GetEntries()
423  << " (" << fHevs->GetBinContent(2) << "/" << fHevs->GetEntries() << ")" << endl;
424 }
425 
426 Bool_t AliAodSkimTask::PythiaInfoFromFile(const char* currFile, Float_t &xsec, Float_t &trials, Int_t &pthard)
427 {
428  TString file(currFile);
429  xsec = 0;
430  trials = 1;
431 
432  if (file.Contains(".zip#")) {
433  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
434  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
435  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
436  file.Replace(pos+1,pos2-pos1,"");
437  } else {
438  // not an archive take the basename....
439  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
440  }
441  AliDebug(1,Form("File name: %s",file.Data()));
442 
443  // Get the pt hard bin
444  TString strPthard(file);
445 
446  strPthard.Remove(strPthard.Last('/'));
447  strPthard.Remove(strPthard.Last('/'));
448  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
449  strPthard.Remove(0,strPthard.Last('/')+1);
450  if (strPthard.IsDec()) {
451  pthard = strPthard.Atoi();
452  }
453  else {
454  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
455  pthard = -1;
456  }
457 
458  // 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
459  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
460 
461  if (!fxsec) {
462  // next trial fetch the histgram file
463  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
464  if (!fxsec) {
465  // not a severe condition but inciate that we have no information
466  return kFALSE;
467  } else {
468  // find the tlist we want to be independtent of the name so use the Tkey
469  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
470  if (!key) {
471  fxsec->Close();
472  return kFALSE;
473  }
474  TList *list = dynamic_cast<TList*>(key->ReadObj());
475  if (!list) {
476  fxsec->Close();
477  return kFALSE;
478  }
479  xsec = static_cast<TProfile*>(list->FindObject("h1Xsec"))->GetBinContent(1);
480  trials = static_cast<TH1F*>(list->FindObject("h1Trials"))->GetBinContent(1);
481  fxsec->Close();
482  }
483  } else { // no tree pyxsec.root
484  TTree *xtree = static_cast<TTree*>(fxsec->Get("Xsection"));
485  if (!xtree) {
486  fxsec->Close();
487  return kFALSE;
488  }
489  UInt_t ntrials = 0;
490  Double_t xsection = 0;
491  xtree->SetBranchAddress("xsection",&xsection);
492  xtree->SetBranchAddress("ntrials",&ntrials);
493  xtree->GetEntry(0);
494  trials = ntrials;
495  xsec = xsection;
496  fxsec->Close();
497  }
498  return kTRUE;
499 }
500 
501 const char *AliAodSkimTask::Str() const
502 {
503  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",
504  fClusMinE,
505  fCutMC,
506  fYCutMC,
507  fDoCopyHeader,
508  fDoCopyVZERO,
509  fDoCopyTZERO,
510  fDoCopyVertices,
511  fDoCopyTOF,
512  fDoCopyTracklets,
513  fDoCopyTracks,
514  fDoCopyTrigger,
515  fDoCopyPTrigger,
516  fDoCopyCells,
517  fDoCopyPCells,
518  fDoCopyClusters,
519  fDoCopyDiMuons,
520  fDoCopyTrdTracks,
521  fDoCopyV0s,
522  fDoCopyCascades,
523  fDoCopyZDC,
524  fDoCopyConv,
525  fDoCopyMC,
526  fDoCopyMCHeader);
527 }
528 
virtual Bool_t KeepTrack(AliAODTrack *t)
double Double_t
Definition: External.C:58
TSystem * gSystem
void Terminate(Option_t *option)
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()
Use to skim AOD files.
void UserCreateOutputObjects()
virtual void CleanTrack(AliAODTrack *t)
Bool_t PythiaInfoFromFile(const char *currFile, Float_t &xsec, Float_t &trials, Int_t &pthard)
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