1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Analysis task to produce trees of lightweight events
17 // evgeny.kryshen@cern.ch
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"
42 #include "AliAODTracklets.h"
50 #include "TClonesArray.h"
51 ClassImp(AliAnalysisTaskCFTree)
53 //-----------------------------------------------------------------------------
54 AliAnalysisTaskCFTree::AliAnalysisTaskCFTree(const char* name) :
55 AliAnalysisTask(name,""),
62 fHybridConstrainedMask(0),
63 fTPConlyConstrainedMask(0),
66 fEventStatistics(0x0),
79 fSelectBit(AliVEvent::kMB),
82 fCentralityMethod("V0M"),
83 fTrackFilterBit(0xffffffff),
86 fSharedClusterCut(0.4),
88 fFoundFractionCut(0.8),
93 Info("AliAnalysisTaskCFTree","Calling Constructor");
94 DefineInput(0,TChain::Class());
95 DefineOutput(0,TList::Class());
96 DefineOutput(1,TTree::Class());
98 //-----------------------------------------------------------------------------
101 //-----------------------------------------------------------------------------
102 void AliAnalysisTaskCFTree::ConnectInputData(Option_t* /*option*/){
103 fInputHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
104 fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
105 fInputHandler->SetNeedField();
107 //-----------------------------------------------------------------------------
110 //-----------------------------------------------------------------------------
111 void AliAnalysisTaskCFTree::CreateOutputObjects(){
112 fListOfHistos = new TList();
113 fListOfHistos->SetOwner();
114 fEventStatistics = new TH1I("fEventStatistics","",10,0,10);
115 fEventStatistics->SetBit(TH1::kCanRebin);
117 fListOfHistos->Add(fEventStatistics);
119 if (fStoreTracks) fParticles = new TClonesArray("AliCFParticle",2000);
120 if (fStoreTracklets) fTracklets = new TClonesArray("AliCFParticle",100);
121 if (fStoreMuons) fMuons = new TClonesArray("AliCFParticle",100);
122 // create file-resident tree
123 TDirectory *owd = gDirectory;
125 fTree = new TTree("events","events");
127 fTree->Branch("cent",&fCentrality);
128 fTree->Branch("zvtx",&fZvtx);
129 fTree->Branch("field",&fField);
130 fTree->Branch("run",&fRunNumber);
131 fTree->Branch("period",&fPeriod);
132 fTree->Branch("orbit",&fOrbit);
133 fTree->Branch("bc",&fBc);
134 fTree->Branch("mask",&fSelectMask);
135 if (fStoreTracks) fTree->Branch("particles",&fParticles);
136 if (fStoreTracklets) fTree->Branch("tracklets",&fTracklets);
137 if (fStoreMuons) fTree->Branch("muons",&fMuons);
138 PostData(0,fListOfHistos);
141 //-----------------------------------------------------------------------------
144 //-----------------------------------------------------------------------------
145 void AliAnalysisTaskCFTree::Exec(Option_t *){
146 if (fParticles) fParticles->Clear();
147 if (fTracklets) fTracklets->Clear();
148 if (fMuons) fMuons->Clear();
150 fEventStatistics->Fill("before cuts",1);
152 AliVEvent* event = fInputHandler->GetEvent();
154 fEventStatistics->Fill("after event check",1);
156 fPIDResponse = fInputHandler->GetPIDResponse();
158 fSelectMask = fInputHandler->IsEventSelected();
159 if (!(fSelectMask & fSelectBit)) return;
160 fEventStatistics->Fill("after trigger selection",1);
161 fRunNumber = event->GetRunNumber();
162 fPeriod = event->GetPeriodNumber();
163 fOrbit = event->GetOrbitNumber();
164 fBc = event->GetBunchCrossNumber();
165 fField = event->GetMagneticField();
166 fCentrality = event->GetCentrality()->GetCentralityPercentile(fCentralityMethod);
168 if (fMode==0){ // Data analysis
169 const AliVVertex* vertex = event->GetPrimaryVertex();
171 fZvtx = vertex->GetZ();
172 TString name(vertex->GetName());
173 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex")) return; // Reject TPC only vertex
174 fEventStatistics->Fill("after rejection of TPC only vertex",1);
175 if (TMath::Abs(fZvtx) >= fZVertexCut || vertex->GetNContributors()<fnTracksVertex) return;
176 fEventStatistics->Fill("after vertex selection",1);
179 for (Int_t ipart=0;ipart<event->GetNumberOfTracks();ipart++){
180 AliVTrack* track = (AliVTrack*) event->GetTrack(ipart);
181 if (!track) continue;
182 UInt_t mask = GetFilterMap(track);
183 if (!(mask & fTrackFilterBit)) continue;
185 if (track->InheritsFrom("AliAODTrack")) AddTrack(track,mask,0);
186 else if (track->InheritsFrom("AliESDtrack")) {
187 if (mask) AddTrack(track,mask,1);
188 if (mask & fHybridConstrainedMask) AddTrack(track,mask,2);
189 if (mask & fTPConlyConstrainedMask) AddTrack(track,mask,3);
194 AliAODEvent* aod = (AliAODEvent*) event;
196 if (fStoreTracklets){
197 AliAODTracklets* tracklets = aod->GetTracklets();
198 Int_t nTracklets = tracklets->GetNumberOfTracklets();
199 for (Int_t i=0;i<nTracklets;i++){
200 Float_t phi = tracklets->GetPhi(i);
201 Float_t theta = tracklets->GetTheta(i);
202 Float_t dphi = tracklets->GetDeltaPhi(i);
203 new ((*fTracklets)[fTracklets->GetEntriesFast()]) AliCFParticle(dphi,-TMath::Log(TMath::Tan(theta/2)),phi,0,0);
208 for (Int_t iTrack = 0; iTrack < aod->GetNTracks(); iTrack++) {
209 AliAODTrack* track = aod->GetTrack(iTrack);
210 if (!track->IsMuonTrack()) continue;
211 Float_t pt = track->Pt();
212 Float_t eta = track->Eta();
213 Float_t phi = track->Phi();
214 Short_t charge = track->Charge();
215 Float_t dca = track->DCA();
216 Float_t chi2 = track->Chi2perNDF();
217 Float_t rabs = track->GetRAtAbsorberEnd();
218 Int_t mask = track->GetMatchTrigger();
219 if (rabs < 17.6 || rabs > 89.5) continue;
220 if (eta < -4 || eta > -2.5) continue;
221 AliCFParticle* part = new ((*fMuons)[fMuons->GetEntriesFast()]) AliCFParticle(pt,eta,phi,charge,mask,3);
230 else { // MC analysis
231 AliMCEvent* mcEvent = fMcHandler ? fMcHandler->MCEvent() : 0;
232 TClonesArray* mcTracks = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
233 if (!mcEvent && !mcTracks) { printf("No mc object found\n"); return; }
234 fEventStatistics->Fill("after mc objects check",1);
238 if (mcEvent) { // ESD
239 AliHeader* header = (AliHeader*) mcEvent->Header();
240 AliGenCocktailEventHeader* cocktailHeader = dynamic_cast<AliGenCocktailEventHeader*> (header->GenEventHeader());
241 AliGenEventHeader* mcHeader = dynamic_cast<AliGenEventHeader*> (cocktailHeader ? cocktailHeader->GetHeaders()->First() : header->GenEventHeader());
242 if (!mcHeader) { printf("mc header not found\n"); return; }
243 nProduced = mcEvent->GetNumberOfTracks();
244 nPrimGen = mcHeader->NProduced();
245 fZvtx = mcEvent->GetPrimaryVertex()->GetZ();
247 AliAODMCHeader* mcHeader = (AliAODMCHeader*) event->FindListObject(AliAODMCHeader::StdBranchName());
248 if (!mcHeader) { printf("AliAODMCHeader not found\n"); return; }
249 if (mcHeader->GetCocktailHeaders()) {
250 nProduced = mcTracks->GetEntriesFast();
251 AliGenEventHeader* header0 = mcHeader->GetCocktailHeader(0);
252 if (!header0) { printf("first header expected but not found\n"); return; }
253 nPrimGen = header0->NProduced();
254 } else nPrimGen = nProduced;
255 fZvtx = mcHeader->GetVtxZ();
257 fEventStatistics->Fill("after mc header check",1);
259 if (TMath::Abs(fZvtx)>fZVertexCut) return;
260 fEventStatistics->Fill("after MC vertex cut",1);
262 const AliVVertex* vertex = event->GetPrimaryVertex();
264 fEventStatistics->Fill("after check for vertex object",1);
265 TString name(vertex->GetName());
266 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex")) fZvtx=1000; // Reject TPC only vertex
267 fEventStatistics->Fill("after rejection of TPC only vertex",1);
268 if (vertex->GetNContributors()<fnTracksVertex) fZvtx=-1000;
269 fEventStatistics->Fill("after check for vertex contributors",1);
272 UInt_t* masks = new UInt_t[nProduced];
273 for (Int_t i=0;i<nProduced;i++) masks[i]=0;
275 // Loop over reconstructed tracks to set masks
276 for (Int_t ipart=0;ipart<event->GetNumberOfTracks();ipart++){
277 AliVTrack* part = (AliVTrack*) event->GetTrack(ipart);
279 Int_t label = TMath::Abs(part->GetLabel());
280 if (label>=nProduced) continue;
281 masks[label] |= GetFilterMap(part);
284 // Loop over mc tracks to be stored
285 for (Int_t ipart=0;ipart<nProduced;ipart++){
286 AliVParticle* part = 0;
287 Bool_t isPrimary = 0;
288 if (mcTracks) { // AOD analysis
289 AliAODMCParticle* particle = (AliAODMCParticle*) mcTracks->At(ipart);
290 if (!particle) continue;
291 // skip injected signals
292 AliAODMCParticle* mother = particle;
293 while (!mother->IsPrimary()) mother = (AliAODMCParticle*) mcTracks->At(mother->GetMother());
294 if (mother->GetLabel()>=nPrimGen) continue;
297 isPrimary = particle->IsPhysicalPrimary();
298 } else { // ESD analysis
299 AliMCParticle* particle = (AliMCParticle*) mcEvent->GetTrack(ipart);
300 if (!particle) continue;
301 // skip injected signals
302 AliMCParticle* mother = particle;
303 while (mother->GetMother()>=0) mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
304 if (mother->GetLabel()>=nPrimGen) continue;
307 isPrimary = mcEvent->IsPhysicalPrimary(ipart);
310 // store only primaries and all reconstructed (non-zero mask)
311 Int_t mask = masks[ipart];
312 if (isPrimary) mask |= (1 << 30);
320 AliAODEvent* aod = (AliAODEvent*) event;
322 if (fStoreTracklets) {
323 AliAODTracklets* tracklets = aod->GetTracklets();
324 Int_t nTracklets = tracklets->GetNumberOfTracklets();
325 for (Int_t i=0;i<nTracklets;i++){
326 Float_t phi = tracklets->GetPhi(i);
327 Float_t theta = tracklets->GetTheta(i);
328 Float_t dphi = tracklets->GetDeltaPhi(i);
329 AliCFParticle* tracklet = new ((*fTracklets)[fTracklets->GetEntriesFast()]) AliCFParticle(dphi,-TMath::Log(TMath::Tan(theta/2)),phi,0,0,4);
330 Int_t label1 = tracklets->GetLabel(i,0);
331 Int_t label2 = tracklets->GetLabel(i,1);
332 if (label1!=label2) continue;
333 AliAODMCParticle* particle = (AliAODMCParticle*) mcTracks->At(label1);
334 if (!particle) continue;
335 Short_t charge = particle->Charge();
336 Float_t ptMC = particle->Pt();
337 Float_t etaMC = particle->Eta();
338 Float_t phiMC = particle->Phi();
339 Float_t pdg = particle->PdgCode();
340 tracklet->SetCharge(charge);
341 tracklet->SetAt(ptMC,0);
342 tracklet->SetAt(etaMC,1);
343 tracklet->SetAt(phiMC,2);
344 tracklet->SetAt(pdg,3);
349 for (Int_t iTrack = 0; iTrack < aod->GetNTracks(); iTrack++) {
350 AliAODTrack* track = aod->GetTrack(iTrack);
351 if (!track->IsMuonTrack()) continue;
352 Float_t pt = track->Pt();
353 Float_t eta = track->Eta();
354 Float_t phi = track->Phi();
355 Short_t charge = track->Charge();
356 Float_t dca = track->DCA();
357 Float_t chi2 = track->Chi2perNDF();
358 Float_t rabs = track->GetRAtAbsorberEnd();
359 Int_t mask = track->GetMatchTrigger();
360 if (rabs < 17.6 || rabs > 89.5) continue;
361 if (eta < -4 || eta > -2.5) continue;
362 AliCFParticle* part = new ((*fMuons)[fMuons->GetEntriesFast()]) AliCFParticle(pt,eta,phi,charge,mask,3);
371 PostData(0,fListOfHistos);
374 //-----------------------------------------------------------------------------
377 //-----------------------------------------------------------------------------
378 UInt_t AliAnalysisTaskCFTree::GetFilterMap(AliVParticle* particle){
380 if (particle->InheritsFrom("AliAODTrack")) {
381 AliAODTrack* part = (AliAODTrack*) particle;
382 mask = part->GetFilterMap();
383 Double_t nCrossedRaws = part->GetTPCNCrossedRows();
384 Double_t nFindableClusters = part->GetTPCNclsF();
385 Double_t nSharedClusters = part->GetTPCnclsS();
386 Double_t nClusters = part->GetTPCncls();
387 Bool_t itsRefit = part->GetStatus() & AliVTrack::kITSrefit;
388 if (nCrossedRaws/nFindableClusters > fFoundFractionCut) mask |= (1 << 26);
389 if (nCrossedRaws>fCrossedRowsCut) mask |= (1 << 27);
390 if (itsRefit) mask |= (1 << 28);
391 if (nSharedClusters/nClusters<=fSharedClusterCut) mask |= (1 << 29);
392 if (part->GetLabel()<0) mask |= (1 << 31);
393 } else if (particle->InheritsFrom("AliESDtrack")){
394 AliESDtrack* part = (AliESDtrack*) particle;
395 if (!fTrackFilter) AliFatal("Track filter undefined");
396 mask |= fTrackFilter->IsSelected(part);
401 //-----------------------------------------------------------------------------
404 //-----------------------------------------------------------------------------
405 AliCFParticle* AliAnalysisTaskCFTree::AddTrack(AliVParticle* part, UInt_t mask, UInt_t flag){
407 // skip neutral mc particles
408 Char_t charge = part->Charge();
409 if (charge==0) return NULL;
412 Float_t pt=0,eta=0,phi=0;
413 if (flag==0 || flag==1){ // AOD, MC or Global ESD tracks
417 if (flag==1) mask &= (~fHybridConstrainedMask) & (~fTPConlyConstrainedMask);
419 else if (flag==2) { // Hybrid constrained tracks (ESD)
420 AliESDtrack* track = (AliESDtrack*) part;
421 const AliExternalTrackParam* param = track->GetConstrainedParam();
425 mask &= fHybridConstrainedMask;
427 else if (flag==3) { // TPC only constrained tracks (ESD)
428 AliESDtrack* track = (AliESDtrack*) part;
429 AliESDtrack tpcTrack;
430 if (!track->FillTPCOnlyTrack(tpcTrack)) return NULL;
431 AliExternalTrackParam param;
432 const AliESDVertex* vtxSPD = ((AliESDEvent*) fInputHandler->GetEvent())->GetPrimaryVertexSPD();
433 if (!tpcTrack.RelateToVertexTPC(vtxSPD,fField,1.e30,¶m)) return NULL;
437 mask &= fTPConlyConstrainedMask;
441 if (pt < fPtMin || TMath::Abs(eta) > fTrackEtaCut) return NULL;
444 AliCFParticle* cfpart = new ((*fParticles)[fParticles->GetEntriesFast()]) AliCFParticle(pt,eta,phi,charge,mask,3);
446 AliVTrack* track = dynamic_cast<AliVTrack*> (part);
447 if (!track) return cfpart;
448 Float_t ncl = track->GetTPCsignalN();
449 Float_t dedx = track->GetTPCsignalTunedOnData(); if (dedx<=0) dedx = track->GetTPCsignal();
451 if (track->GetStatus()&AliESDtrack::kTOFpid){
453 track->GetIntegratedTimes(tof);
454 beta = tof[0]/track->GetTOFsignal();
456 cfpart->SetAt(ncl,0);
457 cfpart->SetAt(dedx,1);
458 cfpart->SetAt(beta,2);