AliPhysics  d565ceb (d565ceb)
AliAodSkimTask.cxx
Go to the documentation of this file.
1 #include <Riostream.h>
2 #include <TChain.h>
3 #include <TGraph.h>
4 #include <TH1F.h>
5 #include <TList.h>
6 #include <TMath.h>
7 #include <TProfile.h>
8 #include <TFile.h>
9 #include <TSystem.h>
10 #include <TKey.h>
11 #include <AliAODMCHeader.h>
12 #include <AliAODEvent.h>
13 #include <AliAODHandler.h>
14 #include <AliAODInputHandler.h>
15 #include <AliAODMCParticle.h>
16 #include <AliAnalysisManager.h>
17 #include "AliAodSkimTask.h"
18 #include <AliLog.h>
19 
20 using namespace std;
21 ClassImp(AliAodSkimTask)
22 
23 AliAodSkimTask::AliAodSkimTask() : AliAnalysisTaskSE(), fClusMinE(-1), fCutMC(1), fTrials(0), fPyxsec(0), fPytrials(0),
24  fPypthardbin(0), fAOD(0), fAODMcHeader(0), fOutputList(0)
25 {
26 }
27 
28 AliAodSkimTask::AliAodSkimTask(const char* name) : AliAnalysisTaskSE(name), fClusMinE(-1), fCutMC(1), fTrials(0), fPyxsec(0), fPytrials(0),
29  fPypthardbin(0), fAOD(0), fAODMcHeader(0), fOutputList(0)
30 {
31  DefineInput(0, TChain::Class());
32  DefineOutput(1, TList::Class());
33 }
34 
36 {
37  if (fOutputList) {
38  delete fOutputList;
39  }
40 }
41 
43 {
44 
45  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
46  AliAODHandler *oh = (AliAODHandler*)man->GetOutputEventHandler();
47  if (oh) {
48  TFile *fout = oh->GetTree()->GetCurrentFile();
49  fout->SetCompressionLevel(2);
50  }
51 
52  fOutputList = new TList;
53  fOutputList->SetOwner();
54  fHevs = new TH1F("hEvs","",2,-0.5,1.5);
55  fOutputList->Add(fHevs);
56  fHclus = new TH1F("hClus",";E (GeV)",200,0,100);
57  fOutputList->Add(fHclus);
58  PostData(1, fOutputList);
59 }
60 
62 {
63  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
64  if (!fAOD)
65  return;
66 
67  Bool_t store = kFALSE;
68 
69  if (fClusMinE>0) {
70  TClonesArray *cls = fAOD->GetCaloClusters();
71  for (Int_t i=0; i<cls->GetEntriesFast(); ++i) {
72  AliAODCaloCluster *clus = static_cast<AliAODCaloCluster*>(cls->At(i));
73  if (!clus->IsEMCAL())
74  continue;
75  Double_t e = clus->E();
76  fHclus->Fill(e);
77  if (e>fClusMinE) {
78  store = kTRUE;
79  }
80  }
81  } else {
82  store = kTRUE;
83  }
84 
85  if (!store) {
86  ++fTrials;
87  fHevs->Fill(0);
88  return;
89  }
90  fHevs->Fill(1);
91 
92  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
93  AliAODHandler *oh = (AliAODHandler*)man->GetOutputEventHandler();
94  if (oh)
95  oh->SetFillAOD(kFALSE);
96 
97  if (1) {
98  oh->SetFillAOD(kTRUE);
99  AliAODEvent *eout = dynamic_cast<AliAODEvent*>(oh->GetAOD());
100  AliAODEvent *evin = dynamic_cast<AliAODEvent*>(InputEvent());
101  TTree *tout = oh->GetTree();
102  if (tout) {
103  TList *lout = tout->GetUserInfo();
104  if (lout->FindObject("alirootVersion")==0) {
105  TList *lin = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetUserInfo();
106  for (Int_t jj=0;jj<lin->GetEntries()-1;++jj) {
107  lout->Add(lin->At(jj)->Clone(lin->At(jj)->GetName()));
108  }
109  }
110  }
111 
112  if (1) {
113  AliAODHeader *out = (AliAODHeader*)eout->GetHeader();
114  AliAODHeader *in = (AliAODHeader*)evin->GetHeader();
115  *out = *in;
116  out->SetUniqueID(fTrials);
117  }
118  if (1) {
119  AliAODVZERO *out = eout->GetVZEROData();
120  AliAODVZERO *in = evin->GetVZEROData();
121  *out = *in;
122  }
123  if (1) {
124  AliAODTZERO *out = eout->GetTZEROData();
125  AliAODTZERO *in = evin->GetTZEROData();
126  *out = *in;
127  }
128  if (1) {
129  TClonesArray *out = eout->GetVertices();
130  TClonesArray *in = evin->GetVertices();
131  for (Int_t i=0;i<in->GetEntriesFast();++i) {
132  AliAODVertex *v = static_cast<AliAODVertex*>(in->At(i));
133  new ((*out)[i]) AliAODVertex(*v);
134  }
135  }
136  if (1) {
137  AliTOFHeader *out = const_cast<AliTOFHeader*>(eout->GetTOFHeader());
138  const AliTOFHeader *in = evin->GetTOFHeader();
139  *out = *in;
140  }
141  if (1) {
142  TClonesArray *out = eout->GetTracks();
143  TClonesArray *in = evin->GetTracks();
144  for (Int_t i=0;i<in->GetEntriesFast();++i) {
145  AliAODTrack *t = static_cast<AliAODTrack*>(in->At(i));
146  new ((*out)[i]) AliAODTrack(*t);
147  }
148  }
149  if (1) {
150  AliAODCaloTrigger *out = eout->GetCaloTrigger("EMCAL");
151  AliAODCaloTrigger *in = evin->GetCaloTrigger("EMCAL");
152  *out = *in;
153  }
154  if (1) {
155  AliAODCaloCells *out = eout->GetEMCALCells();
156  AliAODCaloCells *in = evin->GetEMCALCells();
157  *out = *in;
158  }
159  if (1) {
160  TClonesArray *out = eout->GetCaloClusters();
161  TClonesArray *in = evin->GetCaloClusters();
162  for (Int_t i=0;i<in->GetEntriesFast();++i) {
163  AliAODCaloCluster *c = static_cast<AliAODCaloCluster*>(in->At(i));
164  new ((*out)[i]) AliAODCaloCluster(*c);
165  }
166  }
167 
168  if (1) {
169  TClonesArray *out = static_cast<TClonesArray*>(eout->FindListObject(AliAODMCParticle::StdBranchName()));
170  TClonesArray *in = static_cast<TClonesArray*>(evin->FindListObject(AliAODMCParticle::StdBranchName()));
171  if (in && !out) {
172  fgAODMCParticles = new TClonesArray("AliAODMCParticle",1000);
173  fgAODMCParticles->SetName(AliAODMCParticle::StdBranchName());
174  oh->AddBranch("TClonesArray", &fgAODMCParticles);
175  out = static_cast<TClonesArray*>(eout->FindListObject(AliAODMCParticle::StdBranchName()));
176  }
177  if (in && out) {
178  if (fCutMC) {
179  for (Int_t i=0;i<in->GetEntriesFast();++i) {
180  AliAODMCParticle *mc = static_cast<AliAODMCParticle*>(in->At(i));
181  if (TMath::Abs(mc->Y())>1.2)
182  new ((*out)[i]) AliAODMCParticle;
183  else
184  new ((*out)[i]) AliAODMCParticle(*mc);
185  }
186  }
187  }
188  }
189 
190  if (1) {
191  AliAODMCHeader *out = static_cast<AliAODMCHeader*>(eout->FindListObject(AliAODMCHeader::StdBranchName()));
192  AliAODMCHeader *in = static_cast<AliAODMCHeader*>(evin->FindListObject(AliAODMCHeader::StdBranchName()));
193  if (in && !out) {
194  fAODMcHeader = new AliAODMCHeader();
195  fAODMcHeader->SetName(AliAODMCHeader::StdBranchName());
196  oh->AddBranch("AliAODMCHeader",&fAODMcHeader);
197  out = static_cast<AliAODMCHeader*>(eout->FindListObject(AliAODMCHeader::StdBranchName()));
198  }
199  if (in && out) {
200  *out = *in;
201  if ((in->GetCrossSection()==0) && (fPyxsec>0)) {
202  out->SetCrossSection(fPyxsec);
203  out->SetTrials(fPytrials);
204  out->SetPtHard(fPypthardbin);
205  }
206  }
207  }
208  fTrials = 0;
209  PostData(1, fOutputList);
210  }
211 }
212 
214 {
215  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
216  if (!tree) {
217  AliError(Form("%s - UserNotify: No current tree!",GetName()));
218  return kFALSE;
219  }
220 
221  Float_t xsection = 0;
222  Float_t trials = 0;
223  Int_t pthardbin = 0;
224 
225  TFile *curfile = tree->GetCurrentFile();
226  if (!curfile) {
227  AliError(Form("%s - UserNotify: No current file!",GetName()));
228  return kFALSE;
229  }
230 
231  TChain *chain = dynamic_cast<TChain*>(tree);
232  if (chain) tree = chain->GetTree();
233  Int_t nevents = tree->GetEntriesFast();
234 
235  Int_t slevel = gErrorIgnoreLevel;
236  gErrorIgnoreLevel = kFatal;
237  Bool_t res = PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
238  gErrorIgnoreLevel=slevel;
239 
240  if (res) {
241  cout << "AliAodSkimTask " << GetName() << " found " << xsection << " " << trials << " " << pthardbin << " " << nevents << endl;
242  fPyxsec = xsection;
243  fPytrials = trials;
244  fPypthardbin = pthardbin;
245  }
246 
247  return res;
248 }
249 
251 {
252  cout << "AliAodSkimTask " << GetName() << " terminated with accepted fraction of events: " << fHevs->GetBinContent(2)/fHevs->GetEntries() << endl;
253 }
254 
255 Bool_t AliAodSkimTask::PythiaInfoFromFile(const char* currFile, Float_t &xsec, Float_t &trials, Int_t &pthard)
256 {
257  TString file(currFile);
258  xsec = 0;
259  trials = 1;
260 
261  if (file.Contains(".zip#")) {
262  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
263  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
264  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
265  file.Replace(pos+1,pos2-pos1,"");
266  } else {
267  // not an archive take the basename....
268  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
269  }
270  AliDebug(1,Form("File name: %s",file.Data()));
271 
272  // Get the pt hard bin
273  TString strPthard(file);
274 
275  strPthard.Remove(strPthard.Last('/'));
276  strPthard.Remove(strPthard.Last('/'));
277  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
278  strPthard.Remove(0,strPthard.Last('/')+1);
279  if (strPthard.IsDec()) {
280  pthard = strPthard.Atoi();
281  }
282  else {
283  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
284  pthard = -1;
285  }
286 
287  // 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
288  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
289 
290  if (!fxsec) {
291  // next trial fetch the histgram file
292  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
293  if (!fxsec) {
294  // not a severe condition but inciate that we have no information
295  return kFALSE;
296  } else {
297  // find the tlist we want to be independtent of the name so use the Tkey
298  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
299  if (!key) {
300  fxsec->Close();
301  return kFALSE;
302  }
303  TList *list = dynamic_cast<TList*>(key->ReadObj());
304  if (!list) {
305  fxsec->Close();
306  return kFALSE;
307  }
308  xsec = static_cast<TProfile*>(list->FindObject("h1Xsec"))->GetBinContent(1);
309  trials = static_cast<TH1F*>(list->FindObject("h1Trials"))->GetBinContent(1);
310  fxsec->Close();
311  }
312  } else { // no tree pyxsec.root
313  TTree *xtree = static_cast<TTree*>(fxsec->Get("Xsection"));
314  if (!xtree) {
315  fxsec->Close();
316  return kFALSE;
317  }
318  UInt_t ntrials = 0;
319  Double_t xsection = 0;
320  xtree->SetBranchAddress("xsection",&xsection);
321  xtree->SetBranchAddress("ntrials",&ntrials);
322  xtree->GetEntry(0);
323  trials = ntrials;
324  xsec = xsection;
325  fxsec->Close();
326  }
327  return kTRUE;
328 }
Int_t fPypthardbin
pythia trials
double Double_t
Definition: External.C:58
TSystem * gSystem
TH1F * fHevs
output list
TCanvas * c
Definition: TestFitELoss.C:172
Float_t fPyxsec
events seen since last acceptance
void Terminate(Option_t *option)
AliAODEvent * fAOD
pythia pthard bin
Float_t fPytrials
pythia xsection
Double_t fClusMinE
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
void UserExec(Option_t *option)
virtual ~AliAodSkimTask()
void UserCreateOutputObjects()
TList * fOutputList
MC header.
TH1F * fHclus
events processed/accepted
AliAODMCHeader * fAODMcHeader
input event
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