macro dir
[u/mrichter/AliRoot.git] / PWG / EMCAL / AliEsdSkimTask.cxx
1 // $Id$
2 //
3 // Task to skim ESD files.
4 //
5 // Author: C.Loizides
6
7 #include "AliEsdSkimTask.h"
8 #include <TClonesArray.h>
9 #include <TFile.h>
10 #include <TTree.h>
11 #include "AliAnalysisManager.h"
12 #include "AliCentrality.h"
13 #include "AliESDEvent.h"
14 #include "AliESDFMD.h"
15 #include "AliESDtrackCuts.h"
16 #include "AliEsdTrackExt.h"
17 #include "AliEventplane.h"
18 #include "AliInputEventHandler.h"
19 #include "AliMultiplicity.h"
20 #include "AliPicoTrack.h"
21 #include "AliVEventHandler.h"
22
23 //_________________________________________________________________________________________________
24 AliEsdSkimTask::AliEsdSkimTask(const char *opt) :
25   AliAnalysisTaskSE(opt), fEvent(0), fTree(0), fCuts(0),
26   fDoZDC(1), fDoV0(1), fDoT0(1), fDoTPCv(1), fDoSPDv(1), fDoPriv(1),
27   fDoEmCs(1), fDoPCs(0), fDoEmT(1), fDoPT(0), fDoTracks(0), fDoFmd(1),
28   fDoMult(1), fDoTof(0), fDoPileup(1), fDoClus(0), fEmcNames(""), 
29   fDoMiniTracks(0), fTracks("Tracks"), fPhosClusOnly(0), fEmcalClusOnly(0),
30   fDoSaveBytes(0), fDoCent(1), fDoRP(1), fRemoveCP(0), fResetCov(1), 
31   fDoPicoTracks(0)
32 {
33   // Constructor.
34
35   if (!opt)
36     return;
37
38   fBranchNames = "ESD:AliESDHeader.,AliESDRun.";
39
40   DefineOutput(1, TTree::Class());
41 }
42
43 //_________________________________________________________________________________________________
44 void AliEsdSkimTask::UserExec(Option_t */*opt*/) 
45 {
46   // Process event.
47
48   AliESDEvent *esdin = dynamic_cast<AliESDEvent*>(InputEvent());
49   if (!esdin)
50     return;
51
52   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
53
54   fEvent->ResetStdContent();
55
56   TList* objsin  = esdin->GetList();
57   TList* objsout = fEvent->GetList();
58
59   AliESDHeader *header = dynamic_cast<AliESDHeader*>(objsout->FindObject("AliESDHeader"));
60   if (header) {
61     am->LoadBranch("AliESDHeader.");
62     if (!fDoSaveBytes) {
63       *header = *esdin->GetHeader();
64     } else {
65       AliESDHeader *hin = esdin->GetHeader();
66       header->SetTriggerMask(hin->GetTriggerMask());
67       header->SetOrbitNumber(hin->GetOrbitNumber());
68       header->SetTimeStamp(hin->GetTimeStamp());
69       header->SetEventType(hin->GetEventType());
70       header->SetEventSpecie(hin->GetEventSpecie());
71       header->SetEventNumberInFile(hin->GetEventNumberInFile());
72       header->SetBunchCrossNumber(hin->GetBunchCrossNumber());
73       header->SetPeriodNumber(hin->GetPeriodNumber());
74       header->SetTriggerCluster(hin->GetTriggerCluster());
75       header->SetL0TriggerInputs(hin->GetL0TriggerInputs());
76       header->SetL1TriggerInputs(hin->GetL1TriggerInputs());
77       header->SetL2TriggerInputs(hin->GetL2TriggerInputs());
78       for (Int_t i=0;i<24;++i) {
79         const char *name = hin->GetTriggerInputName(i,0);
80         if (name)
81           header->SetActiveTriggerInputs(name,i);
82       }
83       for (Int_t i=0;i<24;++i) {
84         const char *name = hin->GetTriggerInputName(i,1);
85         if (name)
86           header->SetActiveTriggerInputs(name,i+24);
87       }
88       for (Int_t i=0;i<12;++i) {
89         const char *name = hin->GetTriggerInputName(i,2);
90         if (name)
91           header->SetActiveTriggerInputs(name,i+48);
92       }
93     }
94     header->SetUniqueID(((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected());
95     if (am->GetInputEventHandler()) {
96       TTree *tree = am->GetInputEventHandler()->GetTree();
97       if (tree) {
98         TFile *cfile = tree->GetCurrentFile();
99         if (cfile) {
100           TString cname(cfile->GetName());
101           header->SetTitle(cname);
102         }
103       }
104     }
105   }
106   AliESDRun *run = dynamic_cast<AliESDRun*>(objsout->FindObject("AliESDRun"));
107   if (run) {
108     am->LoadBranch("AliESDRun.");
109     *run = *esdin->GetESDRun();
110   }
111   AliCentrality *cent = dynamic_cast<AliCentrality*>(objsout->FindObject("Centrality"));
112   if (cent) {
113     cent->Reset();
114     AliCentrality *centin = esdin->GetCentrality();
115     cent->SetQuality(centin->GetQuality());
116     cent->SetCentralityV0M(centin->GetCentralityPercentileUnchecked("V0M"));
117     cent->SetCentralityFMD(centin->GetCentralityPercentileUnchecked("FMD"));
118     cent->SetCentralityTRK(centin->GetCentralityPercentileUnchecked("TRK"));
119     cent->SetCentralityTKL(centin->GetCentralityPercentileUnchecked("TKL"));
120     cent->SetCentralityCL0(centin->GetCentralityPercentileUnchecked("CL0"));
121     cent->SetCentralityCL1(centin->GetCentralityPercentileUnchecked("CL1"));
122     cent->SetCentralityV0MvsFMD(centin->GetCentralityPercentileUnchecked("V0MvsFMD"));
123     cent->SetCentralityTKLvsV0M(centin->GetCentralityPercentileUnchecked("TKLvsV0M"));
124     cent->SetCentralityZEMvsZDC(centin->GetCentralityPercentileUnchecked("ZEMvsZDC"));
125   }
126   AliEventplane *ep = dynamic_cast<AliEventplane*>(objsout->FindObject("Eventplane"));
127   if (ep) {
128     ep->Reset();
129     AliEventplane *epin = esdin->GetEventplane();
130     if (!fDoSaveBytes) {
131       *ep = *epin;
132     } else {
133       if (epin->GetQVector()) {
134         ep->SetQVector(new TVector2(*epin->GetQVector()));
135         ep->SetEventplaneQ(epin->GetEventplane("Q"));
136         ep->SetQsub(new TVector2(*epin->GetQsub1()),new TVector2(*epin->GetQsub2()));
137         ep->SetQsubRes(epin->GetQsubRes());
138       }
139     }
140   }
141   if (fDoZDC) {
142     AliESDZDC *zdc = dynamic_cast<AliESDZDC*>(objsout->FindObject("AliESDZDC"));
143     if (zdc) {
144       am->LoadBranch("AliESDZDC.");
145       *zdc = *esdin->GetESDZDC();
146     }
147   }
148   if (fDoV0) {
149     AliESDVZERO *v0 = dynamic_cast<AliESDVZERO*>(objsout->FindObject("AliESDVZERO"));
150     if (v0) {
151       am->LoadBranch("AliESDVZERO.");
152       *v0 = *esdin->GetVZEROData();
153     }
154   }
155   if (fDoT0) {
156     AliESDTZERO *t0 = dynamic_cast<AliESDTZERO*>(objsout->FindObject("AliESDTZERO"));
157     if (t0) {
158       am->LoadBranch("AliESDTZERO.");
159       *t0 = *esdin->GetESDTZERO();
160     }
161   }
162   if (fDoTPCv) {
163     AliESDVertex *tpcv = dynamic_cast<AliESDVertex*>(objsout->FindObject("TPCVertex"));
164     if (tpcv) {
165       am->LoadBranch("TPCVertex.");
166       *tpcv = *esdin->GetPrimaryVertexTPC();
167     }
168   }
169   if (fDoSPDv) {
170     AliESDVertex *spdv = dynamic_cast<AliESDVertex*>(objsout->FindObject("SPDVertex"));
171     if (spdv) {
172       am->LoadBranch("SPDVertex.");
173       *spdv = *esdin->GetPrimaryVertexSPD();
174     }
175   }
176   if (fDoPriv) {
177     AliESDVertex *priv = dynamic_cast<AliESDVertex*>(objsout->FindObject("PrimaryVertex"));
178     if (priv) {
179       am->LoadBranch("PrimaryVertex.");
180       *priv = *esdin->GetPrimaryVertexTracks();
181     }
182   }
183   if (fDoEmCs) {
184     AliESDCaloCells *ecells = dynamic_cast<AliESDCaloCells*>(objsout->FindObject("EMCALCells"));
185     if (ecells) {
186       am->LoadBranch("EMCALCells.");
187       *ecells = *esdin->GetEMCALCells();
188     }
189   }
190   if (fDoPCs) {
191     AliESDCaloCells *pcells = dynamic_cast<AliESDCaloCells*>(objsout->FindObject("PHOSCells"));
192     if (pcells) {
193       am->LoadBranch("PHOSCells.");
194       *pcells = *esdin->GetPHOSCells();
195     }
196   }
197   if (fDoEmT) {
198     AliESDCaloTrigger *etrig = dynamic_cast<AliESDCaloTrigger*>(objsout->FindObject("EMCALTrigger"));
199     if (etrig) {
200       am->LoadBranch("EMCALTrigger.");
201       *etrig = *esdin->GetCaloTrigger("EMCAL");
202     }
203   }
204   if (fDoPT) {
205     AliESDCaloTrigger *ptrig = dynamic_cast<AliESDCaloTrigger*>(objsout->FindObject("PHOSTrigger"));
206     if (ptrig) {
207       am->LoadBranch("PHOSTrigger.");
208       *ptrig = *esdin->GetCaloTrigger("PHOS");
209     }
210   }
211   if (fDoFmd) {
212     AliESDFMD *fmd = dynamic_cast<AliESDFMD*>(objsout->FindObject("AliESDFMD"));
213     if (fmd) {
214       am->LoadBranch("AliESDFMD.");
215       if (!fDoSaveBytes) {
216         *fmd = *esdin->GetFMDData();
217       }
218     }
219   }
220   if (fDoMult) {
221     AliMultiplicity *mult = dynamic_cast<AliMultiplicity*>(objsout->FindObject("AliMultiplicity"));
222     if (mult) {
223       am->LoadBranch("AliMultiplicity.");
224       if (!fDoSaveBytes) {
225         *mult = *esdin->GetMultiplicity();
226       } else {
227         const AliMultiplicity *multin = esdin->GetMultiplicity();;
228         mult->SetFiredChips(0, multin->GetNumberOfFiredChips(0));
229         mult->SetFiredChips(1, multin->GetNumberOfFiredChips(1));
230         for (Int_t i=0; i<6; ++i) 
231           mult->SetITSClusters(i,multin->GetNumberOfITSClusters(i));
232       }
233     }
234   }
235   if (fDoTof) {
236     AliTOFHeader *tofh = dynamic_cast<AliTOFHeader*>(objsout->FindObject("AliTOFHeader"));
237     if (tofh) {
238       am->LoadBranch("AliTOFHeader.");
239       *tofh = *esdin->GetTOFHeader();
240     }
241   }
242   if (fDoPileup) {
243     TClonesArray *spup = dynamic_cast<TClonesArray*>(objsout->FindObject("SPDPileupVertices"));
244     if (spup) {
245       am->LoadBranch("SPDPileupVertices");
246       Int_t N = esdin->GetNumberOfPileupVerticesSPD();
247       for (Int_t i=0; i<N; ++i) {
248         const AliESDVertex *vtx = esdin->GetPileupVertexSPD(i);
249         if (vtx)
250           fEvent->AddPileupVertexSPD(vtx);
251       }
252     }
253     TClonesArray *tpup = dynamic_cast<TClonesArray*>(objsout->FindObject("TrkPileupVertices"));
254     if (tpup) {
255       am->LoadBranch("TrkPileupVertices");
256       Int_t N = esdin->GetNumberOfPileupVerticesTracks();
257       for (Int_t i=0; i<N; ++i) {
258         const AliESDVertex *vtx = esdin->GetPileupVertexTracks(i);
259         if (vtx)
260           fEvent->AddPileupVertexTracks(vtx);
261       }
262     }
263   }
264   if (fDoClus) {
265     TClonesArray *clus = dynamic_cast<TClonesArray*>(objsout->FindObject("CaloClusters"));
266     if (clus) {
267       am->LoadBranch("CaloClusters");
268       Int_t N = esdin->GetNumberOfCaloClusters();
269       for (Int_t i=0; i<N; ++i) {
270         AliESDCaloCluster *c = esdin->GetCaloCluster(i);
271         if (!c)
272           continue;
273         if (fPhosClusOnly && c->IsEMCAL())
274           continue;
275         if (fEmcalClusOnly && c->IsPHOS())
276           continue;
277         fEvent->AddCaloCluster(c);
278       }
279     }
280   }
281   TObjArray *namearr = fEmcNames.Tokenize(";");
282   if (namearr) {
283     for (Int_t i=0; i<namearr->GetEntries(); ++i) {
284       TString cname(namearr->At(i)->GetName());
285       if (cname.Length()<=0)
286         continue;
287       TClonesArray *arrin  = dynamic_cast<TClonesArray*>(objsin->FindObject(cname));
288       if (!arrin) {
289         AliFatal(Form("Can not find input clusters with name %s", cname.Data()));
290         return;
291       }
292       TClonesArray *arrout = dynamic_cast<TClonesArray*>(objsout->FindObject(cname));
293       if (!arrout) {
294         AliFatal(Form("Can not find output clusters with name %s", cname.Data()));
295         return;
296       }
297       arrout->Delete();
298       Double_t emin=0.1;
299       if (cname.Contains("FEE"))
300         emin = 1;
301       const Int_t N = arrin->GetEntries();
302       for (Int_t iC=0, nC=0; iC<N; ++iC) {
303         AliESDCaloCluster *c = dynamic_cast<AliESDCaloCluster*>(arrin->At(iC));
304         if (!c)
305           continue;
306         if (c->E()<emin)
307           continue;
308         AliESDCaloCluster *newCluster = new ((*arrout)[nC]) AliESDCaloCluster(*c);
309         newCluster->SetID(nC);
310         ++nC;
311       }
312     }
313     delete namearr;
314   }
315   if (fDoTracks) {
316     if (fTracks == "Tracks")
317       am->LoadBranch("Tracks");
318     TClonesArray *tracksin = dynamic_cast<TClonesArray*>(objsin->FindObject(fTracks));
319     if (!tracksin) {
320       AliFatal(Form("Can not find tracks with name %s", fTracks.Data()));
321       return;
322     }
323     TString trkoutname("Tracks");
324     if (fDoPicoTracks) 
325       trkoutname = "PicoTracks";
326     TClonesArray *tracksout = 0;
327     tracksout = dynamic_cast<TClonesArray*>(objsout->FindObject(trkoutname));
328     if (!tracksout) {
329       AliFatal(Form("Can not find tracks with name %s", trkoutname.Data()));
330       return;
331     }
332     tracksout->Delete();
333     TString classname(tracksin->GetClass()->GetName());
334     if (classname == "AliPicoTrack") {
335       if (!fDoPicoTracks) {
336         AliFatal(Form("Need to enable fDoPicoTracks when reading AliPicoTracks"));
337         return;
338       }
339       const Int_t Ntracks = tracksin->GetEntries();
340       for (Int_t iTracks = 0, nacc = 0; iTracks < Ntracks; ++iTracks) {
341         AliPicoTrack *track = dynamic_cast<AliPicoTrack*>(tracksin->At(iTracks));
342         if (!track)
343           continue;
344         new ((*tracksout)[nacc]) AliPicoTrack(*track);
345         ++nacc;
346       }
347     } else {
348       const Int_t Ntracks = tracksin->GetEntries();
349       for (Int_t iTracks = 0, nacc = 0; iTracks < Ntracks; ++iTracks) {
350         AliESDtrack *track = dynamic_cast<AliESDtrack*>(tracksin->At(iTracks));
351         if (!track)
352           continue;
353         if (fCuts) {
354           if (!fCuts->IsSelected(track))
355             continue;
356         }
357         if (fDoPicoTracks) {
358           AliEsdTrackExt newtrack(*track);
359           Double_t etaemc = 0;
360           Double_t phiemc = 0;
361           if (newtrack.IsEMCAL()) {
362             etaemc = newtrack.GetEmcEta();
363             phiemc = newtrack.GetEmcPhi();
364           }
365           new ((*tracksout)[nacc]) AliPicoTrack(newtrack.Pt(), newtrack.Eta(), newtrack.Phi(), 
366                                                 newtrack.Charge(), newtrack.GetLabel(), 
367                                                 etaemc, phiemc, newtrack.IsEMCAL());
368           ++nacc;
369         } else {
370           AliEsdTrackExt *newtrack = new ((*tracksout)[nacc]) AliEsdTrackExt(*track);
371           if (fDoMiniTracks) {
372             newtrack->MakeMiniTrack(0,fRemoveCP);
373             newtrack->ResetCovariance(fResetCov);
374           } else {
375             newtrack->DeleteParams();
376           }
377           newtrack->SetID(nacc);
378           ++nacc;
379         }
380       }
381     }
382   }
383   fTree->Fill();
384 }
385
386 //_________________________________________________________________________________________________
387 void AliEsdSkimTask::UserCreateOutputObjects() 
388 {
389   // Create output objects.
390
391   TFile *file = OpenFile(1);
392   fTree = new TTree("esdTree", "Tree with skimmed ESD objects");
393   file->SetCompressionLevel(5);
394   fTree->SetDirectory(file);
395   fTree->SetAutoFlush(-10*1024*1024);
396
397   fEvent = new AliESDEvent;
398   fEvent->CreateStdContent();
399   if (fDoCent) {
400     AliCentrality *cent = new AliCentrality;
401     fEvent->AddObject(cent);
402   }
403   if (fDoRP) {
404     AliEventplane *ep = new AliEventplane;
405     fEvent->AddObject(ep);
406   }
407   TObjArray *namearr = fEmcNames.Tokenize(";");
408   if (namearr) {
409     for (Int_t i=0; i<namearr->GetEntries(); ++i) {
410       TString cname(namearr->At(i)->GetName());
411       if (cname.Length()<=0)
412         continue;
413       TClonesArray *arr = new TClonesArray("AliESDCaloCluster",0);
414       arr->SetName(cname);
415       fEvent->AddObject(arr);
416     }
417     delete namearr;
418   }
419   if (fDoTracks) {
420     if (fDoPicoTracks) {
421       TClonesArray *arr = new TClonesArray("AliPicoTrack",0);
422       arr->SetName("PicoTracks");
423       fEvent->AddObject(arr);
424     } else {
425       TObject *obj = fEvent->GetList()->FindObject("Tracks");
426       fEvent->GetList()->Remove(obj);
427       TClonesArray *arr = new TClonesArray("AliEsdTrackExt",0);
428       arr->SetName("Tracks");
429       fEvent->AddObject(arr);
430     }
431   }
432   fEvent->GetStdContent();
433   fEvent->WriteToTree(fTree);
434   fTree->GetUserInfo()->Add(fEvent);
435   PostData(1,fTree);
436 }