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