]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/Correlations/Base/AliAnalysisTaskCFTree.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Unification of AOD and ESD modes
[u/mrichter/AliRoot.git] / PWGCF / Correlations / Base / AliAnalysisTaskCFTree.cxx
CommitLineData
5fe7e785 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16// Analysis task to produce trees of lightweight events
17// evgeny.kryshen@cern.ch
18
19// aliroot
20#include "AliAnalysisTaskCFTree.h"
21#include "AliCFParticle.h"
22#include "AliAnalysisManager.h"
23#include "AliInputEventHandler.h"
24#include "AliESDEvent.h"
25#include "AliAODEvent.h"
26#include "AliMCEvent.h"
27#include "AliHeader.h"
28#include "AliESDHeader.h"
29#include "AliAODHeader.h"
30#include "AliAODMCHeader.h"
31#include "AliGenEventHeader.h"
32#include "AliGenCocktailEventHeader.h"
33#include "AliAODTrack.h"
34#include "AliESDtrack.h"
35#include "AliMCParticle.h"
36#include "AliAODMCParticle.h"
37#include "AliExternalTrackParam.h"
38#include "AliCentrality.h"
39#include "AliAnalysisFilter.h"
40#include "AliPIDResponse.h"
41#include "AliTPCPIDResponse.h"
4c66fa46 42#include "AliAODTracklets.h"
5fe7e785 43// root
44#include "TMath.h"
45#include "TFile.h"
46#include "TList.h"
47#include "TH1I.h"
48#include "TChain.h"
49#include "TTree.h"
50#include "TClonesArray.h"
51ClassImp(AliAnalysisTaskCFTree)
52
53//-----------------------------------------------------------------------------
54AliAnalysisTaskCFTree::AliAnalysisTaskCFTree(const char* name) :
55 AliAnalysisTask(name,""),
56 fDebug(0),
57 fMode(0),
4c66fa46 58 fIsAOD(1),
5fe7e785 59 fInputHandler(0x0),
60 fMcHandler(0x0),
61 fTrackFilter(0x0),
62 fHybridConstrainedMask(0),
63 fTPConlyConstrainedMask(0),
64 fPIDResponse(0x0),
65 fListOfHistos(0x0),
66 fEventStatistics(0x0),
67 fTree(0x0),
68 fParticles(0x0),
4c66fa46 69 fTracklets(0x0),
70 fMuons(0x0),
5fe7e785 71 fField(0),
72 fCentrality(0),
73 fZvtx(0),
74 fRunNumber(0),
75 fPeriod(0),
76 fOrbit(),
77 fBc(),
78 fSelectMask(0),
79 fSelectBit(AliVEvent::kMB),
80 fZVertexCut(10.),
81 fnTracksVertex(1),
82 fCentralityMethod("V0M"),
83 fTrackFilterBit(0xffffffff),
84 fTrackEtaCut(1.0),
85 fPtMin(0.15),
86 fSharedClusterCut(0.4),
87 fCrossedRowsCut(100),
4c66fa46 88 fFoundFractionCut(0.8),
c0ff46c8 89 fDphiCut(1.e9),
8e77b3fd 90 fStoreTracks(0),
4c66fa46 91 fStoreTracklets(0),
92 fStoreMuons(0)
5fe7e785 93{
94 Info("AliAnalysisTaskCFTree","Calling Constructor");
95 DefineInput(0,TChain::Class());
96 DefineOutput(0,TList::Class());
97 DefineOutput(1,TTree::Class());
98}
99//-----------------------------------------------------------------------------
100
101
102//-----------------------------------------------------------------------------
103void AliAnalysisTaskCFTree::ConnectInputData(Option_t* /*option*/){
104 fInputHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
105 fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
aafaf2f3 106 if (fInputHandler) fInputHandler->SetNeedField();
5fe7e785 107}
108//-----------------------------------------------------------------------------
109
110
111//-----------------------------------------------------------------------------
112void AliAnalysisTaskCFTree::CreateOutputObjects(){
113 fListOfHistos = new TList();
114 fListOfHistos->SetOwner();
115 fEventStatistics = new TH1I("fEventStatistics","",10,0,10);
116 fEventStatistics->SetBit(TH1::kCanRebin);
117
118 fListOfHistos->Add(fEventStatistics);
119
8e77b3fd 120 if (fStoreTracks) fParticles = new TClonesArray("AliCFParticle",2000);
4c66fa46 121 if (fStoreTracklets) fTracklets = new TClonesArray("AliCFParticle",100);
122 if (fStoreMuons) fMuons = new TClonesArray("AliCFParticle",100);
5fe7e785 123 // create file-resident tree
124 TDirectory *owd = gDirectory;
125 OpenFile(1);
126 fTree = new TTree("events","events");
127 owd->cd();
128 fTree->Branch("cent",&fCentrality);
129 fTree->Branch("zvtx",&fZvtx);
130 fTree->Branch("field",&fField);
131 fTree->Branch("run",&fRunNumber);
132 fTree->Branch("period",&fPeriod);
133 fTree->Branch("orbit",&fOrbit);
134 fTree->Branch("bc",&fBc);
135 fTree->Branch("mask",&fSelectMask);
8e77b3fd 136 if (fStoreTracks) fTree->Branch("particles",&fParticles);
137 if (fStoreTracklets) fTree->Branch("tracklets",&fTracklets);
138 if (fStoreMuons) fTree->Branch("muons",&fMuons);
5fe7e785 139 PostData(0,fListOfHistos);
140 PostData(1,fTree);
141}
142//-----------------------------------------------------------------------------
143
144
145//-----------------------------------------------------------------------------
146void AliAnalysisTaskCFTree::Exec(Option_t *){
8e77b3fd 147 if (fParticles) fParticles->Clear();
148 if (fTracklets) fTracklets->Clear();
149 if (fMuons) fMuons->Clear();
5fe7e785 150
151 fEventStatistics->Fill("before cuts",1);
152
153 AliVEvent* event = fInputHandler->GetEvent();
154 if (!event) return;
155 fEventStatistics->Fill("after event check",1);
156
157 fPIDResponse = fInputHandler->GetPIDResponse();
158
159 fSelectMask = fInputHandler->IsEventSelected();
160 if (!(fSelectMask & fSelectBit)) return;
161 fEventStatistics->Fill("after trigger selection",1);
162 fRunNumber = event->GetRunNumber();
163 fPeriod = event->GetPeriodNumber();
164 fOrbit = event->GetOrbitNumber();
165 fBc = event->GetBunchCrossNumber();
166 fField = event->GetMagneticField();
167 fCentrality = event->GetCentrality()->GetCentralityPercentile(fCentralityMethod);
168
169 if (fMode==0){ // Data analysis
170 const AliVVertex* vertex = event->GetPrimaryVertex();
171 if (!vertex) return;
172 fZvtx = vertex->GetZ();
173 TString name(vertex->GetName());
174 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex")) return; // Reject TPC only vertex
175 fEventStatistics->Fill("after rejection of TPC only vertex",1);
176 if (TMath::Abs(fZvtx) >= fZVertexCut || vertex->GetNContributors()<fnTracksVertex) return;
177 fEventStatistics->Fill("after vertex selection",1);
178
8e77b3fd 179 if (fParticles) {
180 for (Int_t ipart=0;ipart<event->GetNumberOfTracks();ipart++){
181 AliVTrack* track = (AliVTrack*) event->GetTrack(ipart);
182 if (!track) continue;
183 UInt_t mask = GetFilterMap(track);
184 if (!(mask & fTrackFilterBit)) continue;
5fe7e785 185
8e77b3fd 186 if (track->InheritsFrom("AliAODTrack")) AddTrack(track,mask,0);
187 else if (track->InheritsFrom("AliESDtrack")) {
188 if (mask) AddTrack(track,mask,1);
189 if (mask & fHybridConstrainedMask) AddTrack(track,mask,2);
190 if (mask & fTPConlyConstrainedMask) AddTrack(track,mask,3);
191 }
5fe7e785 192 }
193 }
4c66fa46 194 if (fIsAOD){
195 AliAODEvent* aod = (AliAODEvent*) event;
196
197 if (fStoreTracklets){
198 AliAODTracklets* tracklets = aod->GetTracklets();
199 Int_t nTracklets = tracklets->GetNumberOfTracklets();
200 for (Int_t i=0;i<nTracklets;i++){
201 Float_t phi = tracklets->GetPhi(i);
202 Float_t theta = tracklets->GetTheta(i);
203 Float_t dphi = tracklets->GetDeltaPhi(i);
c0ff46c8 204 if (TMath::Abs(dphi)>fDphiCut) continue;
4c66fa46 205 new ((*fTracklets)[fTracklets->GetEntriesFast()]) AliCFParticle(dphi,-TMath::Log(TMath::Tan(theta/2)),phi,0,0);
206 }
207 }
208
209 if (fStoreMuons){
210 for (Int_t iTrack = 0; iTrack < aod->GetNTracks(); iTrack++) {
211 AliAODTrack* track = aod->GetTrack(iTrack);
212 if (!track->IsMuonTrack()) continue;
213 Float_t pt = track->Pt();
214 Float_t eta = track->Eta();
215 Float_t phi = track->Phi();
216 Short_t charge = track->Charge();
217 Float_t dca = track->DCA();
218 Float_t chi2 = track->Chi2perNDF();
219 Float_t rabs = track->GetRAtAbsorberEnd();
220 Int_t mask = track->GetMatchTrigger();
221 if (rabs < 17.6 || rabs > 89.5) continue;
222 if (eta < -4 || eta > -2.5) continue;
223 AliCFParticle* part = new ((*fMuons)[fMuons->GetEntriesFast()]) AliCFParticle(pt,eta,phi,charge,mask,3);
224 part->SetAt(dca,0);
225 part->SetAt(chi2,1);
226 part->SetAt(rabs,2);
227 }
228 }
229 }
5fe7e785 230 }
4c66fa46 231
5fe7e785 232 else { // MC analysis
2bb5fec4 233 AliMCEvent* mcEvent = fMcHandler ? fMcHandler->MCEvent() : 0;
234 TClonesArray* mcTracks = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
5fe7e785 235 if (!mcEvent && !mcTracks) { printf("No mc object found\n"); return; }
236 fEventStatistics->Fill("after mc objects check",1);
237
2bb5fec4 238 Int_t nPrimGen = 0;
239 Int_t nProduced = 0;
240 if (mcEvent) { // ESD
5fe7e785 241 AliHeader* header = (AliHeader*) mcEvent->Header();
242 AliGenCocktailEventHeader* cocktailHeader = dynamic_cast<AliGenCocktailEventHeader*> (header->GenEventHeader());
243 AliGenEventHeader* mcHeader = dynamic_cast<AliGenEventHeader*> (cocktailHeader ? cocktailHeader->GetHeaders()->First() : header->GenEventHeader());
e428489e 244 if (!mcHeader) { printf("mc header not found\n"); return; }
5fe7e785 245 nProduced = mcEvent->GetNumberOfTracks();
246 nPrimGen = mcHeader->NProduced();
247 fZvtx = mcEvent->GetPrimaryVertex()->GetZ();
2bb5fec4 248 } else { // AOD
5fe7e785 249 AliAODMCHeader* mcHeader = (AliAODMCHeader*) event->FindListObject(AliAODMCHeader::StdBranchName());
5a502471 250 if (!mcHeader) { printf("AliAODMCHeader not found\n"); return; }
251 if (mcHeader->GetCocktailHeaders()) {
2bb5fec4 252 nProduced = mcTracks->GetEntriesFast();
5a502471 253 AliGenEventHeader* header0 = mcHeader->GetCocktailHeader(0);
254 if (!header0) { printf("first header expected but not found\n"); return; }
255 nPrimGen = header0->NProduced();
256 } else nPrimGen = nProduced;
5fe7e785 257 fZvtx = mcHeader->GetVtxZ();
2bb5fec4 258 }
5fe7e785 259 fEventStatistics->Fill("after mc header check",1);
260
261 if (TMath::Abs(fZvtx)>fZVertexCut) return;
262 fEventStatistics->Fill("after MC vertex cut",1);
263
264 const AliVVertex* vertex = event->GetPrimaryVertex();
265 if (!vertex) return;
266 fEventStatistics->Fill("after check for vertex object",1);
267 TString name(vertex->GetName());
268 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex")) fZvtx=1000; // Reject TPC only vertex
269 fEventStatistics->Fill("after rejection of TPC only vertex",1);
270 if (vertex->GetNContributors()<fnTracksVertex) fZvtx=-1000;
271 fEventStatistics->Fill("after check for vertex contributors",1);
272
8e77b3fd 273 if (fParticles) {
274 UInt_t* masks = new UInt_t[nProduced];
275 for (Int_t i=0;i<nProduced;i++) masks[i]=0;
276
277 // Loop over reconstructed tracks to set masks
278 for (Int_t ipart=0;ipart<event->GetNumberOfTracks();ipart++){
279 AliVTrack* part = (AliVTrack*) event->GetTrack(ipart);
280 if (!part) continue;
281 Int_t label = TMath::Abs(part->GetLabel());
282 if (label>=nProduced) continue;
283 masks[label] |= GetFilterMap(part);
284 }
285
286 // Loop over mc tracks to be stored
287 for (Int_t ipart=0;ipart<nProduced;ipart++){
288 AliVParticle* part = 0;
289 Bool_t isPrimary = 0;
290 if (mcTracks) { // AOD analysis
291 AliAODMCParticle* particle = (AliAODMCParticle*) mcTracks->At(ipart);
292 if (!particle) continue;
293 // skip injected signals
294 AliAODMCParticle* mother = particle;
295 while (!mother->IsPrimary()) mother = (AliAODMCParticle*) mcTracks->At(mother->GetMother());
296 if (mother->GetLabel()>=nPrimGen) continue;
297 part = particle;
298 // check for primary
299 isPrimary = particle->IsPhysicalPrimary();
300 } else { // ESD analysis
301 AliMCParticle* particle = (AliMCParticle*) mcEvent->GetTrack(ipart);
302 if (!particle) continue;
303 // skip injected signals
304 AliMCParticle* mother = particle;
305 while (mother->GetMother()>=0) mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
306 if (mother->GetLabel()>=nPrimGen) continue;
307 part = particle;
308 // check for primary
309 isPrimary = mcEvent->IsPhysicalPrimary(ipart);
310 }
311
312 // store only primaries and all reconstructed (non-zero mask)
313 Int_t mask = masks[ipart];
314 if (isPrimary) mask |= (1 << 30);
315 if (!mask) continue;
316 AddTrack(part,mask);
317 }
318 delete[] masks;
5fe7e785 319 }
320
8e77b3fd 321 if (fIsAOD){
322 AliAODEvent* aod = (AliAODEvent*) event;
323
aafaf2f3 324 if (fTracklets) {
8e77b3fd 325 AliAODTracklets* tracklets = aod->GetTracklets();
326 Int_t nTracklets = tracklets->GetNumberOfTracklets();
327 for (Int_t i=0;i<nTracklets;i++){
328 Float_t phi = tracklets->GetPhi(i);
329 Float_t theta = tracklets->GetTheta(i);
330 Float_t dphi = tracklets->GetDeltaPhi(i);
c0ff46c8 331 if (TMath::Abs(dphi)>fDphiCut) continue;
8e77b3fd 332 AliCFParticle* tracklet = new ((*fTracklets)[fTracklets->GetEntriesFast()]) AliCFParticle(dphi,-TMath::Log(TMath::Tan(theta/2)),phi,0,0,4);
333 Int_t label1 = tracklets->GetLabel(i,0);
334 Int_t label2 = tracklets->GetLabel(i,1);
335 if (label1!=label2) continue;
aafaf2f3 336 if (!mcTracks) continue;
8e77b3fd 337 AliAODMCParticle* particle = (AliAODMCParticle*) mcTracks->At(label1);
338 if (!particle) continue;
339 Short_t charge = particle->Charge();
340 Float_t ptMC = particle->Pt();
341 Float_t etaMC = particle->Eta();
342 Float_t phiMC = particle->Phi();
343 Float_t pdg = particle->PdgCode();
344 tracklet->SetCharge(charge);
345 tracklet->SetAt(ptMC,0);
346 tracklet->SetAt(etaMC,1);
347 tracklet->SetAt(phiMC,2);
348 tracklet->SetAt(pdg,3);
349 }
350 }
351
aafaf2f3 352 if (fMuons){
8e77b3fd 353 for (Int_t iTrack = 0; iTrack < aod->GetNTracks(); iTrack++) {
354 AliAODTrack* track = aod->GetTrack(iTrack);
355 if (!track->IsMuonTrack()) continue;
356 Float_t pt = track->Pt();
357 Float_t eta = track->Eta();
358 Float_t phi = track->Phi();
359 Short_t charge = track->Charge();
360 Float_t dca = track->DCA();
361 Float_t chi2 = track->Chi2perNDF();
362 Float_t rabs = track->GetRAtAbsorberEnd();
363 Int_t mask = track->GetMatchTrigger();
364 if (rabs < 17.6 || rabs > 89.5) continue;
365 if (eta < -4 || eta > -2.5) continue;
366 AliCFParticle* part = new ((*fMuons)[fMuons->GetEntriesFast()]) AliCFParticle(pt,eta,phi,charge,mask,3);
367 part->SetAt(dca,0);
368 part->SetAt(chi2,1);
369 part->SetAt(rabs,2);
370 }
5fe7e785 371 }
5fe7e785 372 }
5fe7e785 373 }
374 fTree->Fill();
375 PostData(0,fListOfHistos);
376 PostData(1,fTree);
377}
378//-----------------------------------------------------------------------------
379
380
381//-----------------------------------------------------------------------------
382UInt_t AliAnalysisTaskCFTree::GetFilterMap(AliVParticle* particle){
383 UInt_t mask = 0;
384 if (particle->InheritsFrom("AliAODTrack")) {
385 AliAODTrack* part = (AliAODTrack*) particle;
386 mask = part->GetFilterMap();
387 Double_t nCrossedRaws = part->GetTPCNCrossedRows();
388 Double_t nFindableClusters = part->GetTPCNclsF();
389 Double_t nSharedClusters = part->GetTPCnclsS();
390 Double_t nClusters = part->GetTPCncls();
391 Bool_t itsRefit = part->GetStatus() & AliVTrack::kITSrefit;
392 if (nCrossedRaws/nFindableClusters > fFoundFractionCut) mask |= (1 << 26);
393 if (nCrossedRaws>fCrossedRowsCut) mask |= (1 << 27);
394 if (itsRefit) mask |= (1 << 28);
395 if (nSharedClusters/nClusters<=fSharedClusterCut) mask |= (1 << 29);
396 if (part->GetLabel()<0) mask |= (1 << 31);
397 } else if (particle->InheritsFrom("AliESDtrack")){
398 AliESDtrack* part = (AliESDtrack*) particle;
399 if (!fTrackFilter) AliFatal("Track filter undefined");
400 mask |= fTrackFilter->IsSelected(part);
401 }
402
403 return mask;
404}
405//-----------------------------------------------------------------------------
406
407
408//-----------------------------------------------------------------------------
409AliCFParticle* AliAnalysisTaskCFTree::AddTrack(AliVParticle* part, UInt_t mask, UInt_t flag){
410
411 // skip neutral mc particles
412 Char_t charge = part->Charge();
413 if (charge==0) return NULL;
414
415 // set pt,eta,phi
416 Float_t pt=0,eta=0,phi=0;
417 if (flag==0 || flag==1){ // AOD, MC or Global ESD tracks
418 pt = part->Pt();
419 eta = part->Eta();
420 phi = part->Phi();
421 if (flag==1) mask &= (~fHybridConstrainedMask) & (~fTPConlyConstrainedMask);
422 }
423 else if (flag==2) { // Hybrid constrained tracks (ESD)
424 AliESDtrack* track = (AliESDtrack*) part;
425 const AliExternalTrackParam* param = track->GetConstrainedParam();
426 pt = param->Pt();
427 eta = param->Eta();
428 phi = param->Phi();
429 mask &= fHybridConstrainedMask;
430 }
431 else if (flag==3) { // TPC only constrained tracks (ESD)
432 AliESDtrack* track = (AliESDtrack*) part;
433 AliESDtrack tpcTrack;
434 if (!track->FillTPCOnlyTrack(tpcTrack)) return NULL;
435 AliExternalTrackParam param;
436 const AliESDVertex* vtxSPD = ((AliESDEvent*) fInputHandler->GetEvent())->GetPrimaryVertexSPD();
437 if (!tpcTrack.RelateToVertexTPC(vtxSPD,fField,1.e30,&param)) return NULL;
438 pt = param.Pt();
439 eta = param.Eta();
440 phi = param.Phi();
441 mask &= fTPConlyConstrainedMask;
442 }
443
444 // kinematic cuts
445 if (pt < fPtMin || TMath::Abs(eta) > fTrackEtaCut) return NULL;
446
e428489e 447
448 AliCFParticle* cfpart = new ((*fParticles)[fParticles->GetEntriesFast()]) AliCFParticle(pt,eta,phi,charge,mask,3);
449
450 AliVTrack* track = dynamic_cast<AliVTrack*> (part);
451 if (!track) return cfpart;
452 Float_t ncl = track->GetTPCsignalN();
453 Float_t dedx = track->GetTPCsignalTunedOnData(); if (dedx<=0) dedx = track->GetTPCsignal();
454 Float_t beta = -99;
455 if (track->GetStatus()&AliESDtrack::kTOFpid){
456 Double_t tof[5];
457 track->GetIntegratedTimes(tof);
458 beta = tof[0]/track->GetTOFsignal();
459 }
460 cfpart->SetAt(ncl,0);
461 cfpart->SetAt(dedx,1);
462 cfpart->SetAt(beta,2);
463 return cfpart;
5fe7e785 464}