]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/STRANGENESS/Cascades/AliAnalysisTaskExtractCascade.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
facilitate LEGO train running
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Cascades / AliAnalysisTaskExtractCascade.cxx
CommitLineData
76029adc 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// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
17//
18// Modified version of AliAnalysisTaskCheckCascade.cxx.
19// This is a 'hybrid' output version, in that it uses a classic TTree
20// ROOT object to store the candidates, plus a couple of histograms filled on
21// a per-event basis for storing variables too numerous to put in a tree.
22//
76029adc 23//
24// --- Algorithm Description
25// 1. Loop over primaries in stack to acquire generated charged Xi
9a8f3aee 26// 2. Loop over stack to find Cascades, fill TH3Fs "PrimRawPt"s for Efficiency
76029adc 27// 3. Perform Physics Selection
28// 4. Perform Primary Vertex |z|<10cm selection
9a8f3aee 29// 5. Perform Primary Vertex NoTPCOnly vertexing selection
76029adc 30// 6. Perform Pileup Rejection
31// 7. Analysis Loops:
32// 7a. Fill TH3Fs "PrimAnalysisPt" for control purposes only
76029adc 33//
34// Please Report Any Bugs!
35//
36// --- David Dobrigkeit Chinellato
37// (david.chinellato@gmail.com)
38//
39// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
40
41class TTree;
42class TParticle;
43class TVector3;
44
45//class AliMCEventHandler;
46//class AliMCEvent;
47//class AliStack;
48
49class AliESDVertex;
50class AliAODVertex;
51class AliESDv0;
52class AliAODv0;
53
54#include <Riostream.h>
55#include "TList.h"
56#include "TH1.h"
57#include "TH2.h"
58#include "TH3.h"
59#include "TFile.h"
60#include "THnSparse.h"
61#include "TVector3.h"
62#include "TCanvas.h"
63#include "TMath.h"
64#include "TLegend.h"
65#include "AliLog.h"
66
67#include "AliESDEvent.h"
68#include "AliAODEvent.h"
69#include "AliV0vertexer.h"
70#include "AliCascadeVertexer.h"
71#include "AliESDpid.h"
72#include "AliESDtrack.h"
73#include "AliESDtrackCuts.h"
74#include "AliInputEventHandler.h"
75#include "AliAnalysisManager.h"
76#include "AliMCEventHandler.h"
77#include "AliMCEvent.h"
78#include "AliStack.h"
79
80#include "AliCFContainer.h"
81#include "AliMultiplicity.h"
82#include "AliAODMCParticle.h"
83#include "AliESDcascade.h"
84#include "AliAODcascade.h"
85#include "AliESDUtils.h"
86#include "AliGenEventHeader.h"
87
88#include "AliAnalysisTaskExtractCascade.h"
89
90using std::cout;
91using std::endl;
92
93ClassImp(AliAnalysisTaskExtractCascade)
94
95AliAnalysisTaskExtractCascade::AliAnalysisTaskExtractCascade()
96 : AliAnalysisTaskSE(), fListHist(0), fTreeCascade(0), fPIDResponse(0), fESDtrackCuts(0),
97 fkIsNuclear ( kFALSE ),
98 fkLowEnergyPP ( kFALSE ),
99
100//------------------------------------------------
101// HISTOGRAMS
102// --- Filled on an Event-by-event basis
103//------------------------------------------------
104 fHistV0MultiplicityBeforeTrigSel(0),
105 fHistV0MultiplicityForTrigEvt(0),
106 fHistV0MultiplicityForSelEvt(0),
107 fHistV0MultiplicityForSelEvtNoTPCOnly(0),
108 fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
109 fHistMultiplicityBeforeTrigSel(0),
110 fHistMultiplicityForTrigEvt(0),
111 fHistMultiplicity(0),
112 fHistMultiplicityNoTPCOnly(0),
113 fHistMultiplicityNoTPCOnlyNoPileup(0),
114 fHistPVx(0),
115 fHistPVy(0),
116 fHistPVz(0),
117 fHistPVxAnalysis(0),
118 fHistPVyAnalysis(0),
119 fHistPVzAnalysis(0)
120{
121 // Dummy Constructor
122}
123
124AliAnalysisTaskExtractCascade::AliAnalysisTaskExtractCascade(const char *name)
125 : AliAnalysisTaskSE(name), fListHist(0), fTreeCascade(0), fPIDResponse(0), fESDtrackCuts(0),
126 fkIsNuclear ( kFALSE ),
127 fkLowEnergyPP ( kFALSE ),
128
129//------------------------------------------------
130// HISTOGRAMS
131// --- Filled on an Event-by-event basis
132//------------------------------------------------
133 fHistV0MultiplicityBeforeTrigSel(0),
134 fHistV0MultiplicityForTrigEvt(0),
135 fHistV0MultiplicityForSelEvt(0),
136 fHistV0MultiplicityForSelEvtNoTPCOnly(0),
137 fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
138 fHistMultiplicityBeforeTrigSel(0),
139 fHistMultiplicityForTrigEvt(0),
140 fHistMultiplicity(0),
141 fHistMultiplicityNoTPCOnly(0),
142 fHistMultiplicityNoTPCOnlyNoPileup(0),
143 fHistPVx(0),
144 fHistPVy(0),
145 fHistPVz(0),
146 fHistPVxAnalysis(0),
147 fHistPVyAnalysis(0),
148 fHistPVzAnalysis(0)
149{
150 // Constructor
151 // Output slot #0 writes into a TList container (Cascade)
152 DefineOutput(1, TList::Class());
153 DefineOutput(2, TTree::Class());
154}
155
156
157AliAnalysisTaskExtractCascade::~AliAnalysisTaskExtractCascade()
158{
159//------------------------------------------------
160// DESTRUCTOR
161//------------------------------------------------
162
163 if (fListHist){
164 delete fListHist;
165 fListHist = 0x0;
166 }
167 if (fTreeCascade){
168 delete fTreeCascade;
169 fTreeCascade = 0x0;
170 }
171 //cleanup esd track cuts object too...
172 if (fESDtrackCuts){
173 delete fESDtrackCuts;
174 fESDtrackCuts = 0x0;
175 }
176
177}
178
179//________________________________________________________________________
180void AliAnalysisTaskExtractCascade::UserCreateOutputObjects()
181{
182 OpenFile(2);
183 // Called once
184
185//------------------------------------------------
186
187 fTreeCascade = new TTree("fTreeCascade","CascadeCandidates");
188
189//------------------------------------------------
190// fTreeCascade Branch definitions - Cascade Tree
191//------------------------------------------------
192
193//------------------------------------------------
194// fTreeCascade Branch definitions
195//------------------------------------------------
196
197//-----------BASIC-INFO---------------------------
198/* 1*/ fTreeCascade->Branch("fTreeCascVarCharge",&fTreeCascVarCharge,"fTreeCascVarCharge/I");
199/* 2*/ fTreeCascade->Branch("fTreeCascVarMassAsXi",&fTreeCascVarMassAsXi,"fTreeCascVarMassAsXi/F");
200/* 3*/ fTreeCascade->Branch("fTreeCascVarMassAsOmega",&fTreeCascVarMassAsOmega,"fTreeCascVarMassAsOmega/F");
201/* 4*/ fTreeCascade->Branch("fTreeCascVarPt",&fTreeCascVarPt,"fTreeCascVarPt/F");
202/* 5*/ fTreeCascade->Branch("fTreeCascVarRapXi",&fTreeCascVarRapXi,"fTreeCascVarRapXi/F");
203/* 6*/ fTreeCascade->Branch("fTreeCascVarRapOmega",&fTreeCascVarRapOmega,"fTreeCascVarRapOmega/F");
204/* 7*/ fTreeCascade->Branch("fTreeCascVarNegEta",&fTreeCascVarNegEta,"fTreeCascVarNegEta/F");
205/* 8*/ fTreeCascade->Branch("fTreeCascVarPosEta",&fTreeCascVarPosEta,"fTreeCascVarPosEta/F");
206/* 9*/ fTreeCascade->Branch("fTreeCascVarBachEta",&fTreeCascVarBachEta,"fTreeCascVarBachEta/F");
207//-----------INFO-FOR-CUTS------------------------
208/*10*/ fTreeCascade->Branch("fTreeCascVarDCACascDaughters",&fTreeCascVarDCACascDaughters,"fTreeCascVarDCACascDaughters/F");
209/*11*/ fTreeCascade->Branch("fTreeCascVarDCABachToPrimVtx",&fTreeCascVarDCABachToPrimVtx,"fTreeCascVarDCABachToPrimVtx/F");
210/*12*/ fTreeCascade->Branch("fTreeCascVarDCAV0Daughters",&fTreeCascVarDCAV0Daughters,"fTreeCascVarDCAV0Daughters/F");
211/*13*/ fTreeCascade->Branch("fTreeCascVarDCAV0ToPrimVtx",&fTreeCascVarDCAV0ToPrimVtx,"fTreeCascVarDCAV0ToPrimVtx/F");
212/*14*/ fTreeCascade->Branch("fTreeCascVarDCAPosToPrimVtx",&fTreeCascVarDCAPosToPrimVtx,"fTreeCascVarDCAPosToPrimVtx/F");
213/*15*/ fTreeCascade->Branch("fTreeCascVarDCANegToPrimVtx",&fTreeCascVarDCANegToPrimVtx,"fTreeCascVarDCANegToPrimVtx/F");
214/*16*/ fTreeCascade->Branch("fTreeCascVarCascCosPointingAngle",&fTreeCascVarCascCosPointingAngle,"fTreeCascVarCascCosPointingAngle/F");
215/*17*/ fTreeCascade->Branch("fTreeCascVarCascRadius",&fTreeCascVarCascRadius,"fTreeCascVarCascRadius/F");
216/*18*/ fTreeCascade->Branch("fTreeCascVarV0Mass",&fTreeCascVarV0Mass,"fTreeCascVarV0Mass/F");
217/*19*/ fTreeCascade->Branch("fTreeCascVarV0CosPointingAngle",&fTreeCascVarV0CosPointingAngle,"fTreeCascVarV0CosPointingAngle/F");
218/*20*/ fTreeCascade->Branch("fTreeCascVarV0Radius",&fTreeCascVarV0Radius,"fTreeCascVarV0Radius/F");
219/*21*/ fTreeCascade->Branch("fTreeCascVarLeastNbrClusters",&fTreeCascVarLeastNbrClusters,"fTreeCascVarLeastNbrClusters/I");
220//-----------MULTIPLICITY-INFO--------------------
221/*22*/ fTreeCascade->Branch("fTreeCascVarMultiplicity",&fTreeCascVarMultiplicity,"fTreeCascVarMultiplicity/I");
222//-----------DECAY-LENGTH-INFO--------------------
223/*23*/ fTreeCascade->Branch("fTreeCascVarDistOverTotMom",&fTreeCascVarDistOverTotMom,"fTreeCascVarDistOverTotMom/F");
224//------------------------------------------------
225/*24*/ fTreeCascade->Branch("fTreeCascVarNegNSigmaPion",&fTreeCascVarNegNSigmaPion,"fTreeCascVarNegNSigmaPion/F");
226/*25*/ fTreeCascade->Branch("fTreeCascVarNegNSigmaProton",&fTreeCascVarNegNSigmaProton,"fTreeCascVarNegNSigmaProton/F");
227/*26*/ fTreeCascade->Branch("fTreeCascVarPosNSigmaPion",&fTreeCascVarPosNSigmaPion,"fTreeCascVarPosNSigmaPion/F");
228/*27*/ fTreeCascade->Branch("fTreeCascVarPosNSigmaProton",&fTreeCascVarPosNSigmaProton,"fTreeCascVarPosNSigmaProton/F");
229/*28*/ fTreeCascade->Branch("fTreeCascVarBachNSigmaPion",&fTreeCascVarBachNSigmaPion,"fTreeCascVarBachNSigmaPion/F");
230/*29*/ fTreeCascade->Branch("fTreeCascVarBachNSigmaKaon",&fTreeCascVarBachNSigmaKaon,"fTreeCascVarBachNSigmaKaon/F");
231
232//------------------------------------------------
233// Particle Identification Setup
234//------------------------------------------------
235
236 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
237 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
238 fPIDResponse = inputHandler->GetPIDResponse();
239
240// Multiplicity
241
242 if(! fESDtrackCuts ){
243 fESDtrackCuts = new AliESDtrackCuts();
244 }
245
246//------------------------------------------------
247// V0 Multiplicity Histograms
248//------------------------------------------------
249
250 // Create histograms
251 OpenFile(1);
252 fListHist = new TList();
253 fListHist->SetOwner(); // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
254
255
256 if(! fHistV0MultiplicityBeforeTrigSel) {
257 fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel",
258 "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events",
259 25, 0, 25);
260 fListHist->Add(fHistV0MultiplicityBeforeTrigSel);
261 }
262
263 if(! fHistV0MultiplicityForTrigEvt) {
264 fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt",
265 "V0s per event (for triggered evt);Nbr of V0s/Evt;Events",
266 25, 0, 25);
267 fListHist->Add(fHistV0MultiplicityForTrigEvt);
268 }
269
270 if(! fHistV0MultiplicityForSelEvt) {
271 fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt",
272 "V0s per event;Nbr of V0s/Evt;Events",
273 25, 0, 25);
274 fListHist->Add(fHistV0MultiplicityForSelEvt);
275 }
276
277 if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {
278 fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly",
279 "V0s per event;Nbr of V0s/Evt;Events",
280 25, 0, 25);
281 fListHist->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);
282 }
283 if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {
284 fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup",
285 "V0s per event;Nbr of V0s/Evt;Events",
286 25, 0, 25);
287 fListHist->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);
288 }
289
290//------------------------------------------------
291// Track Multiplicity Histograms
292//------------------------------------------------
293
294 if(! fHistMultiplicityBeforeTrigSel) {
295 fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel",
296 "Tracks per event;Nbr of Tracks;Events",
297 200, 0, 200);
298 fListHist->Add(fHistMultiplicityBeforeTrigSel);
299 }
300 if(! fHistMultiplicityForTrigEvt) {
301 fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt",
302 "Tracks per event;Nbr of Tracks;Events",
303 200, 0, 200);
304 fListHist->Add(fHistMultiplicityForTrigEvt);
305 }
306 if(! fHistMultiplicity) {
307 fHistMultiplicity = new TH1F("fHistMultiplicity",
308 "Tracks per event;Nbr of Tracks;Events",
309 200, 0, 200);
310 fListHist->Add(fHistMultiplicity);
311 }
312 if(! fHistMultiplicityNoTPCOnly) {
313 fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly",
314 "Tracks per event;Nbr of Tracks;Events",
315 200, 0, 200);
316 fListHist->Add(fHistMultiplicityNoTPCOnly);
317 }
318 if(! fHistMultiplicityNoTPCOnlyNoPileup) {
319 fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup",
320 "Tracks per event;Nbr of Tracks;Events",
321 200, 0, 200);
322 fListHist->Add(fHistMultiplicityNoTPCOnlyNoPileup);
323 }
324
325//----------------------------------
326// Primary Vertex Position Histos
327//----------------------------------
328
329 if(! fHistPVx) {
330 fHistPVx = new TH1F("fHistPVx",
331 "PV x position;Nbr of Evts;x",
332 2000, -0.5, 0.5);
333 fListHist->Add(fHistPVx);
334 }
335 if(! fHistPVy) {
336 fHistPVy = new TH1F("fHistPVy",
337 "PV y position;Nbr of Evts;y",
338 2000, -0.5, 0.5);
339 fListHist->Add(fHistPVy);
340 }
341 if(! fHistPVz) {
342 fHistPVz = new TH1F("fHistPVz",
343 "PV z position;Nbr of Evts;z",
344 400, -20, 20);
345 fListHist->Add(fHistPVz);
346 }
347
348 if(! fHistPVxAnalysis) {
349 fHistPVxAnalysis = new TH1F("fHistPVxAnalysis",
350 "PV x position;Nbr of Evts;x",
351 2000, -0.5, 0.5);
352 fListHist->Add(fHistPVxAnalysis);
353 }
354 if(! fHistPVyAnalysis) {
355 fHistPVyAnalysis = new TH1F("fHistPVyAnalysis",
356 "PV y position;Nbr of Evts;y",
357 2000, -0.5, 0.5);
358 fListHist->Add(fHistPVyAnalysis);
359 }
360 if(! fHistPVzAnalysis) {
361 fHistPVzAnalysis = new TH1F("fHistPVzAnalysis",
362 "PV z position;Nbr of Evts;z",
363 400, -20, 20);
364 fListHist->Add(fHistPVzAnalysis);
365 }
366
367 //List of Histograms: Normal
368 PostData(1, fListHist);
369
370 //TTree Object: Saved to base directory. Should cache to disk while saving.
371 //(Important to avoid excessive memory usage, particularly when merging)
372 PostData(2, fTreeCascade);
373
374}// end UserCreateOutputObjects
375
376
377//________________________________________________________________________
378void AliAnalysisTaskExtractCascade::UserExec(Option_t *)
379{
380 // Main loop
381 // Called for each event
382
383 AliESDEvent *lESDevent = 0x0;
384
385 Int_t lNumberOfV0s = -1;
386 Double_t lTrkgPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
387 Double_t lBestPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
388 Double_t lMagneticField = -10.;
389
390 // Connect to the InputEvent
391 // After these lines, we should have an ESD/AOD event + the number of V0s in it.
392
393 // Appropriate for ESD analysis!
394
395 lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
396 if (!lESDevent) {
397 AliWarning("ERROR: lESDevent not available \n");
398 return;
399 }
400
401/* --- Acquisition of exact event ID
402 fTreeVariableRunNumber = lESDevent->GetRunNumber();
403 fTreeVariableEventNumber =
404 ( ( ((ULong64_t)lESDevent->GetPeriodNumber() ) << 36 ) |
405 ( ((ULong64_t)lESDevent->GetOrbitNumber () ) << 12 ) |
406 ((ULong64_t)lESDevent->GetBunchCrossNumber() ) );
407*/
408
409//------------------------------------------------
410// Multiplicity Information Acquistion
411//------------------------------------------------
412
413 //REVISED multiplicity estimator after 'multiplicity day' (2011)
414 Int_t lMultiplicity = -100;
415
416 //testing purposes
417 if(fkIsNuclear == kFALSE) lMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(lESDevent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
418
419 //---> If this is a nuclear collision, then go nuclear on "multiplicity" variable...
420 //---> Warning: Experimental
421 if(fkIsNuclear == kTRUE){
422 AliCentrality* centrality;
423 centrality = lESDevent->GetCentrality();
424 lMultiplicity = ( ( Int_t ) ( centrality->GetCentralityPercentile( "V0M" ) ) );
425 if (centrality->GetQuality()>1) {
426 PostData(1, fListHist);
427 PostData(2, fTreeCascade);
428 return;
429 }
430 }
431
432 //Set variable for filling tree afterwards!
433 //---> pp case......: GetReferenceMultiplicity
434 //---> Pb-Pb case...: Centrality by V0M
435
436 fTreeCascVarMultiplicity = lMultiplicity;
437
438 fHistV0MultiplicityBeforeTrigSel->Fill ( lESDevent->GetNumberOfV0s() );
439 fHistMultiplicityBeforeTrigSel->Fill ( lMultiplicity );
440
441//------------------------------------------------
442// Physics Selection
443//------------------------------------------------
444
445 UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
446 Bool_t isSelected = 0;
447 isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
448
449 //pp at 2.76TeV: special case, ignore FastOnly
450 if ( (fkLowEnergyPP == kTRUE) && (maskIsSelected& AliVEvent::kFastOnly) ){
451 PostData(1, fListHist);
452 PostData(2, fTreeCascade);
453 return;
454 }
455 //Standard Min-Bias Selection
456 if ( ! isSelected ) {
457 PostData(1, fListHist);
458 PostData(2, fTreeCascade);
459 return;
460 }
461
462//------------------------------------------------
463// After Trigger Selection
464//------------------------------------------------
465
466 lNumberOfV0s = lESDevent->GetNumberOfV0s();
467
468 //Set variable for filling tree afterwards!
469 fHistV0MultiplicityForTrigEvt->Fill(lNumberOfV0s);
470 fHistMultiplicityForTrigEvt->Fill ( lMultiplicity );
471
472//------------------------------------------------
473// Getting: Primary Vertex + MagField Info
474//------------------------------------------------
475
476 const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();
477 // get the vtx stored in ESD found with tracks
478 lPrimaryTrackingESDVtx->GetXYZ( lTrkgPrimaryVtxPos );
479
480 const AliESDVertex *lPrimaryBestESDVtx = lESDevent->GetPrimaryVertex();
481 // get the best primary vertex available for the event
482 // As done in AliCascadeVertexer, we keep the one which is the best one available.
483 // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
484 // This one will be used for next calculations (DCA essentially)
485 lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
486
487 Double_t lPrimaryVtxPosition[3];
488 const AliVVertex *primaryVtx = lESDevent->GetPrimaryVertex();
489 lPrimaryVtxPosition[0] = primaryVtx->GetX();
490 lPrimaryVtxPosition[1] = primaryVtx->GetY();
491 lPrimaryVtxPosition[2] = primaryVtx->GetZ();
492 fHistPVx->Fill( lPrimaryVtxPosition[0] );
493 fHistPVy->Fill( lPrimaryVtxPosition[1] );
494 fHistPVz->Fill( lPrimaryVtxPosition[2] );
495
496//------------------------------------------------
497// Primary Vertex Z position: SKIP
498//------------------------------------------------
499
500 if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0 ) {
501 AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
502 PostData(1, fListHist);
503 PostData(2, fTreeCascade);
504 return;
505 }
506
507 lMagneticField = lESDevent->GetMagneticField( );
508 fHistV0MultiplicityForSelEvt ->Fill( lNumberOfV0s );
509 fHistMultiplicity->Fill(lMultiplicity);
510
511//------------------------------------------------
512// SKIP: Events with well-established PVtx
513//------------------------------------------------
514
515 const AliESDVertex *lPrimaryTrackingESDVtxCheck = lESDevent->GetPrimaryVertexTracks();
516 const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
517 if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtxCheck->GetStatus() ){
518 AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
519 PostData(1, fListHist);
520 PostData(2, fTreeCascade);
521 return;
522 }
523 fHistV0MultiplicityForSelEvtNoTPCOnly ->Fill( lNumberOfV0s );
524 fHistMultiplicityNoTPCOnly->Fill(lMultiplicity);
525
526//------------------------------------------------
527// Pileup Rejection Studies
528//------------------------------------------------
529
530 // FIXME : quality selection regarding pile-up rejection
531 if(lESDevent->IsPileupFromSPD() && !fkIsNuclear ){// minContributors=3, minZdist=0.8, nSigmaZdist=3., nSigmaDiamXY=2., nSigmaDiamZ=5. -> see http://alisoft.cern.ch/viewvc/trunk/STEER/AliESDEvent.h?root=AliRoot&r1=41914&r2=42199&pathrev=42199
532 AliWarning("Pb / This is tagged as Pileup from SPD... return !");
533 PostData(1, fListHist);
534 PostData(2, fTreeCascade);
535 return;
536 }
537 fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( lNumberOfV0s );
538 fHistMultiplicityNoTPCOnlyNoPileup->Fill(lMultiplicity);
539
540 //Do control histograms without the IsFromVertexerZ events, but consider them in analysis...
541 if( ! (lESDevent->GetPrimaryVertex()->IsFromVertexerZ() ) ){
542 fHistPVxAnalysis->Fill( lPrimaryVtxPosition[0] );
543 fHistPVyAnalysis->Fill( lPrimaryVtxPosition[1] );
544 fHistPVzAnalysis->Fill( lPrimaryVtxPosition[2] );
545 }
546
547//------------------------------------------------
548// MAIN CASCADE LOOP STARTS HERE
549//------------------------------------------------
550// Code Credit: Antonin Maire (thanks^100)
551// ---> This is an adaptation
552
553 Long_t ncascades = 0;
554 ncascades = lESDevent->GetNumberOfCascades();
555
556 for (Int_t iXi = 0; iXi < ncascades; iXi++){
557 //------------------------------------------------
558 // Initializations
559 //------------------------------------------------
560 //Double_t lTrkgPrimaryVtxRadius3D = -500.0;
561 //Double_t lBestPrimaryVtxRadius3D = -500.0;
562
563 // - 1st part of initialisation : variables needed to store AliESDCascade data members
564 Double_t lEffMassXi = 0. ;
565 //Double_t lChi2Xi = -1. ;
566 Double_t lDcaXiDaughters = -1. ;
567 Double_t lXiCosineOfPointingAngle = -1. ;
568 Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 };
569 Double_t lXiRadius = -1000. ;
570
571 // - 2nd part of initialisation : Nbr of clusters within TPC for the 3 daughter cascade tracks
572 Int_t lPosTPCClusters = -1; // For ESD only ...//FIXME : wait for availability in AOD
573 Int_t lNegTPCClusters = -1; // For ESD only ...
574 Int_t lBachTPCClusters = -1; // For ESD only ...
575
576 // - 3rd part of initialisation : about V0 part in cascades
577 Double_t lInvMassLambdaAsCascDghter = 0.;
578 //Double_t lV0Chi2Xi = -1. ;
579 Double_t lDcaV0DaughtersXi = -1.;
580
581 Double_t lDcaBachToPrimVertexXi = -1., lDcaV0ToPrimVertexXi = -1.;
582 Double_t lDcaPosToPrimVertexXi = -1.;
583 Double_t lDcaNegToPrimVertexXi = -1.;
584 Double_t lV0CosineOfPointingAngleXi = -1. ;
585 Double_t lPosV0Xi[3] = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
586 Double_t lV0RadiusXi = -1000.0;
587 Double_t lV0quality = 0.;
588
589 // - 4th part of initialisation : Effective masses
590 Double_t lInvMassXiMinus = 0.;
591 Double_t lInvMassXiPlus = 0.;
592 Double_t lInvMassOmegaMinus = 0.;
593 Double_t lInvMassOmegaPlus = 0.;
594
595 // - 6th part of initialisation : extra info for QA
596 Double_t lXiMomX = 0. , lXiMomY = 0., lXiMomZ = 0.;
597 Double_t lXiTransvMom = 0. ;
598 Double_t lXiTransvMomMC= 0. ;
599 Double_t lXiTotMom = 0. ;
600
601 Double_t lBachMomX = 0., lBachMomY = 0., lBachMomZ = 0.;
602 //Double_t lBachTransvMom = 0.;
603 //Double_t lBachTotMom = 0.;
604
605 fTreeCascVarNegNSigmaPion = -100;
606 fTreeCascVarNegNSigmaProton = -100;
607 fTreeCascVarPosNSigmaPion = -100;
608 fTreeCascVarPosNSigmaProton = -100;
609 fTreeCascVarBachNSigmaPion = -100;
610 fTreeCascVarBachNSigmaKaon = -100;
611
612 Short_t lChargeXi = -2;
613 //Double_t lV0toXiCosineOfPointingAngle = 0. ;
614
615 Double_t lRapXi = -20.0, lRapOmega = -20.0, lRapMC = -20.0; // lEta = -20.0, lTheta = 360., lPhi = 720. ;
616 //Double_t lAlphaXi = -200., lPtArmXi = -200.0;
617
618 // -------------------------------------
619 // II.ESD - Calculation Part dedicated to Xi vertices (ESD)
620
621 AliESDcascade *xi = lESDevent->GetCascade(iXi);
622 if (!xi) continue;
623
624
625 // - II.Step 1 : around primary vertex
626 //-------------
627 //lTrkgPrimaryVtxRadius3D = TMath::Sqrt( lTrkgPrimaryVtxPos[0] * lTrkgPrimaryVtxPos[0] +
628 // lTrkgPrimaryVtxPos[1] * lTrkgPrimaryVtxPos[1] +
629 // lTrkgPrimaryVtxPos[2] * lTrkgPrimaryVtxPos[2] );
630
631 //lBestPrimaryVtxRadius3D = TMath::Sqrt( lBestPrimaryVtxPos[0] * lBestPrimaryVtxPos[0] +
632 // lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
633 // lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
634
635 // - II.Step 2 : Assigning the necessary variables for specific AliESDcascade data members (ESD)
636 //-------------
637 lV0quality = 0.;
638 xi->ChangeMassHypothesis(lV0quality , 3312); // default working hypothesis : cascade = Xi- decay
639
640 lEffMassXi = xi->GetEffMassXi();
641 //lChi2Xi = xi->GetChi2Xi();
642 lDcaXiDaughters = xi->GetDcaXiDaughters();
643 lXiCosineOfPointingAngle = xi->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0],
644 lBestPrimaryVtxPos[1],
645 lBestPrimaryVtxPos[2] );
646 // Take care : the best available vertex should be used (like in AliCascadeVertexer)
647
648 xi->GetXYZcascade( lPosXi[0], lPosXi[1], lPosXi[2] );
649 lXiRadius = TMath::Sqrt( lPosXi[0]*lPosXi[0] + lPosXi[1]*lPosXi[1] );
650
651 // - II.Step 3 : around the tracks : Bach + V0 (ESD)
652 // ~ Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
653 //-------------
654
655 UInt_t lIdxPosXi = (UInt_t) TMath::Abs( xi->GetPindex() );
656 UInt_t lIdxNegXi = (UInt_t) TMath::Abs( xi->GetNindex() );
657 UInt_t lBachIdx = (UInt_t) TMath::Abs( xi->GetBindex() );
658 // Care track label can be negative in MC production (linked with the track quality)
659 // However = normally, not the case for track index ...
660
661 // FIXME : rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
662 if(lBachIdx == lIdxNegXi) {
663 AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
664 }
665 if(lBachIdx == lIdxPosXi) {
666 AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
667 }
668
669 AliESDtrack *pTrackXi = lESDevent->GetTrack( lIdxPosXi );
670 AliESDtrack *nTrackXi = lESDevent->GetTrack( lIdxNegXi );
671 AliESDtrack *bachTrackXi = lESDevent->GetTrack( lBachIdx );
672
673 if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
674 AliWarning("ERROR: Could not retrieve one of the 3 ESD daughter tracks of the cascade ...");
675 continue;
676 }
677
678 fTreeCascVarPosEta = pTrackXi->Eta();
679 fTreeCascVarNegEta = nTrackXi->Eta();
680 fTreeCascVarBachEta = bachTrackXi->Eta();
681
682 //------------------------------------------------
683 // TPC dEdx information
684 //------------------------------------------------
685 fTreeCascVarNegNSigmaPion = fPIDResponse->NumberOfSigmasTPC( nTrackXi, AliPID::kPion );
686 fTreeCascVarNegNSigmaProton = fPIDResponse->NumberOfSigmasTPC( nTrackXi, AliPID::kProton );
687 fTreeCascVarPosNSigmaPion = fPIDResponse->NumberOfSigmasTPC( pTrackXi, AliPID::kPion );
688 fTreeCascVarPosNSigmaProton = fPIDResponse->NumberOfSigmasTPC( pTrackXi, AliPID::kProton );
689 fTreeCascVarBachNSigmaPion = fPIDResponse->NumberOfSigmasTPC( bachTrackXi, AliPID::kPion );
690 fTreeCascVarBachNSigmaKaon = fPIDResponse->NumberOfSigmasTPC( bachTrackXi, AliPID::kKaon );
691
692 //------------------------------------------------
693 // TPC Number of clusters info
694 // --- modified to save the smallest number
695 // --- of TPC clusters for the 3 tracks
696 //------------------------------------------------
697
698 lPosTPCClusters = pTrackXi->GetTPCNcls();
699 lNegTPCClusters = nTrackXi->GetTPCNcls();
700 lBachTPCClusters = bachTrackXi->GetTPCNcls();
701
702 // 1 - Poor quality related to TPCrefit
703 ULong_t pStatus = pTrackXi->GetStatus();
704 ULong_t nStatus = nTrackXi->GetStatus();
705 ULong_t bachStatus = bachTrackXi->GetStatus();
706 if ((pStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
707 if ((nStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
708 if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / Bach. track has no TPCrefit ... continue!"); continue; }
709 // 2 - Poor quality related to TPC clusters: lowest cut of 70 clusters
710 if(lPosTPCClusters < 70) { AliWarning("Pb / V0 Pos. track has less than 70 TPC clusters ... continue!"); continue; }
711 if(lNegTPCClusters < 70) { AliWarning("Pb / V0 Neg. track has less than 70 TPC clusters ... continue!"); continue; }
712 if(lBachTPCClusters < 70) { AliWarning("Pb / Bach. track has less than 70 TPC clusters ... continue!"); continue; }
713 Int_t leastnumberofclusters = 1000;
714 if( lPosTPCClusters < leastnumberofclusters ) leastnumberofclusters = lPosTPCClusters;
715 if( lNegTPCClusters < leastnumberofclusters ) leastnumberofclusters = lNegTPCClusters;
716 if( lBachTPCClusters < leastnumberofclusters ) leastnumberofclusters = lBachTPCClusters;
717
718 lInvMassLambdaAsCascDghter = xi->GetEffMass();
719 // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
720 lDcaV0DaughtersXi = xi->GetDcaV0Daughters();
721 //lV0Chi2Xi = xi->GetChi2V0();
722
723 lV0CosineOfPointingAngleXi = xi->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0],
724 lBestPrimaryVtxPos[1],
725 lBestPrimaryVtxPos[2] );
726
727 lDcaV0ToPrimVertexXi = xi->GetD( lBestPrimaryVtxPos[0],
728 lBestPrimaryVtxPos[1],
729 lBestPrimaryVtxPos[2] );
730
731 lDcaBachToPrimVertexXi = TMath::Abs( bachTrackXi->GetD( lBestPrimaryVtxPos[0],
732 lBestPrimaryVtxPos[1],
733 lMagneticField ) );
734 // Note : AliExternalTrackParam::GetD returns an algebraic value ...
735
736 xi->GetXYZ( lPosV0Xi[0], lPosV0Xi[1], lPosV0Xi[2] );
737 lV0RadiusXi = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0] + lPosV0Xi[1]*lPosV0Xi[1] );
738
739 lDcaPosToPrimVertexXi = TMath::Abs( pTrackXi ->GetD( lBestPrimaryVtxPos[0],
740 lBestPrimaryVtxPos[1],
741 lMagneticField ) );
742
743 lDcaNegToPrimVertexXi = TMath::Abs( nTrackXi ->GetD( lBestPrimaryVtxPos[0],
744 lBestPrimaryVtxPos[1],
745 lMagneticField ) );
746
747 // - II.Step 4 : around effective masses (ESD)
748 // ~ change mass hypotheses to cover all the possibilities : Xi-/+, Omega -/+
749
750 if( bachTrackXi->Charge() < 0 ) {
751 lV0quality = 0.;
752 xi->ChangeMassHypothesis(lV0quality , 3312);
753 // Calculate the effective mass of the Xi- candidate.
754 // pdg code 3312 = Xi-
755 lInvMassXiMinus = xi->GetEffMassXi();
756
757 lV0quality = 0.;
758 xi->ChangeMassHypothesis(lV0quality , 3334);
759 // Calculate the effective mass of the Xi- candidate.
760 // pdg code 3334 = Omega-
761 lInvMassOmegaMinus = xi->GetEffMassXi();
762
763 lV0quality = 0.;
764 xi->ChangeMassHypothesis(lV0quality , 3312); // Back to default hyp.
765 }// end if negative bachelor
766
767
768 if( bachTrackXi->Charge() > 0 ){
769 lV0quality = 0.;
770 xi->ChangeMassHypothesis(lV0quality , -3312);
771 // Calculate the effective mass of the Xi+ candidate.
772 // pdg code -3312 = Xi+
773 lInvMassXiPlus = xi->GetEffMassXi();
774
775 lV0quality = 0.;
776 xi->ChangeMassHypothesis(lV0quality , -3334);
777 // Calculate the effective mass of the Xi+ candidate.
778 // pdg code -3334 = Omega+
779 lInvMassOmegaPlus = xi->GetEffMassXi();
780
781 lV0quality = 0.;
782 xi->ChangeMassHypothesis(lV0quality , -3312); // Back to "default" hyp.
783 }// end if positive bachelor
784 // - II.Step 6 : extra info for QA (ESD)
785 // miscellaneous pieces of info that may help regarding data quality assessment.
786 //-------------
787
788 xi->GetPxPyPz( lXiMomX, lXiMomY, lXiMomZ );
789 lXiTransvMom = TMath::Sqrt( lXiMomX*lXiMomX + lXiMomY*lXiMomY );
790 lXiTotMom = TMath::Sqrt( lXiMomX*lXiMomX + lXiMomY*lXiMomY + lXiMomZ*lXiMomZ );
791
792 xi->GetBPxPyPz( lBachMomX, lBachMomY, lBachMomZ );
793 //lBachTransvMom = TMath::Sqrt( lBachMomX*lBachMomX + lBachMomY*lBachMomY );
794 //lBachTotMom = TMath::Sqrt( lBachMomX*lBachMomX + lBachMomY*lBachMomY + lBachMomZ*lBachMomZ );
795
796 lChargeXi = xi->Charge();
797
798 //lV0toXiCosineOfPointingAngle = xi->GetV0CosineOfPointingAngle( lPosXi[0], lPosXi[1], lPosXi[2] );
799
800 lRapXi = xi->RapXi();
801 lRapOmega = xi->RapOmega();
802 //lEta = xi->Eta();
803 //lTheta = xi->Theta() *180.0/TMath::Pi();
804 //lPhi = xi->Phi() *180.0/TMath::Pi();
805 //lAlphaXi = xi->AlphaXi();
806 //lPtArmXi = xi->PtArmXi();
807
808 //------------------------------------------------
809 // Set Variables for adding to tree
810 //------------------------------------------------
811
812/* 1*/ fTreeCascVarCharge = lChargeXi;
813/* 2*/ if(lInvMassXiMinus!=0) fTreeCascVarMassAsXi = lInvMassXiMinus;
814/* 2*/ if(lInvMassXiPlus!=0) fTreeCascVarMassAsXi = lInvMassXiPlus;
815/* 3*/ if(lInvMassOmegaMinus!=0) fTreeCascVarMassAsOmega = lInvMassOmegaMinus;
816/* 3*/ if(lInvMassOmegaPlus!=0) fTreeCascVarMassAsOmega = lInvMassOmegaPlus;
817/* 4*/ fTreeCascVarPt = lXiTransvMom;
818/* 4*/ fTreeCascVarPtMC = lXiTransvMomMC;
819/* 5*/ fTreeCascVarRapXi = lRapXi ;
820/* 5*/ fTreeCascVarRapMC = lRapMC ;
821/* 6*/ fTreeCascVarRapOmega = lRapOmega ;
822/* 7*/ fTreeCascVarDCACascDaughters = lDcaXiDaughters;
823/* 8*/ fTreeCascVarDCABachToPrimVtx = lDcaBachToPrimVertexXi;
824/* 9*/ fTreeCascVarDCAV0Daughters = lDcaV0DaughtersXi;
825/*10*/ fTreeCascVarDCAV0ToPrimVtx = lDcaV0ToPrimVertexXi;
826/*11*/ fTreeCascVarDCAPosToPrimVtx = lDcaPosToPrimVertexXi;
827/*12*/ fTreeCascVarDCANegToPrimVtx = lDcaNegToPrimVertexXi;
828/*13*/ fTreeCascVarCascCosPointingAngle = lXiCosineOfPointingAngle;
829/*14*/ fTreeCascVarCascRadius = lXiRadius;
830/*15*/ fTreeCascVarV0Mass = lInvMassLambdaAsCascDghter;
831/*16*/ fTreeCascVarV0CosPointingAngle = lV0CosineOfPointingAngleXi;
832/*17*/ fTreeCascVarV0Radius = lV0RadiusXi;
833/*20*/ fTreeCascVarLeastNbrClusters = leastnumberofclusters;
834/*21*/ fTreeCascVarMultiplicity = lMultiplicity; //multiplicity, whatever that may be
835
836/*23*/ fTreeCascVarDistOverTotMom = TMath::Sqrt(
837 TMath::Power( lPosXi[0] - lBestPrimaryVtxPos[0] , 2) +
838 TMath::Power( lPosXi[1] - lBestPrimaryVtxPos[1] , 2) +
839 TMath::Power( lPosXi[2] - lBestPrimaryVtxPos[2] , 2)
840 );
841/*23*/ fTreeCascVarDistOverTotMom /= (lXiTotMom+1e-13);
842
843//All vars not specified here: specified elsewhere!
844
845//------------------------------------------------
846// Fill Tree!
847//------------------------------------------------
848
849// The conditional is meant to decrease excessive
850// memory usage! Be careful when loosening the
851// cut!
852
853 //Xi Mass window: 150MeV wide
854 //Omega mass window: 150MeV wide
855
856 if( (fTreeCascVarMassAsXi<1.32+0.075&&fTreeCascVarMassAsXi>1.32-0.075) ||
857 (fTreeCascVarMassAsOmega<1.68+0.075&&fTreeCascVarMassAsOmega>1.68-0.075) ){
858 fTreeCascade->Fill();
859 }
860
861//------------------------------------------------
862// Fill tree over.
863//------------------------------------------------
864
865 }// end of the Cascade loop (ESD or AOD)
866
867 // Post output data.
868 PostData(1, fListHist);
869 PostData(2, fTreeCascade);
870}
871
872//________________________________________________________________________
873void AliAnalysisTaskExtractCascade::Terminate(Option_t *)
874{
875 // Draw result to the screen
876 // Called once at the end of the query
877
878 TList *cRetrievedList = 0x0;
879 cRetrievedList = (TList*)GetOutputData(1);
880 if(!cRetrievedList){
881 Printf("ERROR - AliAnalysisTaskExtractCascade : ouput data container list not available\n");
882 return;
883 }
884
885 fHistV0MultiplicityForTrigEvt = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistV0MultiplicityForTrigEvt") );
886 if (!fHistV0MultiplicityForTrigEvt) {
887 Printf("ERROR - AliAnalysisTaskExtractCascade : fHistV0MultiplicityForTrigEvt not available");
888 return;
889 }
890
891 TCanvas *canCheck = new TCanvas("AliAnalysisTaskExtractCascade","V0 Multiplicity",10,10,510,510);
892 canCheck->cd(1)->SetLogy();
893
894 fHistV0MultiplicityForTrigEvt->SetMarkerStyle(22);
895 fHistV0MultiplicityForTrigEvt->DrawCopy("E");
896}
897
898//----------------------------------------------------------------------------
899
900Double_t AliAnalysisTaskExtractCascade::MyRapidity(Double_t rE, Double_t rPz) const
901{
902 // Local calculation for rapidity
903 Double_t ReturnValue = -100;
904 if( (rE-rPz+1.e-13) != 0 && (rE+rPz) != 0 ){
905 ReturnValue = 0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
906 }
907 return ReturnValue;
908}