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 **************************************************************************/
17 // the tree is represented as reduced events
20 // M.Fasel <M.Fasel@gsi.de>
26 #include <TParticle.h>
30 #include "AliAnalysisManager.h"
31 #include "AliAODMCHeader.h"
32 #include "AliAnalysisUtils.h"
33 #include "AliCentrality.h"
34 #include "AliESDEvent.h"
35 #include "AliESDtrack.h"
36 #include "AliESDVertex.h"
37 #include "AliInputEventHandler.h"
38 #include "AliHFEcuts.h"
39 #include "AliHFEextraCuts.h"
40 #include "AliHFEpidTPC.h"
41 #include "AliHFEreducedEvent.h"
42 #include "AliHFEreducedTrack.h"
43 #include "AliHFEreducedMCParticle.h"
44 #include "AliHFEsignalCuts.h"
46 #include "AliMCEvent.h"
47 #include "AliMCParticle.h"
48 #include "AliMCEventHandler.h"
49 #include "AliMultiplicity.h"
50 #include "AliPIDResponse.h"
51 #include "AliVEvent.h"
52 #include "AliVTrack.h"
53 #include "AliVCluster.h"
54 #include "AliVVZERO.h"
56 #include "AliTRDTriggerAnalysis.h"
57 #include "TTreeStream.h"
59 #include "AliHFEreducedEventCreatorESD.h"
61 ClassImp(AliHFEreducedEventCreatorESD)
63 AliHFEreducedEventCreatorESD::AliHFEreducedEventCreatorESD():
72 fTRDTriggerAnalysis(NULL),
77 fRemoveFirstEvent(kFALSE)
80 // Default constructor
84 AliHFEreducedEventCreatorESD::AliHFEreducedEventCreatorESD(const char *name):
85 AliAnalysisTaskSE(name),
93 fTRDTriggerAnalysis(NULL),
98 fRemoveFirstEvent(kFALSE)
101 // Default constructor
103 fTPCpid = new AliHFEpidTPC("QAtpcPID");
104 fAnalysisUtils = new AliAnalysisUtils;
105 fTRDTriggerAnalysis = new AliTRDTriggerAnalysis();
106 DefineOutput(1, TTree::Class());
109 AliHFEreducedEventCreatorESD::~AliHFEreducedEventCreatorESD(){
111 // Default destructor
113 if(fAnalysisUtils) delete fAnalysisUtils;
114 if(fTPCpid) delete fTPCpid;
115 if(fTRDTriggerAnalysis) delete fTRDTriggerAnalysis;
116 if(fHFEevent) delete fHFEevent;
117 if(fSignalCuts) delete fSignalCuts;
118 if(fTrackCuts) delete fTrackCuts;
121 void AliHFEreducedEventCreatorESD::UserCreateOutputObjects(){
123 // Create debug tree, signal cuts and track cuts
128 // printf("testa\n");
129 fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
132 fTrackCuts = new AliHFEcuts("fTrackCuts", "Basic HFE track cuts");
133 fTrackCuts->CreateStandardCuts();
135 fTrackCuts->SetMinNClustersTPC(fNclustersTPC);
136 fTrackCuts->SetMinRatioTPCclusters(0);
137 fTrackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
138 fTrackCuts->SetMinNClustersTPCPID(fNclustersTPCPID);
139 fTrackCuts->SetMinNClustersITS(fNclustersITS);
141 fTrackCuts->SetUseMixedVertex(kTRUE);
142 fTrackCuts->SetVertexRange(10.);
144 fTrackCuts->Initialize();
147 fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
149 fHFEevent = new AliHFEreducedEvent;
151 fHFEtree = new TTree("HFEtree", "HFE event tree");
152 fHFEtree->Branch("HFEevent", "AliHFEreducedEvent", fHFEevent,128000,0);
153 PostData(1, fHFEtree);
156 void AliHFEreducedEventCreatorESD::UserExec(Option_t *){
158 // User Exec: Fill debug Tree
162 AliPIDResponse *pid = NULL;
163 AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
165 // printf("testb\n");
166 pid = handler->GetPIDResponse();
168 AliError("No Handler");
171 // printf("testc\n");
172 AliError("No PID response");
176 AliError("No Input event");
180 if(fRemoveFirstEvent){
181 if(fAnalysisUtils->IsFirstEventInChunk(fInputEvent)) return;
184 AliDebug(1, "Event Selected");
186 AliESDtrack copyTrack;
189 Bool_t mcthere = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) != NULL;
191 fTrackCuts->SetMCEvent(fMCEvent);
192 fSignalCuts->SetMCEvent(fMCEvent);
195 fTrackCuts->SetRecEvent(fInputEvent);
197 if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
198 AliDebug(1, "Event rejected by the event cuts\n");
202 // reject pile up in case of pp
203 AliESDEvent *event = dynamic_cast<AliESDEvent *>(fInputEvent);
205 TString beamtype = event->GetBeamType();
206 //printf("beamtype %s\n",(const char*)beamtype);
207 if (strstr(beamtype,"p-p")) {
208 //printf("Reject\n");
209 if(fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5)){
210 AliDebug(1, "Event flagged as pileup\n");
217 fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
219 fExtraCuts->SetRecEventInfo(fInputEvent);
222 // Make Reduced Event
223 //AliHFEreducedEvent hfeevent;
224 fHFEevent->~AliHFEreducedEvent();
225 new(fHFEevent)AliHFEreducedEvent();
228 fHFEevent->SetRunNumber(fInputEvent->GetRunNumber());
231 UInt_t trigger = fInputHandler->IsEventSelected();
232 if(trigger & AliVEvent::kMB) fHFEevent->SetMBTrigger();
233 if(trigger & AliVEvent::kCentral) fHFEevent->SetCentralTrigger();
234 if(trigger & AliVEvent::kSemiCentral) fHFEevent->SetCentralTrigger();
235 if(trigger & AliVEvent::kEMCEJE) fHFEevent->SetEMCALTrigger();
236 fTRDTriggerAnalysis->CalcTriggers(event);
237 if(fTRDTriggerAnalysis->IsFired(AliTRDTriggerAnalysis::kHSE)) fHFEevent->SetTRDSETrigger();
238 if(fTRDTriggerAnalysis->IsFired(AliTRDTriggerAnalysis::kHQU)) fHFEevent->SetTRDDQTrigger();
240 // Get Primary Vertex
241 const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
244 fHFEevent->SetVX(vtx[0]);
245 fHFEevent->SetVY(vtx[1]);
246 fHFEevent->SetVZ(vtx[2]);
247 Int_t ncontrib(vertex->GetNContributors());
248 fHFEevent->SetNContribVertex(ncontrib);
250 vertex->GetCovarianceMatrix(vcov);
251 fHFEevent->SetVertexResolution(TMath::Sqrt(vcov[5]));
252 // Get Primary Vertex from SPD
253 const AliVVertex *vertexSPD = event->GetPrimaryVertexSPD();
255 memset(vtx, 0, sizeof(Double_t) *3);
256 vertexSPD->GetXYZ(vtx);
257 fHFEevent->SetVXSPD(vtx[0]);
258 fHFEevent->SetVYSPD(vtx[1]);
259 fHFEevent->SetVZSPD(vtx[2]);
260 fHFEevent->SetNContribVertexSPD(vertexSPD->GetNContributors());
261 memset(vcov, 0, sizeof(Double_t)*6);
262 vertex->GetCovarianceMatrix(vcov);
263 fHFEevent->SetVertexResolutionSPD(TMath::Sqrt(vcov[5]));
267 AliCentrality *hicent = fInputEvent->GetCentrality();
268 fHFEevent->SetCentrality(
269 hicent->GetCentralityPercentile("V0M"),
270 hicent->GetCentralityPercentile("V0A"),
271 hicent->GetCentralityPercentile("V0C"),
272 hicent->GetCentralityPercentile("TKL"),
273 hicent->GetCentralityPercentile("TRK"),
274 hicent->GetCentralityPercentile("ZNA")
277 // Get VZERO Information
278 AliVVZERO *vzeroinfo = fInputEvent->GetVZEROData();
279 if(vzeroinfo) fHFEevent->SetV0Multiplicity(vzeroinfo->GetMTotV0A(), vzeroinfo->GetMTotV0C());
281 // Get ZDC Information
282 AliVZDC *zdcinfo = fInputEvent->GetZDCData();
283 if(zdcinfo) fHFEevent->SetZDCEnergy(zdcinfo->GetZNAEnergy(), zdcinfo->GetZNCEnergy(), zdcinfo->GetZPAEnergy(), zdcinfo->GetZPCEnergy());
285 // Set SPD multiplicity
286 const AliMultiplicity *mult = event->GetMultiplicity();
287 if(mult) fHFEevent->SetSPDMultiplicity(mult->GetNumberOfTracklets());
290 // Loop on MC tracks only
292 AliMCParticle *mctrack(NULL);
296 for(Int_t itrack = 0; itrack < fMCEvent->GetNumberOfTracks(); itrack++) {
297 mctrack = (AliMCParticle *)(fMCEvent->GetTrack(itrack));
298 if(!mctrack) continue;
299 AliHFEreducedMCParticle hfemcpart;
300 if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) hfemcpart.SetSignal();
302 hfemcpart.SetSignedPt(mctrack->Pt(), mctrack->Charge() > 0.);
303 hfemcpart.SetP(mctrack->P());
304 hfemcpart.SetEta(mctrack->Eta());
305 hfemcpart.SetPhi(mctrack->Phi());
306 hfemcpart.SetPdg(mctrack->PdgCode());
308 // Get Production Vertex in radial direction
309 hfemcpart.SetProductionVertex(mctrack->Xv(),mctrack->Yv(),mctrack->Zv());
311 // Get Mother PDG code of the particle
312 Int_t motherlabel = TMath::Abs(mctrack->GetMother());
313 if(motherlabel >= 0 && motherlabel < fMCEvent->GetNumberOfTracks()){
314 AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherlabel));
315 if(mother) hfemcpart.SetMotherPdg(mother->PdgCode());
320 if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
321 else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
322 else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
323 else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
324 else if(TMath::Abs(mctrack->PdgCode()) == 11) source = 4;
326 hfemcpart.SetSource(source);
328 fHFEevent->AddMCParticle(&hfemcpart);
333 // Loop on reconstructed tracks
335 TArrayI arraytrack(fInputEvent->GetNumberOfTracks());
338 AliESDtrack *track = 0x0;
339 for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
341 track = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(itrack));
343 // Cut track (Only basic track cuts)
344 // printf("testv\n");
345 if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
350 AliHFEreducedTrack hfetrack;
351 hfetrack.SetSignedPt(track->Pt(), track->Charge() > 0);
352 hfetrack.SetP(track->P());
353 hfetrack.SetEta(track->Eta());
354 hfetrack.SetPhi(track->Phi());
355 hfetrack.SetTPCmomentum(track->GetTPCmomentum());
358 hfetrack.SetTrackID(track->GetID());
361 ULong_t status = track->GetStatus();
362 if((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) hfetrack.SetITSrefit();
363 if((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) hfetrack.SetTPCrefit();
364 if((status & AliVTrack::kTOFpid) == AliVTrack::kTOFpid) hfetrack.SetTOFpid();
365 //if((status & AliVTrack::kTOFmismatch) == AliVTrack::kTOFmismatch) hfetrack.SetTOFmismatch();
366 if(IsTOFmismatch(track, pid)) hfetrack.SetTOFmismatch(); // New version suggested by Pietro Antonioli
367 if(track->IsEMCAL()) hfetrack.SetEMCALpid();
368 // no filter bits available for ESDs
371 // Fill Monte-Carlo Information
372 Int_t label = TMath::Abs(track->GetLabel());
373 if(label && label < fMCEvent->GetNumberOfTracks())
374 mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
375 if(!mctrack) continue;
376 if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) hfetrack.SetMCSignal();
378 hfetrack.SetMCSignedPt(mctrack->Pt(),mctrack->Charge() > 0.);
379 hfetrack.SetMCP(mctrack->P());
380 hfetrack.SetMCEta(mctrack->Eta());
381 hfetrack.SetMCPhi(mctrack->Phi());
382 hfetrack.SetMCPDG(mctrack->PdgCode());
384 // Get Production Vertex in radial direction
385 hfetrack.SetMCProdVtx(mctrack->Xv(),mctrack->Yv(),mctrack->Zv());
387 // Get Mother PDG code of the particle
388 Int_t motherlabel = TMath::Abs(mctrack->GetMother());
389 if(motherlabel >= 0 && motherlabel < fMCEvent->GetNumberOfTracks()){
390 AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherlabel));
391 if(mother) hfetrack.SetMCMotherPdg(mother->PdgCode());
396 if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
397 else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
398 else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
399 else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
400 else if(TMath::Abs(mctrack->PdgCode()) == 11) source = 4;
402 hfetrack.SetMCSource(source);
406 Float_t dcaxy = -999.,
408 fExtraCuts->GetImpactParameters((AliVTrack *)track,dcaxy,dcaz);
409 hfetrack.SetDCA(dcaxy, dcaz);
411 // Different number of clusters definitions
412 Int_t nclustersITS(track->GetITSclusters(NULL)),
413 nclustersTPC(track->GetTPCNcls()),
414 nclustersTPCall(track->GetTPCClusterMap().CountBits()),
415 nclustersTPCshared(0);
416 UChar_t nfindableTPC = track->GetTPCNclsF();
417 const TBits &sharedTPC = track->GetTPCSharedMap();
418 for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
419 hfetrack.SetITSnclusters(nclustersITS);
420 hfetrack.SetTPCnclusters(nclustersTPC);
421 hfetrack.SetTRDnclusters(track->GetTRDncls());
422 hfetrack.SetTPCnclustersPID(track->GetTPCsignalN());
423 hfetrack.SetTPCcrossedRows(track->GetTPCCrossedRows());
424 hfetrack.SetTPCnclustersAll(nclustersTPCall);
425 hfetrack.SetTPCsharedClusters(nclustersTPCshared);
426 hfetrack.SetTPCclusterRatio(nfindableTPC ? static_cast<Float_t>(nclustersTPC)/static_cast<Float_t>(nfindableTPC) : 0);
427 hfetrack.SetTPCclusterRatioAll(nfindableTPC ? static_cast<Float_t>(nclustersTPCall)/static_cast<Float_t>(nfindableTPC) : 0);
428 UChar_t itsPixel = track->GetITSClusterMap();
429 for(int ily = 0; ily < 6; ily++)
430 if(TESTBIT(itsPixel, ily)) hfetrack.SetITScluster(ily);
432 // TRD related quantities (Yvonne)
433 Int_t nslices = track->GetNumberOfTRDslices();
434 hfetrack.SetTRDntrackletsPID(track->GetTRDntrackletsPID());
435 hfetrack.SetTRDnslices(nslices);
436 hfetrack.SetTRDchi2(track->GetTRDchi2());
438 for(Int_t iplane = 0; iplane < 6; iplane++){
440 for(Int_t isl = 0; isl < nslices; isl++){
441 if(track->GetTRDntrackletsPID()>0){
442 if(track->GetTRDslice(iplane, isl)>0.001) nslicetemp++;
445 if(nslicetemp > 0) hfetrack.SetTRDstatus(iplane);
451 if(fExtraCuts->IsKinkDaughter(track)) kink = 1;
454 Int_t kinkmotherpass = 0;
455 if(track->GetKinkIndex(0) != 0 && !kink) kinkmotherpass = 1;
458 Int_t id(track->GetID());
459 for(Int_t l=0; l < counterdc; l++){
460 Int_t iTrack2 = arraytrack.At(l);
462 hfetrack.SetDoubleCounted();
466 // Add the id at this place
467 arraytrack.AddAt(id,counterdc);
471 hfetrack.SetTPCdEdx(track->GetTPCsignal());
472 hfetrack.SetTPCsigmaEl(pid->NumberOfSigmasTPC(track, AliPID::kElectron));
473 hfetrack.SetTOFsigmaEl(pid->NumberOfSigmasTOF(track, AliPID::kElectron));
474 hfetrack.SetTOFmismatchProbability(pid->GetTOFMismatchProbability(track));
476 copyTrack.~AliESDtrack();
477 new(©Track) AliESDtrack(*track);
478 if(fTPCpid->HasCentralityCorrection()) fTPCpid->ApplyCentralityCorrection(©Track, static_cast<Double_t>(ncontrib),AliHFEpidObject::kESDanalysis);
479 if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(©Track, AliHFEpidObject::kESDanalysis);
480 hfetrack.SetTPCsigmaElCorrected(pid->NumberOfSigmasTPC(©Track, AliPID::kElectron));
481 hfetrack.SetTPCdEdxCorrected(copyTrack.GetTPCsignal());
482 if(track->IsEMCAL()){
484 Double_t emcalEnergyOverP = -1.,
485 showershape[4] = {0.,0.,0.,0.};
486 hfetrack.SetEMCALSigmaEl(pid->NumberOfSigmasEMCAL(track, AliPID::kElectron, emcalEnergyOverP, &showershape[0]));
487 hfetrack.SetEMCALEoverP(emcalEnergyOverP);
488 hfetrack.SetEMCALShowerShape(showershape);
491 // Track finished, add NOW to the Event
492 fHFEevent->AddTrack(&hfetrack);
496 // Fill the debug tree
497 //AliInfo(Form("Number of tracks: %d\n", fHFEevent->GetNumberOfTracks()));
498 //AliInfo(Form("Number of MC particles: %d\n", fHFEevent->GetNumberOfMCParticles()));
502 PostData(1, fHFEtree);
505 void AliHFEreducedEventCreatorESD::Terminate(Option_t *){
509 AliInfo("terminating...\n");
513 Bool_t AliHFEreducedEventCreatorESD::IsTOFmismatch(const AliVTrack *const track, const AliPIDResponse *const pid) const {
517 Double_t probs[AliPID::kSPECIESC];
518 AliPIDResponse::EDetPidStatus status = pid->ComputeTOFProbability(track, AliPID::kSPECIESC, probs);
519 return status == AliPIDResponse::kDetMismatch;