]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliHFEdebugTreeTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix of sigmaZ for crossing tracklets from Alex
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEdebugTreeTask.cxx
CommitLineData
fda75e2d 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//
18// Authors:
19// Markus Fasel <M.Fasel@gsi.de>
20//
21#include <TBits.h>
22#include <TString.h>
23
24#include "AliAnalysisManager.h"
25#include "AliCentrality.h"
26#include "AliESDEvent.h"
27#include "AliESDtrack.h"
28#include "AliESDtrackCuts.h"
29#include "AliESDVertex.h"
30#include "AliHFEcuts.h"
31#include "AliHFEsignalCuts.h"
32#include "AliInputEventHandler.h"
33#include "AliLog.h"
34#include "AliMCEvent.h"
35#include "AliMCEventHandler.h"
36#include "AliMCParticle.h"
d44e3c84 37#include "AliStack.h"
fda75e2d 38#include "AliPIDResponse.h"
11ff28c5 39#include "AliTrackReference.h"
fda75e2d 40#include "AliVEvent.h"
11ff28c5 41#include "AliHFEpidTPC.h"
a8ef1999 42#include "AliHFEpidTRD.h"
11ff28c5 43#include "AliHFEmcQA.h"
fda75e2d 44#include "TTreeStream.h"
45
46#include "AliHFEdebugTreeTask.h"
47
48ClassImp(AliHFEdebugTreeTask)
49
50AliHFEdebugTreeTask::AliHFEdebugTreeTask():
51 AliAnalysisTaskSE(),
52 fTrackCuts(NULL),
53 fSignalCuts(NULL),
a8ef1999 54 fTRDpid(NULL),
11ff28c5 55 fTPCpid(NULL),
56 fExtraCuts(NULL),
fda75e2d 57 fNclustersTPC(70),
58 fNclustersTPCPID(0),
59 fNclustersITS(2),
60 fFilename("HFEtree.root"),
11ff28c5 61 fDebugTree(NULL),
62 fNparents(-1)
fda75e2d 63{
64
65}
66
67AliHFEdebugTreeTask::AliHFEdebugTreeTask(const char *name):
68 AliAnalysisTaskSE(name),
69 fTrackCuts(NULL),
70 fSignalCuts(NULL),
a8ef1999 71 fTRDpid(NULL),
11ff28c5 72 fTPCpid(NULL),
73 fExtraCuts(NULL),
fda75e2d 74 fNclustersTPC(70),
75 fNclustersTPCPID(0),
76 fNclustersITS(2),
77 fFilename("HFEtree.root"),
11ff28c5 78 fDebugTree(NULL),
79 fNparents(-1)
fda75e2d 80{
11ff28c5 81 fTRDpid = new AliHFEpidTRD("QAtrdPID");
82 fTPCpid = new AliHFEpidTPC("QAtpcPID");
fda75e2d 83}
84
85AliHFEdebugTreeTask::~AliHFEdebugTreeTask(){
a8ef1999 86 if(fDebugTree) delete fDebugTree;
87 if(fTRDpid) delete fTRDpid;
11ff28c5 88 if(fTPCpid) delete fTPCpid;
fda75e2d 89}
90
91void AliHFEdebugTreeTask::UserCreateOutputObjects(){
92 //
93 // Create debug tree, signal cuts and track cuts
94 //
95 fDebugTree = new TTreeSRedirector(fFilename.Data());
96
97 fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
98
99 fTrackCuts = new AliHFEcuts("fTrackCuts", "Basic HFE track cuts");
100 fTrackCuts->CreateStandardCuts();
101 // Track cuts
102 fTrackCuts->SetMinNClustersTPC(fNclustersTPC);
103 fTrackCuts->SetMinRatioTPCclusters(0);
104 fTrackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
105 fTrackCuts->SetMinNClustersTPCPID(fNclustersTPCPID);
106 fTrackCuts->SetMinNClustersITS(fNclustersITS);
107 // Event cuts
108 fTrackCuts->SetUseMixedVertex(kTRUE);
109 fTrackCuts->SetVertexRange(10.);
110 fTrackCuts->Initialize();
a8ef1999 111
11ff28c5 112 fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
a8ef1999 113
fda75e2d 114}
115
116void AliHFEdebugTreeTask::UserExec(Option_t *){
117 //
118 // User Exec: Fill debug Tree
119 //
120
121 // Get PID response
122 AliPIDResponse *pid = NULL;
123 AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
124 if(handler){
125 pid = handler->GetPIDResponse();
126 } else {
127 AliError("No Handler");
128 }
129 if(!pid){
317ae21b 130 AliError("No PID response");
131 return;
132 }
133 if(!fInputEvent) {
134 AliError("No Input event");
fda75e2d 135 return;
136 }
11ff28c5 137
138 AliESDtrack copyTrack;
fda75e2d 139
140 fTrackCuts->SetRecEvent(fInputEvent);
141
142 // Cut event
143 if(fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5)){
0e30407a 144 AliDebug(1, "Event flagged as pileup\n");
fda75e2d 145 return;
146 }
147 if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
148 AliDebug(1, "Event rejected by the event cuts\n");
149 return;
150 }
151
152 // Get run number
153 Int_t run = fInputEvent->GetRunNumber();
154
155 // Derive trigger
156 UInt_t trigger = fInputHandler->IsEventSelected();
157 Bool_t isMBTrigger = trigger & AliVEvent::kMB;
158 Bool_t isCentralTrigger = trigger & AliVEvent::kCentral;
159 Bool_t isSemicentralTrigger = trigger & AliVEvent::kSemiCentral;
160 Bool_t isEMCALTrigger = trigger & AliVEvent::kEMCEJE;
161
162 // Check if MC information is available
163 Bool_t mcthere = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) != NULL;
164 if(mcthere){
165 fTrackCuts->SetMCEvent(fMCEvent);
166 fSignalCuts->SetMCEvent(fMCEvent);
167 }
168
d44e3c84 169 // MC get stack
170 AliStack* stack = 0x0;
171 if(mcthere){
172 stack = fMCEvent->Stack();
173 if(!stack) AliError("No Stack");
174 }
175
fda75e2d 176 // Get Primary Vertex
177 const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
11ff28c5 178 Double_t vtx[3];
179 vertex->GetXYZ(vtx);
180 Double_t ncontrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
fda75e2d 181
182 // Get centrality
183 Float_t centrality = -1.;
317ae21b 184 AliESDEvent *event = (dynamic_cast<AliESDEvent *>(fInputEvent));
185 if(!event) return;
186 TString beamtype = event->GetBeamType();
187 //printf("Beam type %s\n",(const char*)beamtype);
fda75e2d 188 if(!beamtype.CompareTo("Pb-Pb") || !beamtype.CompareTo("A-A")){
189 // Heavy ion run
190 AliDebug(1, "Heavy-Ion event\n");
191 AliCentrality *hicent = fInputEvent->GetCentrality();
192 centrality = hicent->GetCentralityPercentile("V0M");
193 }
194
11ff28c5 195 if(!fExtraCuts){
196 fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
197 }
198 fExtraCuts->SetRecEventInfo(event);
199
200 // Store event selection variables
201 (*fDebugTree) << "EventDebug"
202 << "Centrality=" << centrality
203 << "VertexZ=" << vtx[2]
204 << "NumberOfContributors=" << ncontrib
959ea9d8 205 << "run=" << run
11ff28c5 206 << "\n";
207
fda75e2d 208 // Common variables
209 Double_t charge, eta, phi, momentum, transversemomentum;
210 Int_t source;
211
212 // Monte-Carlo loop
213 if(mcthere){
214 AliMCParticle * mcpart;
215 for(Int_t itrk = 0; itrk < fMCEvent->GetNumberOfTracks(); itrk++){
216 mcpart = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(itrk));
217 if(!mcpart) continue;
218 if(!fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mcpart)) continue;
219 // Kinematics
220 charge = mcpart->Charge() > 0. ? 1. : -1.;
221 momentum = mcpart->P() * charge;
222 transversemomentum = mcpart->Pt() * charge;
223 eta = mcpart->Eta();
224 phi = mcpart->Phi();
225 Int_t pdg = mcpart->Particle()->GetPdgCode();
226
227 // Get Production Vertex in radial direction
228 Double_t vx = mcpart->Particle()->Vx(),
0e30407a 229 vy = mcpart->Particle()->Vy();
fda75e2d 230 Double_t productionVertex = TMath::Sqrt(vx*vx+vy*vy);
231
232 // Get Mother PDG code of the particle
233 Int_t motherPdg = 0;
234 Int_t motherlabel = mcpart->Particle()->GetFirstMother();
235 if(motherlabel >= 0 && motherlabel < fMCEvent->GetNumberOfTracks()){
236 AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherlabel));
237 if(mother) motherPdg = mother->Particle()->GetPdgCode();
238 }
239
240 // derive source
241 source = 5;
242 if(fSignalCuts->IsCharmElectron(mcpart)) source = 0;
243 else if(fSignalCuts->IsBeautyElectron(mcpart)) source = 1;
244 else if(fSignalCuts->IsGammaElectron(mcpart)) source = 2;
245 else if(fSignalCuts->IsNonHFElectron(mcpart)) source = 3;
246 else if(TMath::Abs(pdg) == 11) source = 4;
247 else source = 5;
248
11ff28c5 249 // Momemntum at the inner wall of the TPC
250 Double_t pTPC = 0., ptTPC = 0.;
251 AliTrackReference *ref = FindTrackReference(mcpart, 80, 270, AliTrackReference::kTPC);
252 if(ref){
253 pTPC = ref->P();
254 ptTPC = ref->Pt();
255 }
256
257
258
0e30407a 259 (*fDebugTree) << "MCDebug"
fda75e2d 260 << "centrality=" << centrality
261 << "MBtrigger=" << isMBTrigger
262 << "CentralTrigger=" << isCentralTrigger
263 << "SemicentralTrigger=" << isSemicentralTrigger
264 << "EMCALtrigger=" << isEMCALTrigger
265 << "run=" << run
266 << "p=" << momentum
267 << "pt=" << transversemomentum
268 << "eta=" << eta
269 << "phi=" << phi
270 << "pdg=" << pdg
271 << "ProductionVertex=" << productionVertex
272 << "motherPdg=" << motherPdg
cedf0381 273 << "source=" << source
fda75e2d 274 << "\n";
275 }
276 }
277
278 AliESDtrack *track;
11ff28c5 279 Double_t mcp, mcpt, mcptTPC, mcpTPC; // MC Variables added to the debug tree
fda75e2d 280 for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
281 // fill the tree
282 track = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(itrack));
283 if(!track) continue;
284 // Cut track (Only basic track cuts)
285 if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
286 // Debug streaming of PID-related quantities
0e30407a 287 copyTrack.~AliESDtrack();
11ff28c5 288 new(&copyTrack) AliESDtrack(*track);
289 if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(&copyTrack, AliHFEpidObject::kESDanalysis); // Apply Eta Correction on copy track
fda75e2d 290 Double_t nSigmaTOF = pid->NumberOfSigmasTOF(track, AliPID::kElectron);
11ff28c5 291 Double_t nSigmaTPC = pid->NumberOfSigmasTPC(&copyTrack, AliPID::kElectron);
292 //if(TMath::Abs(nSigmaTOF) > 5) continue;
fda75e2d 293 // we are not interested in tracks which are more than 5 sigma away from the electron hypothesis in either TOF or TPC
cedf0381 294 Double_t tPCdEdx = copyTrack.GetTPCsignal();
fda75e2d 295 // Signal, source and MCPID
296 Bool_t signal = kTRUE;
297 source = 5;
11ff28c5 298 mcp = mcpt = mcpTPC = mcptTPC = 0.;
299
300
301
302 Double_t bgcategory = 0.;
303 Int_t mArr = -1;
304 Int_t mesonID = -999;
305 Double_t xr[3]={-999,-999,-999};
306 Double_t eR=-999;
307 Double_t eZ=-999;
308 Double_t unique=-999;
309 Double_t mesonunique=-999;
310 Double_t mesonR=-999;
311 Double_t mesonZ=-999;
312 Double_t mesonMomPdg=-999;
313 Double_t mesonMomPt=-999;
314 Double_t mesonGMomPdg=-999;
315 Double_t mesonGGMomPdg=-999;
316 Double_t mesonPt = -999;
317 Double_t mceta = -999;
318 Double_t mcphi = -999;
319 Int_t mcpdg;
320
321 if(mcthere){
fda75e2d 322 // Signal
323 AliMCParticle *mctrack;
324 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))){
325 if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
326 }
327
328 // Source
329 if(fSignalCuts->IsCharmElectron(track)) source = 0;
330 else if(fSignalCuts->IsBeautyElectron(track)) source = 1;
331 else if(fSignalCuts->IsGammaElectron(track)) source = 2;
332 else if(fSignalCuts->IsNonHFElectron(track)) source = 3;
333 else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)) source = 4;
334 else source = 5;
11ff28c5 335
3024f297 336 if(!mctrack) continue;
337
11ff28c5 338 // Kinematics
339 mcpt = mctrack->Pt();
340 mcp = mctrack->P();
341 mceta = mctrack->Eta();
342 mcphi = mctrack->Phi();
343 mcpdg = mctrack->Particle()->GetPdgCode();
344
345 AliTrackReference *ref = FindTrackReference(mctrack, 80, 270, AliTrackReference::kTPC);
346 if(ref){
347 mcpTPC = ref->P();
348 mcptTPC = ref->Pt();
349 }
350
0e30407a 351
11ff28c5 352 TParticle *mctrack1 = mctrack->Particle();
353 mesonID=GetElecSourceMC(mctrack1);
354 if(mesonID==AliHFEmcQA::kGammaPi0 || mesonID==AliHFEmcQA::kPi0) mArr=0; //pion
355 else if(mesonID==AliHFEmcQA::kGammaEta || mesonID==AliHFEmcQA::kEta) mArr=1; //eta
356 else if(mesonID==AliHFEmcQA::kGammaOmega || mesonID==AliHFEmcQA::kOmega) mArr=2; //omega
357 else if(mesonID==AliHFEmcQA::kGammaPhi || mesonID==AliHFEmcQA::kPhi) mArr=3; //phi
358 else if(mesonID==AliHFEmcQA::kGammaEtaPrime || mesonID==AliHFEmcQA::kEtaPrime) mArr=4; //etaprime
359 else if(mesonID==AliHFEmcQA::kGammaRho0 || mesonID==AliHFEmcQA::kRho0) mArr=5; //rho
360
361 mctrack->XvYvZv(xr);
0e30407a 362
11ff28c5 363 eR= TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
364 eZ = xr[2];
365 TParticle *mctrackt = mctrack->Particle();
366 unique=mctrackt->GetUniqueID();
367
368 AliMCParticle *mctrackmother = NULL;
cedf0381 369 AliMCParticle *mctrackmother2 = NULL;
11ff28c5 370
371 if(!(mArr<0)){
0e30407a 372 if(mesonID>=AliHFEmcQA::kGammaPi0) { // conversion electron, be careful with the enum odering
373 Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
374 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
375 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
376 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
377 mesonPt = mctrackmother->Pt(); //meson pt
378 bgcategory = 1.;
379 mctrackmother->XvYvZv(xr);
380 mesonR = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
381 mesonZ = xr[2];
382
383 mctrackt = mctrackmother->Particle();
384 if(mctrackt){
cedf0381 385 mesonunique = mctrackt->GetUniqueID();
0e30407a 386 }
cedf0381 387
388 Int_t glabel2=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
389 if((mctrackmother2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel2)))){
390 mesonMomPdg=mctrackmother2->PdgCode();
391 mesonMomPt=mctrackmother2->Pt();
392 }
393
0e30407a 394 if(glabel>fMCEvent->GetNumberOfPrimaries()) {
395 bgcategory = 2.;
396 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
397 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
398 mesonMomPdg=mctrackmother->PdgCode();
399 mesonMomPt=mctrackmother->Pt();
400 if(TMath::Abs(mctrackmother->PdgCode())==310){
401 bgcategory = 3.;
402 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
403 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
404 mesonGMomPdg=mctrackmother->PdgCode();
405 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
406 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
407 mesonGGMomPdg=mctrackmother->PdgCode();
408 }
409 }
410 }
411 }
412 }
413 }
414 }
415 }
416 else{ // nonHFE except for the conversion electron
417 Int_t glabel=TMath::Abs(mctrack->GetMother());
418 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
419 mesonPt = mctrackmother->Pt(); //meson pt
420 bgcategory = -1.;
421 mctrackmother->XvYvZv(xr);
422 mesonR = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
423 mesonZ = xr[2];
424
425 mctrackt = mctrackmother->Particle();
426 if(mctrackt){
427 mesonunique = mctrackt->GetUniqueID();
428 }
429 if(glabel>fMCEvent->GetNumberOfPrimaries()) {
430 bgcategory = -2.;
431 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
432 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
433 mesonMomPdg=mctrackmother->PdgCode();
434 mesonMomPt=mctrackmother->Pt();
435 if(TMath::Abs(mctrackmother->PdgCode())==310){
436 bgcategory = -3.;
437 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
438 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
11ff28c5 439 mesonGMomPdg=mctrackmother->PdgCode();
0e30407a 440 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
441 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
442 mesonGGMomPdg=mctrackmother->PdgCode();
443 }
444 }
445 }
446 }
447 }
448 }
449 }
11ff28c5 450 }
451
452
453
fda75e2d 454 }
455 // Get V0 tag (if available)
456 Int_t v0pid = -1;
457 if(track->TestBit(BIT(14))) v0pid = AliPID::kElectron;
458 else if(track->TestBit(BIT(15))) v0pid = AliPID::kPion;
459 else if(track->TestBit(BIT(16))) v0pid = AliPID::kProton;
460 // Kinematics
461 charge = track->Charge() > 0 ? 1. : -1.;
462 eta = track->Eta();
463 phi = track->Phi();
464 momentum = track->P() * charge;
465 transversemomentum = track->Pt() * charge;
11ff28c5 466 Double_t momentumTPC = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->P() : 0.;
467 Double_t transversemomentumTPC = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->Pt() : 0.;
fda75e2d 468 // ITS number of clusters
469 UChar_t nclustersITS = track->GetITSclusters(NULL);
470 Double_t chi2matching = track->GetChi2TPCConstrainedVsGlobal(dynamic_cast<const AliESDVertex *>(vertex));
471 Double_t chi2PerClusterITS = 0.0;
472 if (nclustersITS != 0) chi2PerClusterITS = track->GetITSchi2() / Float_t(nclustersITS);
473 // TPC number of clusters (different definitions)
76d0b522 474 UChar_t nclustersTPCfit = track->GetTPCNcls();
475 UChar_t nclustersTPCall = 0;
476 const TBits &clusterTPC = track->GetTPCClusterMap();
477 nclustersTPCall = clusterTPC.CountBits();
fda75e2d 478 UChar_t nclustersTPCPID = track->GetTPCsignalN();
479 UChar_t nfindableTPC = track->GetTPCNclsF();
76d0b522 480 Double_t clusterRatioTPCfit = 0.0;
481 if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPCfit = static_cast<Double_t>(nclustersTPCfit)/static_cast<Double_t>(nfindableTPC);
482 Double_t clusterRatioTPCall = 0.0;
483 if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPCall = static_cast<Double_t>(nclustersTPCall)/static_cast<Double_t>(nfindableTPC);
fda75e2d 484 UChar_t nclustersTPCshared = 0;
485 Float_t ncrossedRowsTPC = track->GetTPCCrossedRows();
486 const TBits &sharedTPC = track->GetTPCSharedMap();
487 for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
488 // TRD clusters and tracklets
489 UChar_t nclustersTRD = track->GetTRDncls();
a8ef1999 490 UChar_t ntrackletsTRDPID = track->GetTRDntrackletsPID();
fda75e2d 491 // ITS and TRD acceptance maps
91e50e2b 492 UChar_t hasClusterITS[6], statusITS[6], hasTrackletTRD[6];
fda75e2d 493 UChar_t itsPixel = track->GetITSClusterMap();
cedf0381 494 for(Int_t icl = 0; icl < 6; icl++){
495 hasClusterITS[icl] = TESTBIT(itsPixel, icl) ? 1 : 0;
496 if(CheckITSstatus(track, icl)) statusITS[icl] = 1;
497 else statusITS[icl] = 0;
498 }
a8ef1999 499 Double_t trddEdxSum[6];
500 for(Int_t a=0;a<6;a++) { trddEdxSum[a]= 0.;}
fda75e2d 501 for(Int_t itl = 0; itl < 6; itl++){
0e30407a 502 Int_t nSliceNonZero = 0;
11ff28c5 503 trddEdxSum[itl] = track->GetTRDslice(itl, 0); // in new reconstruction slice 0 contains the total charge
fda75e2d 504 for(Int_t islice = 0; islice < 8; islice++){
505 if(track->GetTRDslice(itl, islice) > 0.001) nSliceNonZero++;
506 }
507 hasTrackletTRD[itl] = nSliceNonZero ? 1 : 0;
508 }
509 // TRD PID
510 Double_t pidprobs[5];
511 track->GetTRDpid(pidprobs);
512 Double_t likeEleTRD = pidprobs[0];
513 Double_t likeEleTRDn = likeEleTRD/(likeEleTRD + pidprobs[2]);
a8ef1999 514 Double_t trdtruncmean1 = fTRDpid->GetTRDSignalV1(track, 0.6);
515 Double_t trdtruncmean2 = fTRDpid->GetTRDSignalV2(track, 0.6);
516
fda75e2d 517 // DCA
518 Float_t b[2] = {0.,0.};
519 Float_t bCov[3] = {0.,0.,0.};
520 track->GetImpactParameters(b,bCov);
521 Double_t dca = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); // impact parameter space
522 Double_t dcaSR=0, dcaSZ=0;
523 if(bCov[0]>0) dcaSR = b[0]/TMath::Sqrt(bCov[0]); // normalised impact parameter xy
524 if(bCov[2]>0) dcaSZ = b[1]/TMath::Sqrt(bCov[2]); // normalised impact parameter z
525 Double_t dcaS = AliESDtrackCuts::GetSigmaToVertex(track); // n_sigma
11ff28c5 526
527 // HFE DCA
528 Double_t hfeb[2] = {-99.,-99.};
529 Double_t hfebCov[3] = {-999.,-999.,-999.};
530 fExtraCuts->GetHFEImpactParameters(track, hfeb, hfebCov);
531
0e30407a 532 Double_t tofdx= -999.0;
533 Double_t tofdz= -999.0;
534 tofdx=track->GetTOFsignalDx();
535 tofdz=track->GetTOFsignalDz();
536
d44e3c84 537 // TOF track status
538 UInt_t status = 0;
539 status = track->GetStatus();
540 Bool_t hasTOFout = status&AliESDtrack::kTOFout;
541 Bool_t hasTOFtime = status&AliESDtrack::kTIME;
542 Bool_t hasTOFpid = status&AliESDtrack::kTOFpid;
543 Bool_t hasgoodTOF = kFALSE;
544 if (hasTOFout && hasTOFtime && hasTOFpid) hasgoodTOF = kTRUE;
545
546 // TRD track status
547 Bool_t hasTRDin = status&AliESDtrack::kTRDin;
548
549
550 // TOF mismatch (particle spectra group)
551 Int_t mismatchlevel=0; // default value; in data always 0
552 if(mcthere){
553 Int_t tofLabel[3];
554 track->GetTOFLabel(tofLabel);
555 if(TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0) mismatchlevel=1;
556 TParticle *matchedTrack = stack->Particle(TMath::Abs(tofLabel[0]));
557 if(TMath::Abs(matchedTrack->GetFirstMother()) == TMath::Abs(track->GetLabel()))
558 {
559 if(mismatchlevel==1) mismatchlevel=3;
560 else mismatchlevel=2;
561 }
562 }
563
564
fda75e2d 565 // Fill Tree
566 (*fDebugTree) << "PIDdebug"
567 << "centrality=" << centrality
568 << "MBtrigger=" << isMBTrigger
569 << "CentralTrigger=" << isCentralTrigger
570 << "SemicentralTrigger=" << isSemicentralTrigger
571 << "EMCALtrigger=" << isEMCALTrigger
572 << "signal=" << signal
573 << "source=" << source
574 << "v0pid=" << v0pid
575 << "run=" << run
576 << "p=" << momentum
11ff28c5 577 << "ptpc=" << momentumTPC
fda75e2d 578 << "pt=" << transversemomentum
11ff28c5 579 << "pttpc=" << transversemomentumTPC
580 << "mcp=" << mcp
581 << "mcpt=" << mcpt
582 << "mcpTPC=" << mcpTPC
583 << "mcptTPC=" << mcptTPC
584 << "mceta=" << mceta
585 << "mcphi=" << mcphi
586 << "mcpdg=" << mcpdg
fda75e2d 587 << "eta=" << eta
588 << "phi=" << phi
a8ef1999 589 << "ntracklets=" << ntrackletsTRDPID
76d0b522 590 << "nclustersTPC=" << nclustersTPCfit
591 << "nclustersTPCall=" << nclustersTPCall
fda75e2d 592 << "nclustersTPCPID=" << nclustersTPCPID
593 << "nclustersTPCshared=" << nclustersTPCshared
594 << "ncrossedRowsTPC=" << ncrossedRowsTPC
76d0b522 595 << "clusterRatioTPC=" << clusterRatioTPCfit
596 << "clusterRatioTPCall=" << clusterRatioTPCall
fda75e2d 597 << "nclustersITS=" << nclustersITS
598 << "nclusters=" << nclustersTRD
599 << "chi2matching=" << chi2matching
cedf0381 600 << "chi2PerClusterITS=" << chi2PerClusterITS
fda75e2d 601 << "its0=" << hasClusterITS[0]
602 << "its1=" << hasClusterITS[1]
603 << "its2=" << hasClusterITS[2]
604 << "its3=" << hasClusterITS[3]
605 << "its4=" << hasClusterITS[4]
606 << "its5=" << hasClusterITS[5]
cedf0381 607 << "statusITS0=" << statusITS[0]
608 << "statusITS1=" << statusITS[1]
609 << "statusITS2=" << statusITS[2]
610 << "statusITS3=" << statusITS[3]
611 << "statusITS4=" << statusITS[4]
612 << "statusITS5=" << statusITS[5]
fda75e2d 613 << "trd0=" << hasTrackletTRD[0]
614 << "trd1=" << hasTrackletTRD[1]
615 << "trd2=" << hasTrackletTRD[2]
616 << "trd3=" << hasTrackletTRD[3]
617 << "trd4=" << hasTrackletTRD[4]
618 << "trd5=" << hasTrackletTRD[5]
cedf0381 619 << "TRDdEdxl0=" << trddEdxSum[0]
620 << "TRDdEdxl1=" << trddEdxSum[1]
621 << "TRDdEdxl2=" << trddEdxSum[2]
622 << "TRDdEdxl3=" << trddEdxSum[3]
623 << "TRDdEdxl4=" << trddEdxSum[4]
624 << "TRDdEdxl5=" << trddEdxSum[5]
fda75e2d 625 << "TOFsigmaEl=" << nSigmaTOF
626 << "TPCsigmaEl=" << nSigmaTPC
627 << "TPCdEdx=" << tPCdEdx
628 << "TRDlikeEl=" << likeEleTRD
629 << "TRDlikeEln=" << likeEleTRDn
cedf0381 630 << "trdtruncmean1=" << trdtruncmean1
a8ef1999 631 << "trdtruncmean2=" << trdtruncmean2
fda75e2d 632 << "dcaR=" << b[0]
633 << "dcaZ=" << b[1]
634 << "dca=" << dca
635 << "dcaSR=" << dcaSR
636 << "dcaSZ=" << dcaSZ
637 << "dcaS=" << dcaS
11ff28c5 638 << "hfedcaR=" << hfeb[0]
639 << "hfedcaZ=" << hfeb[1]
640 << "hfedcacovR=" << hfebCov[0]
641 << "hfedcacovZ=" << hfebCov[2]
fda75e2d 642 << "vx=" << vtx[0]
643 << "vy=" << vtx[1]
0e30407a 644 << "vz=" << vtx[2]
cedf0381 645 << "tofdx=" << tofdx
d44e3c84 646 << "tofdz=" << tofdz
647 << "statusTOFtracking=" << hasgoodTOF
648 << "TOFmismatchlevel=" << mismatchlevel
649 << "statusTRDtracking=" << hasTRDin
cedf0381 650 << "ncontrib=" << ncontrib
651 << "mesonID=" << mesonID
652 << "eR=" << eR
653 << "mesonR=" << mesonR
654 << "eZ=" << eZ
655 << "mesonZ=" << mesonZ
656 << "unique=" << unique
657 << "mesonunique=" << mesonunique
658 << "bgcategory=" << bgcategory
659 << "mesonpt=" << mesonPt
11ff28c5 660 << "mesonMomPdg=" << mesonMomPdg
661 << "mesonGMomPdg=" << mesonGMomPdg
662 << "mesonGGMomPdg=" << mesonGGMomPdg
663 << "mesonMomPt=" << mesonMomPt
fda75e2d 664 << "\n";
665 }
666}
667
668
669void AliHFEdebugTreeTask::SetFileName(const char *filename){ fFilename = filename; }
11ff28c5 670//___________________________________________________________
671AliTrackReference *AliHFEdebugTreeTask::FindTrackReference(AliMCParticle *track, Float_t minRadius, Float_t maxRadius, Int_t detectorID)
672{
673 //
674 // Find the track reference
675 //
676 AliTrackReference *ref = NULL, *reftmp;
677 Float_t radius;
678 for(Int_t iref = 0; iref < track->GetNumberOfTrackReferences(); iref++){
679 reftmp = track->GetTrackReference(iref);
680 if(reftmp->DetectorId() != detectorID) continue;
681 radius = reftmp->R();
682 if(radius >= minRadius && radius < maxRadius){
683 ref = reftmp;
684 break;
685 }
686 if(radius > maxRadius) break;
687 }
688 return ref;
689}
690
691//__________________________________________
692Int_t AliHFEdebugTreeTask::GetElecSourceMC(TParticle * const mcpart)
693{
694 // decay particle's origin
695
696 if(!mcpart){
697 AliDebug(1, "no mcparticle, return\n");
698 return -1;
699 }
700
701 if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return AliHFEmcQA::kMisID;
702
703 Int_t origin = -1;
704 Bool_t isFinalOpenCharm = kFALSE;
705
706 Int_t iLabel = mcpart->GetFirstMother();
707 if (iLabel<0){
708 AliDebug(1, "Stack label is negative, return\n");
709 return -1;
710 }
711
712 AliMCParticle *mctrack = NULL;
713 Int_t tmpMomLabel=0;
714 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
715 TParticle *partMother = mctrack->Particle();
716 TParticle *partMotherCopy = mctrack->Particle();
717 Int_t maPdgcode = partMother->GetPdgCode();
718
0e30407a 719 // if the mother is charmed hadron
720 if ( (int(abs(maPdgcode)/100.)%10) == AliHFEmcQA::kCharm || (int(abs(maPdgcode)/1000.)%10) == AliHFEmcQA::kCharm ) {
11ff28c5 721
0e30407a 722 for (Int_t i=0; i<fNparents; i++){
11ff28c5 723 if (abs(maPdgcode)==fParentSelect[0][i]){
724 isFinalOpenCharm = kTRUE;
725 }
0e30407a 726 }
727 if (!isFinalOpenCharm) return -1;
11ff28c5 728
0e30407a 729 // iterate until you find B hadron as a mother or become top ancester
730 for (Int_t i=1; i<fgkMaxIter; i++){
11ff28c5 731
732 Int_t jLabel = partMother->GetFirstMother();
733 if (jLabel == -1){
734 origin = AliHFEmcQA::kDirectCharm;
735 return origin;
736 }
737 if (jLabel < 0){ // safety protection
738 AliDebug(1, "Stack label is negative, return\n");
739 return -1;
740 }
741
742 // if there is an ancester
743 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
744 TParticle* grandMa = mctrack->Particle();
745 Int_t grandMaPDG = grandMa->GetPdgCode();
746
747 for (Int_t j=0; j<fNparents; j++){
0e30407a 748 if (abs(grandMaPDG)==fParentSelect[1][j]){
749 origin = AliHFEmcQA::kBeautyCharm;
750 return origin;
751 }
11ff28c5 752 }
753
754 partMother = grandMa;
0e30407a 755 } // end of iteration
756 } // end of if
757 else if ( (int(abs(maPdgcode)/100.)%10) == AliHFEmcQA::kBeauty || (int(abs(maPdgcode)/1000.)%10) == AliHFEmcQA::kBeauty ) {
758 for (Int_t i=0; i<fNparents; i++){
11ff28c5 759 if (abs(maPdgcode)==fParentSelect[1][i]){
760 origin = AliHFEmcQA::kDirectBeauty;
761 return origin;
762 }
0e30407a 763 }
764 } // end of if
765 else if ( abs(maPdgcode) == 22 ) { //conversion
766
767 tmpMomLabel = partMotherCopy->GetFirstMother();
768 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
769 partMother = mctrack->Particle();
770 maPdgcode = partMother->GetPdgCode();
771 if ( abs(maPdgcode) == 111 ) {
772 origin = AliHFEmcQA::kGammaPi0;
773 return origin;
774 }
775 else if ( abs(maPdgcode) == 221 ) {
776 origin = AliHFEmcQA::kGammaEta;
777 return origin;
778 }
779 else if ( abs(maPdgcode) == 223 ) {
780 origin = AliHFEmcQA::kGammaOmega;
781 return origin;
782 }
783 else if ( abs(maPdgcode) == 333 ) {
784 origin = AliHFEmcQA::kGammaPhi;
785 return origin;
786 }
787 else if ( abs(maPdgcode) == 331 ) {
788 origin = AliHFEmcQA::kGammaEtaPrime;
789 return origin;
790 }
791 else if ( abs(maPdgcode) == 113 ) {
792 origin = AliHFEmcQA::kGammaRho0;
793 return origin;
794 }
795 else origin = AliHFEmcQA::kElse;
796 //origin = kGamma; // finer category above
797 return origin;
798
799 } // end of if
800 else if ( abs(maPdgcode) == 111 ) {
801 origin = AliHFEmcQA::kPi0;
802 return origin;
803 } // end of if
804 else if ( abs(maPdgcode) == 221 ) {
805 origin = AliHFEmcQA::kEta;
806 return origin;
807 } // end of if
808 else if ( abs(maPdgcode) == 223 ) {
809 origin = AliHFEmcQA::kOmega;
810 return origin;
811 } // end of if
812 else if ( abs(maPdgcode) == 333 ) {
813 origin = AliHFEmcQA::kPhi;
814 return origin;
815 } // end of if
816 else if ( abs(maPdgcode) == 331 ) {
817 origin = AliHFEmcQA::kEtaPrime;
818 return origin;
819 } // end of if
820 else if ( abs(maPdgcode) == 113 ) {
821 origin = AliHFEmcQA::kRho0;
822 return origin;
823 } // end of if
824 else{
11ff28c5 825 origin = AliHFEmcQA::kElse;
0e30407a 826 }
827 return origin;
11ff28c5 828}
cedf0381 829
830//______________________________________________________
831Bool_t AliHFEdebugTreeTask::CheckITSstatus( const AliESDtrack * const esdtrack, Int_t layer) const {
832 //
833 // Check whether ITS area is dead
834 //
835 Int_t itsStatus = 0;
836 Int_t det;
837 Float_t xloc, zloc;
838 esdtrack->GetITSModuleIndexInfo(layer, det, itsStatus, xloc, zloc);
839 Bool_t status;
840 switch(itsStatus){
841 case 2: status = kFALSE; break;
842 case 3: status = kFALSE; break;
843 case 7: status = kFALSE; break;
844 default: status = kTRUE;
845 }
846 return status;
847}
848