]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliHFEreducedEventCreatorAOD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Use configuration script
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEreducedEventCreatorAOD.cxx
CommitLineData
3513afb7 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// Debug tree task
17// the tree is represented as reduced events
18//
19// Authors:
20// M.Fasel <M.Fasel@gsi.de>
21//
22#include <TArrayI.h>
23#include <TBits.h>
24#include <TFile.h>
25#include <TString.h>
26#include <TTree.h>
7e695443 27#include <TMath.h>
3513afb7 28
29#include "AliAnalysisManager.h"
30#include "AliAnalysisUtils.h"
31#include "AliAODEvent.h"
32#include "AliAODMCHeader.h"
33#include "AliAODMCParticle.h"
34#include "AliAODPid.h"
35#include "AliAODTrack.h"
36#include "AliAODTracklets.h"
37#include "AliAODVertex.h"
38#include "AliCentrality.h"
39#include "AliInputEventHandler.h"
40#include "AliHFEcuts.h"
41#include "AliHFEextraCuts.h"
ff8249bd 42#include "AliHFEmcQA.h"
3513afb7 43#include "AliHFEpidTPC.h"
44#include "AliHFEreducedEvent.h"
45#include "AliHFEreducedTrack.h"
46#include "AliHFEreducedMCParticle.h"
47#include "AliHFEsignalCuts.h"
ff8249bd 48#include "AliHFEV0taginfo.h"
3513afb7 49#include "AliLog.h"
50#include "AliMCEvent.h"
51#include "AliPIDResponse.h"
52#include "AliVEvent.h"
53#include "AliVTrack.h"
54#include "AliVCluster.h"
55#include "AliVVZERO.h"
56#include "AliVZDC.h"
57#include "TTreeStream.h"
58
59#include "AliHFEreducedEventCreatorAOD.h"
60
61ClassImp(AliHFEreducedEventCreatorAOD)
62
63AliHFEreducedEventCreatorAOD::AliHFEreducedEventCreatorAOD():
64 AliAnalysisTaskSE(),
65 fAODMCHeader(NULL),
66 fAODArrayMCInfo(NULL),
67 fHFEtree(NULL),
68 fAnalysisUtils(NULL),
69 fHFEevent(NULL),
70 fTrackCuts(NULL),
71 fExtraCuts(NULL),
72 fSignalCuts(NULL),
73 fTPCpid(NULL),
ff8249bd 74 fV0Tagger(NULL),
3513afb7 75 fEventNumber(0),
76 fNclustersTPC(70),
77 fNclustersTPCPID(0),
78 fNclustersITS(2),
7e695443 79 fNbOfTOFSigma(-1.0),
ff8249bd 80 fRemoveFirstEvent(kFALSE),
81 fFlagPileupEvents(kFALSE),
82 fSelectSignalOnly(kFALSE)
3513afb7 83{
84 //
85 // Default constructor
86 //
87}
88
89AliHFEreducedEventCreatorAOD::AliHFEreducedEventCreatorAOD(const char *name):
90 AliAnalysisTaskSE(name),
91 fAODMCHeader(NULL),
92 fAODArrayMCInfo(NULL),
93 fHFEtree(NULL),
94 fAnalysisUtils(NULL),
95 fHFEevent(NULL),
96 fTrackCuts(NULL),
97 fExtraCuts(NULL),
98 fSignalCuts(NULL),
99 fTPCpid(NULL),
ff8249bd 100 fV0Tagger(NULL),
3513afb7 101 fEventNumber(0),
102 fNclustersTPC(70),
103 fNclustersTPCPID(0),
104 fNclustersITS(2),
7e695443 105 fNbOfTOFSigma(-1.0),
ff8249bd 106 fRemoveFirstEvent(kFALSE),
107 fFlagPileupEvents(kFALSE),
108 fSelectSignalOnly(kFALSE)
3513afb7 109{
110 //
111 // Default constructor
112 //
113 fTPCpid = new AliHFEpidTPC("QAtpcPID");
114 fAnalysisUtils = new AliAnalysisUtils;
ff8249bd 115 fV0Tagger = new AliHFEV0taginfo("Tagger");
3513afb7 116 DefineOutput(1, TTree::Class());
117}
118
119AliHFEreducedEventCreatorAOD::~AliHFEreducedEventCreatorAOD(){
120 //
121 // Default destructor
122 //
123 if(fTPCpid) delete fTPCpid;
124 if(fAnalysisUtils) delete fAnalysisUtils;
ff8249bd 125 if(fV0Tagger) delete fV0Tagger;
3513afb7 126 if(fHFEevent) delete fHFEevent;
127}
128
129void AliHFEreducedEventCreatorAOD::UserCreateOutputObjects(){
130 //
131 // Create debug tree, signal cuts and track cuts
132 //
133
134 //printf("test\n");
135
136 // printf("testa\n");
137 fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
138 //printf("testb\n");
139
140 fTrackCuts = new AliHFEcuts("fTrackCuts", "Basic HFE track cuts");
141 fTrackCuts->CreateStandardCuts();
142 fTrackCuts->SetAOD();
143 // Track cuts
144 fTrackCuts->SetMinNClustersTPC(fNclustersTPC);
145 fTrackCuts->SetMinRatioTPCclusters(0);
146 fTrackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
147 fTrackCuts->SetMinNClustersTPCPID(fNclustersTPCPID);
148 fTrackCuts->SetMinNClustersITS(fNclustersITS);
149 // Event cuts
150 fTrackCuts->SetUseMixedVertex(kTRUE);
151 fTrackCuts->SetVertexRange(10.);
152 //printf("testa\n");
153 fTrackCuts->Initialize();
154 //printf("testb\n");
155
156 fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
157
158 fHFEevent = new AliHFEreducedEvent;
159 OpenFile(1);
160 fHFEtree = new TTree("HFEtree", "HFE event tree");
161 fHFEtree->Branch("HFEevent", "AliHFEreducedEvent", fHFEevent,128000,0);
162 PostData(1, fHFEtree);
163}
164
165void AliHFEreducedEventCreatorAOD::UserExec(Option_t *){
166 //
167 // User Exec: Fill debug Tree
168 //
169
170 // Get PID response
171 AliPIDResponse *pid = NULL;
172 AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
173 if(handler){
7e695443 174 // printf("testb\n");
3513afb7 175 pid = handler->GetPIDResponse();
176 } else {
177 AliError("No Handler");
178 }
179 if(!pid){
7e695443 180 // printf("testc\n");
3513afb7 181 AliError("No PID response");
182 return;
183 }
184 if(!fInputEvent) {
185 AliError("No Input event");
186 return;
187 }
188
189 if(fRemoveFirstEvent){
190 if(fAnalysisUtils->IsFirstEventInChunk(fInputEvent)) return;
191 }
192
193 AliAODTrack copyTrack;
194
195 // MC info
196 Bool_t mcthere = kTRUE;
197 AliAODEvent *aodE = dynamic_cast<AliAODEvent *>(fInputEvent);
198 if(!aodE){
7e695443 199 // printf("testd\n");
3513afb7 200 AliError("No AOD Event");
201 return;
202 }
203 fAODMCHeader = dynamic_cast<AliAODMCHeader *>(fInputEvent->FindListObject(AliAODMCHeader::StdBranchName()));
204 if(!fAODMCHeader){
205 mcthere = kFALSE;
7e695443 206 // return;
3513afb7 207 }
208 fAODArrayMCInfo = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
209 if(!fAODArrayMCInfo){
210 mcthere = kFALSE;
7e695443 211 // return;
3513afb7 212 }
213 else {
214 fSignalCuts->SetMCAODInfo(fAODArrayMCInfo);
215 fTrackCuts->SetMCEvent(aodE);
216 }
217
218 // Set for track cuts
219 if(mcthere){
220 // Check Monte-Carlo events for AODs
221 AliDebug(1, Form("Monte-Carlo Event: %p\n", fMCEvent));
222 if(fMCEvent){
223 AliDebug(1, Form("Available! Number of particles: %d\n", fMCEvent->GetNumberOfTracks()));
224 }
225 }
226 fTrackCuts->SetRecEvent(fInputEvent);
227
228 if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
229 AliDebug(1, "Event rejected by the event cuts\n");
230 return;
231 }
232 if(!fExtraCuts){
233 fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
234 }
235 fExtraCuts->SetRecEventInfo(fInputEvent);
236
ff8249bd 237 // Tag all v0s in current event
238 if(fV0Tagger){
239 fV0Tagger->Reset();
240 fV0Tagger->TagV0Tracks(fInputEvent);
241 }
3513afb7 242
243 // Make Reduced Event
244 //AliHFEreducedEvent hfeevent;
245 fHFEevent->~AliHFEreducedEvent();
246 new(fHFEevent)AliHFEreducedEvent();
247
248 // Get run number
249 fHFEevent->SetRunNumber(fInputEvent->GetRunNumber());
250
251 // Derive trigger
0a4bd5d7 252 AliDebug(1, "Get triggers\n");
3513afb7 253 UInt_t trigger = fInputHandler->IsEventSelected();
254 if(trigger & AliVEvent::kMB) fHFEevent->SetMBTrigger();
255 if(trigger & AliVEvent::kCentral) fHFEevent->SetCentralTrigger();
256 if(trigger & AliVEvent::kSemiCentral) fHFEevent->SetCentralTrigger();
257 if(trigger & AliVEvent::kEMCEJE) fHFEevent->SetEMCALTrigger();
258
259 // Get Primary Vertex
0a4bd5d7 260 AliDebug(1, "Get Primary Vertex\n");
3513afb7 261 const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
42d33aff 262 Double_t vtx[3]={0.,0.,0.};
0a4bd5d7 263 Double_t vcov[6];
264 Int_t ncontrib = -1;
265 if(vertex) {
266 AliDebug(1, "Found vertex\n");
267 vertex->GetXYZ(vtx);
268 ncontrib = vertex->GetNContributors();
269 vertex->GetCovarianceMatrix(vcov);
270 }
3513afb7 271 fHFEevent->SetVX(vtx[0]);
272 fHFEevent->SetVY(vtx[1]);
273 fHFEevent->SetVZ(vtx[2]);
3513afb7 274 fHFEevent->SetNContribVertex(ncontrib);
0a4bd5d7 275 fHFEevent->SetVertexResolution(TMath::Sqrt(TMath::Abs(vcov[5])));
ff8249bd 276 //fHFEevent->SetVertexDispersion(vertex->GetDispersion());
7e695443 277 // Get Primary Vertex from SPD
278 const AliVVertex *vertexSPD = aodE->GetPrimaryVertexSPD();
279 if(vertexSPD){
0a4bd5d7 280 AliDebug(1, "Found SPD vertex\n");
281 memset(vtx, 0, sizeof(Double_t) *3);
7e695443 282 vertexSPD->GetXYZ(vtx);
283 fHFEevent->SetVXSPD(vtx[0]);
284 fHFEevent->SetVYSPD(vtx[1]);
285 fHFEevent->SetVZSPD(vtx[2]);
286 fHFEevent->SetNContribVertexSPD(vertexSPD->GetNContributors());
287 memset(vcov, 0, sizeof(Double_t)*6);
0a4bd5d7 288 vertexSPD->GetCovarianceMatrix(vcov);
289 AliDebug(1, Form("Covariance Matrix vcov[5] %f\n",vcov[5]));
290 fHFEevent->SetVertexResolutionSPD(TMath::Sqrt(TMath::Abs(vcov[5])));
ff8249bd 291 //fHFEevent->SetVertexDispersionSPD(vertex->GetDispersion());
7e695443 292 }
3513afb7 293
294 // Get centrality
0a4bd5d7 295 AliDebug(1, "Centrality\n");
3513afb7 296 AliCentrality *hicent = fInputEvent->GetCentrality();
0a4bd5d7 297 if(hicent) fHFEevent->SetCentrality(
3513afb7 298 hicent->GetCentralityPercentile("V0M"),
299 hicent->GetCentralityPercentile("V0A"),
300 hicent->GetCentralityPercentile("V0C"),
301 hicent->GetCentralityPercentile("TKL"),
302 hicent->GetCentralityPercentile("TRK"),
4437a0d2 303 hicent->GetCentralityPercentile("ZNA"),
304 hicent->GetCentralityPercentile("ZNC"),
305 hicent->GetCentralityPercentile("CL0"),
306 hicent->GetCentralityPercentile("CL1"),
307 hicent->GetCentralityPercentile("CND")
3513afb7 308 );
309
310 // Get VZERO Information
0a4bd5d7 311 AliDebug(1, "VZERO info\n");
3513afb7 312 AliVVZERO *vzeroinfo = fInputEvent->GetVZEROData();
313 if(vzeroinfo) fHFEevent->SetV0Multiplicity(vzeroinfo->GetMTotV0A(), vzeroinfo->GetMTotV0C());
314
315 // Get ZDC Information
0a4bd5d7 316 AliDebug(1, "ZDC info\n");
3513afb7 317 AliVZDC *zdcinfo = fInputEvent->GetZDCData();
318 if(zdcinfo) fHFEevent->SetZDCEnergy(zdcinfo->GetZNAEnergy(), zdcinfo->GetZNCEnergy(), zdcinfo->GetZPAEnergy(), zdcinfo->GetZPCEnergy());
319
320 // Set SPD multiplicity
0a4bd5d7 321 AliDebug(1, "SPD multiplicity\n");
3513afb7 322 AliAODTracklets *tls = aodE->GetTracklets();
323 if(tls) fHFEevent->SetSPDMultiplicity(tls->GetNumberOfTracklets());
324
ff8249bd 325 // Flag Pileup
326 if(fFlagPileupEvents){
327 if(fAnalysisUtils->IsPileUpEvent(fInputEvent)) fHFEevent->SetPileupFlag();
328 }
329
3513afb7 330 // Look for kink mother
0a4bd5d7 331 AliDebug(1, "Vertices\n");
3513afb7 332 Int_t numberofvertices = aodE->GetNumberOfVertices();
333 Double_t listofmotherkink[numberofvertices];
334 Int_t numberofmotherkink = 0;
335 //printf("Number of vertices %d\n",numberofvertices);
336 for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
337 AliAODVertex *aodvertex = aodE->GetVertex(ivertex);
338 if(!aodvertex) continue;
339 //printf("Type %d\n",aodvertex->GetType());
340 if(aodvertex->GetType()==AliAODVertex::kKink) {
341 //printf("Find one kink\n");
342 AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
343 if(!mother) continue;
344 Int_t idmother = mother->GetID();
345 listofmotherkink[numberofmotherkink] = idmother;
346 //printf("ID %d\n",idmother);
347 numberofmotherkink++;
348 }
349 }
350 //printf("Number of kink mother in the events %d\n",numberofmotherkink);
351
352 //
353 // Loop on MC tracks only
354 //
355 AliAODMCParticle *mctrack(NULL);
356 // Monte-Carlo info
ff8249bd 357 //Int_t source(5);
3513afb7 358 if(mcthere){
ff8249bd 359 fHFEevent->SetVMC(fAODMCHeader->GetVtxX(),fAODMCHeader->GetVtxY(),fAODMCHeader->GetVtxZ());
0a4bd5d7 360 AliDebug(1, "Loop MC tracks\n");
3513afb7 361 for(Int_t itrack = 0; itrack < fAODArrayMCInfo->GetEntriesFast(); itrack++) {
362 mctrack = (AliAODMCParticle *)(fAODArrayMCInfo->At(itrack));
363 if(!mctrack) continue;
ff8249bd 364 if(TMath::Abs(mctrack->PdgCode()) != 11) continue; // in MC truth list store only electrons
365 if(fSelectSignalOnly && !fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) continue;
3513afb7 366 AliHFEreducedMCParticle hfemcpart;
4437a0d2 367 hfemcpart.SetSignal();
3513afb7 368 // Kinematics
369 hfemcpart.SetSignedPt(mctrack->Pt(), mctrack->Charge() > 0.);
370 hfemcpart.SetP(mctrack->P());
371 hfemcpart.SetEta(mctrack->Eta());
372 hfemcpart.SetPhi(mctrack->Phi());
373 hfemcpart.SetPdg(mctrack->GetPdgCode());
374
375 // Get Production Vertex in radial direction
376 hfemcpart.SetProductionVertex(mctrack->Xv(),mctrack->Yv(),mctrack->Zv());
377
378 // Get Mother PDG code of the particle
379 Int_t motherlabel = TMath::Abs(mctrack->GetMother());
380 if(motherlabel >= 0 && motherlabel < fAODArrayMCInfo->GetEntriesFast()){
381 AliAODMCParticle *mother = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(motherlabel));
ff8249bd 382 if(mother){
383 hfemcpart.SetMotherPdg(mother->GetPdgCode());
384 hfemcpart.SetMotherProductionVertex(mother->Xv(),mother->Yv(),mother->Zv());
385 }
3513afb7 386 }
387
388 // derive source
ff8249bd 389 /*
390 source = 7;
3513afb7 391 if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
ff8249bd 392 else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
393 else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
3513afb7 394 else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
ff8249bd 395 else if(fSignalCuts->IsJpsiElectron(mctrack)) source = 4;
396 else if(fSignalCuts->IsB2JpsiElectron(mctrack)) source = 5;
397 else if(fSignalCuts->IsKe3Electron(mctrack)) source = 6;
398 else source = 7;
3513afb7 399 hfemcpart.SetSource(source);
ff8249bd 400 */
401 hfemcpart.SetSource(static_cast<Int_t>(fSignalCuts->GetSignalSource(mctrack)));
7bdde22f 402 hfemcpart.SetElectronSource(fSignalCuts->GetMCQAObject()->GetElecSource(mctrack,kTRUE));
3513afb7 403
404 fHFEevent->AddMCParticle(&hfemcpart);
405 }
406 }
407
408 //
409 // Loop on reconstructed tracks
410 //
411 TArrayI arraytrack(fInputEvent->GetNumberOfTracks());
412 Int_t counterdc=0;
413
414 AliAODTrack *track = 0x0;
0a4bd5d7 415 AliDebug(1, "Loop reconstructed tracks\n");
3513afb7 416 for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
417 // Run track loop
418 track = dynamic_cast<AliAODTrack *>(fInputEvent->GetTrack(itrack));
419 if(!track) continue;
420 // Cut track (Only basic track cuts)
421 // printf("testv\n");
422 if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
423 //
424 //printf("testu\n");
425
426 // Kinematics
427 AliHFEreducedTrack hfetrack;
428 hfetrack.SetSignedPt(track->Pt(), track->Charge() > 0);
429 hfetrack.SetP(track->P());
430 hfetrack.SetEta(track->Eta());
431 hfetrack.SetPhi(track->Phi());
432 hfetrack.SetTPCmomentum(track->GetDetPid() ? track->GetDetPid()->GetTPCmomentum() : track->P());
433
434 // Track ID
435 hfetrack.SetTrackID(track->GetID());
436
437 // status
438 ULong_t status = track->GetStatus();
439 if((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) hfetrack.SetITSrefit();
440 if((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) hfetrack.SetTPCrefit();
441 if((status & AliVTrack::kTOFpid) == AliVTrack::kTOFpid) hfetrack.SetTOFpid();
442 //if((status & AliVTrack::kTOFmismatch) == AliVTrack::kTOFmismatch) hfetrack.SetTOFmismatch();
443 if(IsTOFmismatch(track, pid)) hfetrack.SetTOFmismatch(); // New version suggested by Pietro Antonioli
444 if(track->IsEMCAL()) hfetrack.SetEMCALpid();
ff8249bd 445
446 // fill counts of v0-identified particles
447 AliPID::EParticleType myv0pid = fV0Tagger->GetV0Info(track->GetID());
448 AliHFEreducedTrack::EV0PID_t v0pid = AliHFEreducedTrack::kV0undef;
449 if(myv0pid == AliPID::kElectron) v0pid = AliHFEreducedTrack::kV0electron;
450 else if(myv0pid == AliPID::kPion) v0pid = AliHFEreducedTrack::kV0pion;
451 else if(myv0pid == AliPID::kProton) v0pid = AliHFEreducedTrack::kV0proton;
452 hfetrack.SetV0PID(v0pid);
453
454 Double_t v0prodR = fV0Tagger->GetV0ProdR(track->GetID());
455 hfetrack.SetV0prodR(v0prodR);
456
3513afb7 457 // Filter
458 for(Int_t k=0; k<20; k++) {
459 Int_t u = 1<<k;
460 if((track->TestFilterBit(u))) {
461 hfetrack.SetFilterBit(k);
462 }
463 }
464
465 if(mcthere){
466 // Fill Monte-Carlo Information
467 Int_t label = TMath::Abs(track->GetLabel());
4437a0d2 468 if(label < fAODArrayMCInfo->GetEntriesFast())
3513afb7 469 mctrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(label));
4437a0d2 470 if(mctrack){
471 AliDebug(2, "Associated MC particle found\n");
472 if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) hfetrack.SetMCSignal();
473 // Kinematics
474 hfetrack.SetMCSignedPt(mctrack->Pt(),mctrack->Charge() > 0.);
475 hfetrack.SetMCP(mctrack->P());
476 hfetrack.SetMCEta(mctrack->Eta());
477 hfetrack.SetMCPhi(mctrack->Phi());
478 hfetrack.SetMCPDG(mctrack->GetPdgCode());
3513afb7 479
4437a0d2 480 // Get Production Vertex in radial direction
481 hfetrack.SetMCProdVtx(mctrack->Xv(),mctrack->Yv(),mctrack->Zv());
3513afb7 482
4437a0d2 483 // Get Mother PDG code of the particle
484 Int_t motherlabel = TMath::Abs(mctrack->GetMother());
485 if(motherlabel >= 0 && motherlabel < fAODArrayMCInfo->GetEntriesFast()){
486 AliAODMCParticle *mother = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(motherlabel));
ff8249bd 487 if(mother){
488 hfetrack.SetMCMotherPdg(mother->GetPdgCode());
489 hfetrack.SetMCMotherProdVtx(mother->Xv(),mother->Yv(),mother->Zv());
490 }
4437a0d2 491 }
3513afb7 492
4437a0d2 493 // derive source
ff8249bd 494 /*
495 source = 7;
496 if(fSignalCuts->IsCharmElectron(track)) source = 0;
497 else if(fSignalCuts->IsBeautyElectron(track)) source = 1;
498 else if(fSignalCuts->IsGammaElectron(track)) source = 2;
499 else if(fSignalCuts->IsNonHFElectron(track)) source = 3;
500 else if(fSignalCuts->IsJpsiElectron(track)) source = 4;
501 else if(fSignalCuts->IsB2JpsiElectron(track)) source = 5;
502 else if(fSignalCuts->IsKe3Electron(track)) source = 6;
503 else source = 7;
4437a0d2 504 hfetrack.SetMCSource(source);
ff8249bd 505 */
506 hfetrack.SetMCSource(static_cast<Int_t>(fSignalCuts->GetSignalSource(track)));
7bdde22f 507 hfetrack.SetMCElectronSource(fSignalCuts->GetMCQAObject()->GetElecSource(mctrack,kTRUE));
4437a0d2 508 } else {
509 AliDebug(2, "Associated MC particle not found\n");
510 }
3513afb7 511 }
512
513 // HFE DCA
514 Float_t dcaxy = -999.,
515 dcaz = -999.;
516 fExtraCuts->GetImpactParameters((AliVTrack *)track,dcaxy,dcaz);
517 hfetrack.SetDCA(dcaxy, dcaz);
4437a0d2 518 Double_t hfeImpactParam(-999.), hfeImpactParamResol(-999.);
519 fExtraCuts->GetHFEImpactParameters((AliVTrack *)track,hfeImpactParam,hfeImpactParamResol);
520 hfetrack.SetHFEImpactParam(hfeImpactParam,hfeImpactParamResol);
3513afb7 521
522 // Different number of clusters definitions
523 Int_t nclustersITS(track->GetITSNcls()),
524 nclustersTPC(track->GetTPCNcls()),
525 nclustersTPCall(track->GetTPCClusterMap().CountBits()),
526 nclustersTPCshared(0);
527 UChar_t nfindableTPC = track->GetTPCNclsF();
528 const TBits &sharedTPC = track->GetTPCSharedMap();
529 for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
530 hfetrack.SetITSnclusters(nclustersITS);
531 hfetrack.SetTPCnclusters(track->GetTPCNcls());
532 hfetrack.SetTRDnclusters(track->GetTRDncls());
533 hfetrack.SetTPCnclustersPID(track->GetTPCsignalN());
534 hfetrack.SetTPCcrossedRows(track->GetTPCNCrossedRows());
535 hfetrack.SetTPCnclustersAll(nclustersTPCall);
536 hfetrack.SetTPCsharedClusters(nclustersTPCshared);
537 hfetrack.SetTPCclusterRatio(nfindableTPC ? static_cast<Float_t>(nclustersTPC)/static_cast<Float_t>(nfindableTPC) : 0);
538 hfetrack.SetTPCclusterRatioAll(nfindableTPC ? static_cast<Float_t>(nclustersTPCall)/static_cast<Float_t>(nfindableTPC) : 0);
539 UChar_t itsPixel = track->GetITSClusterMap();
540 for(int ily = 0; ily < 6; ily++)
541 if(TESTBIT(itsPixel, ily)) hfetrack.SetITScluster(ily);
542
543 // TRD related quantities (Yvonne)
544 hfetrack.SetTRDntrackletsPID(track->GetTRDntrackletsPID());
3513afb7 545 AliAODPid* aodpid= track->GetDetPid();
ff8249bd 546 if(aodpid) hfetrack.SetTRDnslices(aodpid->GetTRDnSlices());
547 hfetrack.SetTRDchi2(track->GetTRDchi2());
548
549 Double_t* arraytrdsignals = 0x0;
550 if(aodpid) arraytrdsignals=aodpid->GetTRDslices();
3513afb7 551 Int_t nslicetemp=0;
552 for(Int_t iplane = 0; iplane < 6; iplane++){
553 nslicetemp=0;
554 for(Int_t b=(iplane*8);b<((iplane*8)+8);b++){
555 if(track->GetTRDntrackletsPID()>0){
58a496d1 556 if(arraytrdsignals && (arraytrdsignals[b]>0.001)) nslicetemp++;
3513afb7 557 }
558 }
559 if(nslicetemp > 0) hfetrack.SetTRDstatus(iplane);
560 }
561
562
563 // Kink
4437a0d2 564 if(fExtraCuts->IsKinkDaughter(track)) hfetrack.SetIsKinkDaughter();
3513afb7 565
566 // kink mother
3513afb7 567 for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
568 if(track->GetID() == listofmotherkink[kinkmother]) {
4437a0d2 569 hfetrack.SetIsKinkMother();
3513afb7 570 continue;
571 }
572 }
573
4437a0d2 574
3513afb7 575 // Double counted
576 Int_t id(track->GetID());
577 for(Int_t l=0; l < counterdc; l++){
578 Int_t iTrack2 = arraytrack.At(l);
579 if(iTrack2==id){
580 hfetrack.SetDoubleCounted();
581 break;
582 }
583 }
584 // Add the id at this place
585 arraytrack.AddAt(id,counterdc);
586 counterdc++;
587
588 // PID
589 hfetrack.SetTPCdEdx(track->GetDetPid() ? track->GetDetPid()->GetTPCsignal() : 0.);
590 hfetrack.SetTPCsigmaEl(pid->NumberOfSigmasTPC(track, AliPID::kElectron));
591 hfetrack.SetTOFsigmaEl(pid->NumberOfSigmasTOF(track, AliPID::kElectron));
592 hfetrack.SetTOFmismatchProbability(pid->GetTOFMismatchProbability(track));
4437a0d2 593 hfetrack.SetITSsigmaEl(pid->NumberOfSigmasITS(track, AliPID::kElectron));
3513afb7 594 // Eta correction
595 copyTrack.~AliAODTrack();
596 new(&copyTrack) AliAODTrack(*track);
597 if(fTPCpid->HasCentralityCorrection()) fTPCpid->ApplyCentralityCorrection(&copyTrack, static_cast<Double_t>(ncontrib),AliHFEpidObject::kAODanalysis);
598 if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(&copyTrack, AliHFEpidObject::kAODanalysis);
599 hfetrack.SetTPCsigmaElCorrected(pid->NumberOfSigmasTPC(&copyTrack, AliPID::kElectron));
600 hfetrack.SetTPCdEdxCorrected(copyTrack.GetDetPid() ? copyTrack.GetDetPid()->GetTPCsignal() : 0.);
601 if(track->IsEMCAL()){
602 // EMCAL cluster
603 Double_t emcalEnergyOverP = -1.,
604 showershape[4] = {0.,0.,0.,0.};
605 hfetrack.SetEMCALSigmaEl(pid->NumberOfSigmasEMCAL(track, AliPID::kElectron, emcalEnergyOverP, &showershape[0]));
606 hfetrack.SetEMCALEoverP(emcalEnergyOverP);
607 hfetrack.SetEMCALShowerShape(showershape);
608 }
609
7e695443 610 // If TOF cut
611 if(fNbOfTOFSigma>0.0){
0a4bd5d7 612 AliDebug(1, "TOF cut\n");
7e695443 613 if(!((status & AliVTrack::kTOFpid) == AliVTrack::kTOFpid)) continue;
614 if(TMath::Abs(pid->NumberOfSigmasTOF(track, AliPID::kElectron))> fNbOfTOFSigma) continue;
615 }
616
3513afb7 617 // Track finished, add NOW to the Event
618 fHFEevent->AddTrack(&hfetrack);
619 //printf("after\n");
620 }
621
622 // Fill the debug tree
623 AliInfo(Form("Number of tracks: %d\n", fHFEevent->GetNumberOfTracks()));
624 AliInfo(Form("Number of MC particles: %d\n", fHFEevent->GetNumberOfMCParticles()));
625 fHFEtree->Fill();
626
627 fEventNumber++;
628 PostData(1, fHFEtree);
629}
630
631void AliHFEreducedEventCreatorAOD::Terminate(Option_t *){
632 //
633 // Terminate
634 //
635 AliInfo("terminating...\n");
636
637}
638
639Bool_t AliHFEreducedEventCreatorAOD::IsTOFmismatch(const AliVTrack *const track, const AliPIDResponse *const pid) const {
640 //
641 // is TOF mismatch
642 //
643 Double_t probs[AliPID::kSPECIESC];
644 AliPIDResponse::EDetPidStatus status = pid->ComputeTOFProbability(track, AliPID::kSPECIESC, probs);
645 return status == AliPIDResponse::kDetMismatch;
646}
647