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