]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/STRANGENESS/LambdaK0/AliAnalysisTaskExtractV0AOD.cxx
Modifications to allow for multiplicity-dependent transverse momentum spectra for the
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / LambdaK0 / AliAnalysisTaskExtractV0AOD.cxx
CommitLineData
a0879a26 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//
23// --- Added bits of code for checking V0s
24// (from AliAnalysisTaskCheckStrange.cxx)
25//
26// --- Algorithm Description
27// 1. Perform Physics Selection
28// 2. Perform Primary Vertex |z|<10cm selection
29// 3. Perform Primary Vertex NoTPCOnly vertexing selection (>0 contrib.)
30// 4. Perform Pileup Rejection
31// 5. Analysis Loops:
32// 5a. Fill TTree object with V0 information, candidates
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 "THnSparse.h"
60#include "TVector3.h"
61#include "TCanvas.h"
62#include "TMath.h"
63#include "TLegend.h"
64#include "AliLog.h"
65#include "AliCentrality.h"
66#include "AliESDEvent.h"
67#include "AliAODEvent.h"
68#include "AliV0vertexer.h"
69#include "AliCascadeVertexer.h"
70#include "AliESDpid.h"
71#include "AliESDtrack.h"
72
73#include "AliESDtrackCuts.h"
74#include "AliInputEventHandler.h"
75#include "AliAnalysisManager.h"
76#include "AliMCEventHandler.h"
77
78#include "AliCFContainer.h"
79#include "AliMultiplicity.h"
80
81#include "AliESDcascade.h"
82#include "AliAODcascade.h"
83#include "AliESDUtils.h"
84#include "AliESDHeader.h"
85#include "AliAODTrack.h"
86#include "AliAnalysisTaskExtractV0AOD.h"
87
88//debugging purposes
89#include "TObjectTable.h"
90
91ClassImp(AliAnalysisTaskExtractV0AOD)
92
93AliAnalysisTaskExtractV0AOD::AliAnalysisTaskExtractV0AOD()
94 : AliAnalysisTaskSE(), fListHistV0(0), fTree(0), fPIDResponse(0),
95 fkIsNuclear ( kFALSE ),
96 fkLowEnergyPP ( kFALSE ),
97 fkUseOnTheFly ( kFALSE ),
98
99//------------------------------------------------
100// HISTOGRAMS
101// --- Filled on an Event-by-event basis
102//------------------------------------------------
103 fHistV0MultiplicityBeforeTrigSel(0),
104 fHistV0MultiplicityForTrigEvt(0),
105 fHistV0MultiplicityForSelEvt(0),
106 fHistV0MultiplicityForSelEvtNoTPCOnly(0),
107 fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
108 fHistMultiplicityBeforeTrigSel(0),
109 fHistMultiplicityForTrigEvt(0),
110 fHistMultiplicity(0),
111 fHistMultiplicityNoTPCOnly(0),
112 fHistMultiplicityNoTPCOnlyNoPileup(0),
113 fHistPVx(0),
114 fHistPVy(0),
115 fHistPVz(0),
116 fHistPVxAnalysis(0),
117 fHistPVyAnalysis(0),
118 fHistPVzAnalysis(0),
119 fHistSwappedV0Counter(0)
120{
121 // Dummy Constructor
122}
123
124AliAnalysisTaskExtractV0AOD::AliAnalysisTaskExtractV0AOD(const char *name)
125 : AliAnalysisTaskSE(name), fListHistV0(0), fTree(0), fPIDResponse(0),
126 fkIsNuclear ( kFALSE ),
127 fkLowEnergyPP ( kFALSE ),
128 fkUseOnTheFly ( kFALSE ),
129
130//------------------------------------------------
131// HISTOGRAMS
132// --- Filled on an Event-by-event basis
133//------------------------------------------------
134 fHistV0MultiplicityBeforeTrigSel(0),
135 fHistV0MultiplicityForTrigEvt(0),
136 fHistV0MultiplicityForSelEvt(0),
137 fHistV0MultiplicityForSelEvtNoTPCOnly(0),
138 fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
139 fHistMultiplicityBeforeTrigSel(0),
140 fHistMultiplicityForTrigEvt(0),
141 fHistMultiplicity(0),
142 fHistMultiplicityNoTPCOnly(0),
143 fHistMultiplicityNoTPCOnlyNoPileup(0),
144 fHistPVx(0),
145 fHistPVy(0),
146 fHistPVz(0),
147 fHistPVxAnalysis(0),
148 fHistPVyAnalysis(0),
149 fHistPVzAnalysis(0),
150 fHistSwappedV0Counter(0)
151{
152 // Constructor
153 // Output slot #0 writes into a TList container (Lambda Histos and fTree)
154 DefineOutput(1, TList::Class());
155 DefineOutput(2, TTree::Class());
156}
157
158
159AliAnalysisTaskExtractV0AOD::~AliAnalysisTaskExtractV0AOD()
160{
161//------------------------------------------------
162// DESTRUCTOR
163//------------------------------------------------
164
165 if (fListHistV0){
166 delete fListHistV0;
167 fListHistV0 = 0x0;
168 }
169 if (fTree){
170 delete fTree;
171 fTree = 0x0;
172 }
173}
174
175
176
177//________________________________________________________________________
178void AliAnalysisTaskExtractV0AOD::UserCreateOutputObjects()
179{
180
181 //Create File-resident Tree, please.
182 OpenFile(2);
183 // Called once
184 fTree = new TTree("fTree","V0Candidates");
185
186//------------------------------------------------
187// fTree Branch definitions
188//------------------------------------------------
189
190//-----------BASIC-INFO---------------------------
191/* 1*/ fTree->Branch("fTreeVariableChi2V0",&fTreeVariableChi2V0,"fTreeVariableChi2V0/F");
192/* 2*/ fTree->Branch("fTreeVariableDcaV0Daughters",&fTreeVariableDcaV0Daughters,"fTreeVariableDcaV0Daughters/F");
193/* 3*/ fTree->Branch("fTreeVariableDcaPosToPrimVertex",&fTreeVariableDcaPosToPrimVertex,"fTreeVariableDcaPosToPrimVertex/F");
194/* 4*/ fTree->Branch("fTreeVariableDcaNegToPrimVertex",&fTreeVariableDcaNegToPrimVertex,"fTreeVariableDcaNegToPrimVertex/F");
195/* 5*/ fTree->Branch("fTreeVariableV0Radius",&fTreeVariableV0Radius,"fTreeVariableV0Radius/F");
196/* 6*/ fTree->Branch("fTreeVariablePt",&fTreeVariablePt,"fTreeVariablePt/F");
197/* 7*/ fTree->Branch("fTreeVariableRapK0Short",&fTreeVariableRapK0Short,"fTreeVariableRapK0Short/F");
198/* 8*/ fTree->Branch("fTreeVariableRapLambda",&fTreeVariableRapLambda,"fTreeVariableRapLambda/F");
199/* 9*/ fTree->Branch("fTreeVariableInvMassK0s",&fTreeVariableInvMassK0s,"fTreeVariableInvMassK0s/F");
200/*10*/ fTree->Branch("fTreeVariableInvMassLambda",&fTreeVariableInvMassLambda,"fTreeVariableInvMassLambda/F");
201/*11*/ fTree->Branch("fTreeVariableInvMassAntiLambda",&fTreeVariableInvMassAntiLambda,"fTreeVariableInvMassAntiLambda/F");
202/*12*/ fTree->Branch("fTreeVariableV0CosineOfPointingAngle",&fTreeVariableV0CosineOfPointingAngle,"fTreeVariableV0CosineOfPointingAngle/F");
203/*13*/ fTree->Branch("fTreeVariableAlphaV0",&fTreeVariableAlphaV0,"fTreeVariableAlphaV0/F");
204/*14*/ fTree->Branch("fTreeVariablePtArmV0",&fTreeVariablePtArmV0,"fTreeVariablePtArmV0/F");
205/*15*/ fTree->Branch("fTreeVariableLeastNbrCrossedRows",&fTreeVariableLeastNbrCrossedRows,"fTreeVariableLeastNbrCrossedRows/I");
206/*16*/ fTree->Branch("fTreeVariableLeastRatioCrossedRowsOverFindable",&fTreeVariableLeastRatioCrossedRowsOverFindable,"fTreeVariableLeastRatioCrossedRowsOverFindable/F");
207//-----------MULTIPLICITY-INFO--------------------
208/*17*/ fTree->Branch("fTreeVariableMultiplicity",&fTreeVariableMultiplicity,"fTreeVariableMultiplicity/I");
209//------------------------------------------------
210/*18*/ fTree->Branch("fTreeVariableDistOverTotMom",&fTreeVariableDistOverTotMom,"fTreeVariableDistOverTotMom/F");
211/*19*/ fTree->Branch("fTreeVariableNSigmasPosProton",&fTreeVariableNSigmasPosProton,"fTreeVariableNSigmasPosProton/F");
212/*20*/ fTree->Branch("fTreeVariableNSigmasPosPion",&fTreeVariableNSigmasPosPion,"fTreeVariableNSigmasPosPion/F");
213/*21*/ fTree->Branch("fTreeVariableNSigmasNegProton",&fTreeVariableNSigmasNegProton,"fTreeVariableNSigmasNegProton/F");
214/*22*/ fTree->Branch("fTreeVariableNSigmasNegPion",&fTreeVariableNSigmasNegPion,"fTreeVariableNSigmasNegPion/F");
215/*23*/ fTree->Branch("fTreeVariableNegEta",&fTreeVariableNegEta,"fTreeVariableNegEta/F");
216/*24*/ fTree->Branch("fTreeVariablePosEta",&fTreeVariablePosEta,"fTreeVariablePosEta/F");
217/*25*/ fTree->Branch("fTreeVariableRunNumber",&fTreeVariableRunNumber,"fTreeVariableRunNumber/I");
218/*26*/ fTree->Branch("fTreeVariableEventNumber",&fTreeVariableEventNumber,"fTreeVariableEventNumber/l");
219
220//------------------------------------------------
221// Particle Identification Setup
222//------------------------------------------------
223
224 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
225 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
226 fPIDResponse = inputHandler->GetPIDResponse();
227
228//Skipped. Will use Local setup.
229
230//------------------------------------------------
231// V0 Multiplicity Histograms
232//------------------------------------------------
233
234 // Create histograms
235 //Create File-resident Tree, please.
236 OpenFile(1);
237
238 fListHistV0 = new TList();
239 fListHistV0->SetOwner(); // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
240
241 if(! fHistV0MultiplicityBeforeTrigSel) {
242 fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel",
243 "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events",
244 25, 0, 25);
245 fListHistV0->Add(fHistV0MultiplicityBeforeTrigSel);
246 }
247
248 if(! fHistV0MultiplicityForTrigEvt) {
249 fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt",
250 "V0s per event (for triggered evt);Nbr of V0s/Evt;Events",
251 25, 0, 25);
252 fListHistV0->Add(fHistV0MultiplicityForTrigEvt);
253 }
254
255 if(! fHistV0MultiplicityForSelEvt) {
256 fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt",
257 "V0s per event;Nbr of V0s/Evt;Events",
258 25, 0, 25);
259 fListHistV0->Add(fHistV0MultiplicityForSelEvt);
260 }
261
262 if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {
263 fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly",
264 "V0s per event;Nbr of V0s/Evt;Events",
265 25, 0, 25);
266 fListHistV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);
267 }
268
269 if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {
270 fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup",
271 "V0s per event;Nbr of V0s/Evt;Events",
272 25, 0, 25);
273 fListHistV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);
274 }
275
276//------------------------------------------------
277// Track Multiplicity Histograms
278//------------------------------------------------
279
280 if(! fHistMultiplicityBeforeTrigSel) {
281 fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel",
282 "Tracks per event;Nbr of Tracks;Events",
283 200, 0, 200);
284 fListHistV0->Add(fHistMultiplicityBeforeTrigSel);
285 }
286 if(! fHistMultiplicityForTrigEvt) {
287 fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt",
288 "Tracks per event;Nbr of Tracks;Events",
289 200, 0, 200);
290 fListHistV0->Add(fHistMultiplicityForTrigEvt);
291 }
292 if(! fHistMultiplicity) {
293 fHistMultiplicity = new TH1F("fHistMultiplicity",
294 "Tracks per event;Nbr of Tracks;Events",
295 200, 0, 200);
296 fListHistV0->Add(fHistMultiplicity);
297 }
298 if(! fHistMultiplicityNoTPCOnly) {
299 fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly",
300 "Tracks per event;Nbr of Tracks;Events",
301 200, 0, 200);
302 fListHistV0->Add(fHistMultiplicityNoTPCOnly);
303 }
304 if(! fHistMultiplicityNoTPCOnlyNoPileup) {
305 fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup",
306 "Tracks per event;Nbr of Tracks;Events",
307 200, 0, 200);
308 fListHistV0->Add(fHistMultiplicityNoTPCOnlyNoPileup);
309 }
310
311 if(! fHistPVx) {
312 fHistPVx = new TH1F("fHistPVx",
313 "PV x position;Nbr of Evts;x",
314 2000, -0.5, 0.5);
315 fListHistV0->Add(fHistPVx);
316 }
317 if(! fHistPVy) {
318 fHistPVy = new TH1F("fHistPVy",
319 "PV y position;Nbr of Evts;y",
320 2000, -0.5, 0.5);
321 fListHistV0->Add(fHistPVy);
322 }
323 if(! fHistPVz) {
324 fHistPVz = new TH1F("fHistPVz",
325 "PV z position;Nbr of Evts;z",
326 400, -20, 20);
327 fListHistV0->Add(fHistPVz);
328 }
329 if(! fHistPVxAnalysis) {
330 fHistPVxAnalysis = new TH1F("fHistPVxAnalysis",
331 "PV x position;Nbr of Evts;x",
332 2000, -0.5, 0.5);
333 fListHistV0->Add(fHistPVxAnalysis);
334 }
335 if(! fHistPVyAnalysis) {
336 fHistPVyAnalysis = new TH1F("fHistPVyAnalysis",
337 "PV y position;Nbr of Evts;y",
338 2000, -0.5, 0.5);
339 fListHistV0->Add(fHistPVyAnalysis);
340 }
341 if(! fHistPVzAnalysis) {
342 fHistPVzAnalysis = new TH1F("fHistPVzAnalysis",
343 "PV z position;Nbr of Evts;z",
344 400, -20, 20);
345 fListHistV0->Add(fHistPVzAnalysis);
346 }
347 if(! fHistSwappedV0Counter) {
348 fHistSwappedV0Counter = new TH1F("fHistSwappedV0Counter",
349 "Swap or not histo;Swapped (1) or not (0); count",
350 2, 0, 2);
351 fListHistV0->Add(fHistSwappedV0Counter);
352 }
353 //Regular output: Histograms
354 PostData(1, fListHistV0);
355 //TTree Object: Saved to base directory. Should cache to disk while saving.
356 //(Important to avoid excessive memory usage, particularly when merging)
357 PostData(2, fTree);
358
359}// end UserCreateOutputObjects
360
361
362//________________________________________________________________________
363void AliAnalysisTaskExtractV0AOD::UserExec(Option_t *)
364{
365
366 // Main loop
367 // Called for each event
368 //gObjectTable->Print();
369 AliAODEvent *lAODevent = 0x0;
370
371 //AliAODEvent *lAODevent = 0x0;
372 Int_t nV0s = -1;
373
374 Double_t lBestPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
375 Double_t lMagneticField = -10.;
376
377 // Connect to the InputEvent
378 // After these lines, we should have an ESD/AOD event + the number of cascades in it.
379
380 // Appropriate for ESD analysis!
381
382 lAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
383 if (!lAODevent) {
384 AliWarning("ERROR: lAODevent not available \n");
385 return;
386 }
387 fTreeVariableRunNumber = lAODevent->GetRunNumber();
388 fTreeVariableEventNumber =
389 ( ( ((ULong64_t)lAODevent->GetPeriodNumber() ) << 36 ) |
390 ( ((ULong64_t)lAODevent->GetOrbitNumber () ) << 12 ) |
391 ((ULong64_t)lAODevent->GetBunchCrossNumber() ) );
392
393 //REVISED multiplicity estimator after 'multiplicity day' (2011)
394 Int_t lMultiplicity = -100;
395
396 if(fkIsNuclear == kFALSE) lMultiplicity = -1;
397
398 //---> If this is a nuclear collision, then go nuclear on "multiplicity" variable...
399 //---> Warning: Experimental
400 if(fkIsNuclear == kTRUE){
401 AliCentrality* centrality;
402 centrality = lAODevent->GetCentrality();
403 lMultiplicity = ( ( Int_t ) ( centrality->GetCentralityPercentile( "V0M" ) ) );
404 if (centrality->GetQuality()>1) {
405 PostData(1, fListHistV0);
406 PostData(2, fTree);
407 return;
408 }
409 }
410
411 //Set variable for filling tree afterwards!
412 //---> pp case......: GetReferenceMultiplicity
413 //---> Pb-Pb case...: Centrality by V0M
414 fTreeVariableMultiplicity = lMultiplicity;
415 fHistMultiplicityBeforeTrigSel->Fill ( lMultiplicity );
416 fHistV0MultiplicityBeforeTrigSel->Fill ( lAODevent->GetNumberOfV0s() );
417
418//------------------------------------------------
419// Physics Selection
420//------------------------------------------------
421
422// new method
423 UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
424 Bool_t isSelected = 0;
425 isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
426
427 //pp at 2.76TeV: special case, ignore FastOnly
428 if ( (fkLowEnergyPP == kTRUE) && (maskIsSelected& AliVEvent::kFastOnly) ){
429 PostData(1, fListHistV0);
430 PostData(2, fTree);
431 return;
432 }
433 //Standard Min-Bias Selection
434 if ( ! isSelected ) {
435 PostData(1, fListHistV0);
436 PostData(2, fTree);
437 return;
438 }
439
440//------------------------------------------------
441// After Trigger Selection
442//------------------------------------------------
443
444 nV0s = lAODevent->GetNumberOfV0s();
445
446 fHistV0MultiplicityForTrigEvt ->Fill( nV0s );
447 fHistMultiplicityForTrigEvt ->Fill( lMultiplicity );
448
449//------------------------------------------------
450// Getting: Primary Vertex + MagField Info
451//------------------------------------------------
452
453 const AliAODVertex *lPrimaryBestAODVtx = lAODevent->GetPrimaryVertex();
454 // get the best primary vertex available for the event
455 // As done in AliCascadeVertexer, we keep the one which is the best one available.
456 // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
457 // This one will be used for next calculations (DCA essentially)
458 lPrimaryBestAODVtx->GetXYZ( lBestPrimaryVtxPos );
459
460 Double_t tPrimaryVtxPosition[3];
461 const AliVVertex *primaryVtx = lAODevent->GetPrimaryVertex();
462 tPrimaryVtxPosition[0] = primaryVtx->GetX();
463 tPrimaryVtxPosition[1] = primaryVtx->GetY();
464 tPrimaryVtxPosition[2] = primaryVtx->GetZ();
465 fHistPVx->Fill( tPrimaryVtxPosition[0] );
466 fHistPVy->Fill( tPrimaryVtxPosition[1] );
467 fHistPVz->Fill( tPrimaryVtxPosition[2] );
468
469//------------------------------------------------
470// Primary Vertex Z position: SKIP
471//------------------------------------------------
472
473 if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0 ) {
474 AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
475 PostData(1, fListHistV0);
476 PostData(2, fTree);
477 return;
478 }
479
480 lMagneticField = lAODevent->GetMagneticField( );
481 fHistV0MultiplicityForSelEvt ->Fill( nV0s );
482 fHistMultiplicity->Fill(lMultiplicity);
483
484//------------------------------------------------
485// Only look at events with well-established PV
486//------------------------------------------------
487
488 const AliAODVertex *lPrimaryTrackingAODVtxCheck = lAODevent->GetPrimaryVertex();
489 const AliAODVertex *lPrimarySPDVtx = lAODevent->GetPrimaryVertexSPD();
490 if (!lPrimarySPDVtx && !lPrimaryTrackingAODVtxCheck ){
491 AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
492 PostData(1, fListHistV0);
493 PostData(2, fTree);
494 return;
495 }
496
497 fHistV0MultiplicityForSelEvtNoTPCOnly ->Fill( nV0s );
498 fHistMultiplicityNoTPCOnly->Fill(lMultiplicity);
499
500//------------------------------------------------
501// Pileup Rejection
502//------------------------------------------------
503
504 // FIXME : quality selection regarding pile-up rejection
505 if(lAODevent->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
506 PostData(1, fListHistV0);
507 PostData(2, fTree);
508 return;
509 }
510 fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( nV0s );
511 fHistMultiplicityNoTPCOnlyNoPileup->Fill(lMultiplicity);
512
513//------------------------------------------------
514// MAIN LAMBDA LOOP STARTS HERE
515//------------------------------------------------
516
517 //if( ! (lAODevent->GetPrimaryVertex()->IsFromVertexerZ() ) ){
518 fHistPVxAnalysis->Fill( tPrimaryVtxPosition[0] );
519 fHistPVyAnalysis->Fill( tPrimaryVtxPosition[1] );
520 fHistPVzAnalysis->Fill( tPrimaryVtxPosition[2] );
521 //}
522
523//Variable definition
524 Int_t lOnFlyStatus = 0;// nv0sOn = 0, nv0sOff = 0;
525 Double_t lChi2V0 = 0;
526 Double_t lDcaV0Daughters = 0, lDcaV0ToPrimVertex = 0;
527 Double_t lDcaPosToPrimVertex = 0, lDcaNegToPrimVertex = 0;
528 Double_t lV0CosineOfPointingAngle = 0;
529 Double_t lV0Radius = 0, lPt = 0;
530 Double_t lRapK0Short = 0, lRapLambda = 0;
531 Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
532 Double_t lAlphaV0 = 0, lPtArmV0 = 0;
533
534 Double_t fMinV0Pt = 0;
535 Double_t fMaxV0Pt = 100;
536
537 Int_t nv0s = 0;
538 nv0s = lAODevent->GetNumberOfV0s();
539
540 //for (Int_t iV0 = 0; iV0 < nv0s; iV0++)
541 for (Int_t iV0 = 0; iV0 < nv0s; iV0++) //extra-crazy test
542 {// This is the begining of the V0 loop
543 AliAODv0 *v0 = lAODevent->GetV0(iV0);
544 if (!v0) continue;
545
546 //Obsolete at AOD level...
547 //---> Fix On-the-Fly candidates, count how many swapped
548 //if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ){
549 // fHistSwappedV0Counter -> Fill( 1 );
550 //}else{
551 // fHistSwappedV0Counter -> Fill( 0 );
552 //}
553 //if ( fkUseOnTheFly ) CheckChargeV0(v0);
554
555 Double_t tDecayVertexV0[3]; v0->GetXYZ(tDecayVertexV0);
556 Double_t tV0mom[3];
557 v0->GetPxPyPz( tV0mom );
558 Double_t lV0TotalMomentum = TMath::Sqrt(
559 tV0mom[0]*tV0mom[0]+tV0mom[1]*tV0mom[1]+tV0mom[2]*tV0mom[2] );
560
561 lV0Radius = TMath::Sqrt(tDecayVertexV0[0]*tDecayVertexV0[0]+tDecayVertexV0[1]*tDecayVertexV0[1]);
562 lPt = v0->Pt();
563 lRapK0Short = v0->RapK0Short();
564 lRapLambda = v0->RapLambda();
565 if ((lPt<fMinV0Pt)||(fMaxV0Pt<lPt)) continue;
566
567 UInt_t lKeyPos = (UInt_t)TMath::Abs(v0->GetPosID());
568 UInt_t lKeyNeg = (UInt_t)TMath::Abs(v0->GetPosID());
569
570 Double_t lMomPos[3]; //v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
571 Double_t lMomNeg[3]; //v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
572 lMomPos[0] = v0->MomPosX();
573 lMomPos[1] = v0->MomPosY();
574 lMomPos[2] = v0->MomPosZ();
575 lMomNeg[0] = v0->MomNegX();
576 lMomNeg[1] = v0->MomNegY();
577 lMomNeg[2] = v0->MomNegZ();
578
579 AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
580 AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
581 if (!pTrack || !nTrack) {
582 Printf("ERROR: Could not retreive one of the daughter track");
583 continue;
584 }
585
586 //Daughter Eta for Eta selection, afterwards
587 fTreeVariableNegEta = nTrack->Eta();
588 fTreeVariablePosEta = pTrack->Eta();
589
590 // Filter like-sign V0 (next: add counter and distribution)
591 if ( pTrack->Charge() == nTrack->Charge()){
592 continue;
593 }
594
595 //Quick test this far!
596
597
598 //________________________________________________________________________
599 // Track quality cuts
600 Float_t lPosTrackCrossedRows = pTrack->GetTPCClusterInfo(2,1);
601 Float_t lNegTrackCrossedRows = nTrack->GetTPCClusterInfo(2,1);
602 fTreeVariableLeastNbrCrossedRows = (Int_t) lPosTrackCrossedRows;
603 if( lNegTrackCrossedRows < fTreeVariableLeastNbrCrossedRows )
604 fTreeVariableLeastNbrCrossedRows = (Int_t) lNegTrackCrossedRows;
605
606 // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)
607 if( !(pTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
608 if( !(nTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
609
610 if ( ( ( pTrack->GetTPCClusterInfo(2,1) ) < 70 ) || ( ( nTrack->GetTPCClusterInfo(2,1) ) < 70 ) ) continue;
611
612 //GetKinkIndex condition
613 //if( pTrack->GetKinkIndex(0)>0 || nTrack->GetKinkIndex(0)>0 ) continue;
614
615 //Findable clusters > 0 condition
616 if( pTrack->GetTPCNclsF()<=0 || nTrack->GetTPCNclsF()<=0 ) continue;
617
618 //Compute ratio Crossed Rows / Findable clusters
619 //Note: above test avoids division by zero!
620 Float_t lPosTrackCrossedRowsOverFindable = lPosTrackCrossedRows / ((double)(pTrack->GetTPCNclsF()));
621 Float_t lNegTrackCrossedRowsOverFindable = lNegTrackCrossedRows / ((double)(nTrack->GetTPCNclsF()));
622
623 fTreeVariableLeastRatioCrossedRowsOverFindable = lPosTrackCrossedRowsOverFindable;
624 if( lNegTrackCrossedRowsOverFindable < fTreeVariableLeastRatioCrossedRowsOverFindable )
625 fTreeVariableLeastRatioCrossedRowsOverFindable = lNegTrackCrossedRowsOverFindable;
626
627 //Lowest Cut Level for Ratio Crossed Rows / Findable = 0.8, set here
628 if ( fTreeVariableLeastRatioCrossedRowsOverFindable < 0.8 ) continue;
629
630 //End track Quality Cuts
631 //________________________________________________________________________
632
633
634 //ESD code: obsolete
635 //lDcaPosToPrimVertex = TMath::Abs(pTrack->GetD(tPrimaryVtxPosition[0],
636 // tPrimaryVtxPosition[1],
637 // lMagneticField) );
638
639 //lDcaNegToPrimVertex = TMath::Abs(nTrack->GetD(tPrimaryVtxPosition[0],
640 // tPrimaryVtxPosition[1],
641 // lMagneticField) );
642
643 lDcaPosToPrimVertex = v0->DcaPosToPrimVertex();
644 lDcaNegToPrimVertex = v0->DcaNegToPrimVertex();
645
646
647 lOnFlyStatus = v0->GetOnFlyStatus();
648 lChi2V0 = v0->Chi2V0();
649 lDcaV0Daughters = v0->DcaV0Daughters();
650 lDcaV0ToPrimVertex = v0->DcaV0ToPrimVertex();
651 lV0CosineOfPointingAngle = v0->CosPointingAngle(tPrimaryVtxPosition);
652 fTreeVariableV0CosineOfPointingAngle=lV0CosineOfPointingAngle;
653
654 // Getting invariant mass infos directly from ESD
655 lInvMassK0s = v0->MassK0Short();
656 lInvMassLambda = v0->MassLambda();
657 lInvMassAntiLambda = v0->MassAntiLambda();
658 lAlphaV0 = v0->AlphaV0();
659 lPtArmV0 = v0->PtArmV0();
660
661 fTreeVariablePt = v0->Pt();
662 fTreeVariableChi2V0 = lChi2V0;
663 fTreeVariableDcaV0ToPrimVertex = lDcaV0ToPrimVertex;
664 fTreeVariableDcaV0Daughters = lDcaV0Daughters;
665 fTreeVariableV0CosineOfPointingAngle = lV0CosineOfPointingAngle;
666 fTreeVariableV0Radius = lV0Radius;
667 fTreeVariableDcaPosToPrimVertex = lDcaPosToPrimVertex;
668 fTreeVariableDcaNegToPrimVertex = lDcaNegToPrimVertex;
669 fTreeVariableInvMassK0s = lInvMassK0s;
670 fTreeVariableInvMassLambda = lInvMassLambda;
671 fTreeVariableInvMassAntiLambda = lInvMassAntiLambda;
672 fTreeVariableRapK0Short = lRapK0Short;
673 fTreeVariableRapLambda = lRapLambda;
674 fTreeVariableAlphaV0 = lAlphaV0;
675 fTreeVariablePtArmV0 = lPtArmV0;
676
677 //Official means of acquiring N-sigmas
678 fTreeVariableNSigmasPosProton = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kProton );
679 fTreeVariableNSigmasPosPion = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kPion );
680 fTreeVariableNSigmasNegProton = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kProton );
681 fTreeVariableNSigmasNegPion = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kPion );
682
683//This requires an Invariant Mass Hypothesis afterwards
684 fTreeVariableDistOverTotMom = TMath::Sqrt(
685 TMath::Power( tDecayVertexV0[0] - lBestPrimaryVtxPos[0] , 2) +
686 TMath::Power( tDecayVertexV0[1] - lBestPrimaryVtxPos[1] , 2) +
687 TMath::Power( tDecayVertexV0[2] - lBestPrimaryVtxPos[2] , 2)
688 );
689 fTreeVariableDistOverTotMom /= (lV0TotalMomentum+1e-10); //avoid division by zero, to be sure
690
691//------------------------------------------------
692// Fill Tree!
693//------------------------------------------------
694
695// The conditionals are meant to decrease excessive
696// memory usage!
697
698//First Selection: Reject OnFly
699 if( (lOnFlyStatus == 0 && fkUseOnTheFly == kFALSE) || (lOnFlyStatus != 0 && fkUseOnTheFly == kTRUE ) ){
700 //Second Selection: rough 20-sigma band, parametric.
701 //K0Short: Enough to parametrize peak broadening with linear function.
702 Double_t lUpperLimitK0Short = (5.63707e-01) + (1.14979e-02)*fTreeVariablePt;
703 Double_t lLowerLimitK0Short = (4.30006e-01) - (1.10029e-02)*fTreeVariablePt;
704 //Lambda: Linear (for higher pt) plus exponential (for low-pt broadening)
705 //[0]+[1]*x+[2]*TMath::Exp(-[3]*x)
706 Double_t lUpperLimitLambda = (1.13688e+00) + (5.27838e-03)*fTreeVariablePt + (8.42220e-02)*TMath::Exp(-(3.80595e+00)*fTreeVariablePt);
707 Double_t lLowerLimitLambda = (1.09501e+00) - (5.23272e-03)*fTreeVariablePt - (7.52690e-02)*TMath::Exp(-(3.46339e+00)*fTreeVariablePt);
708 //Do Selection
709 if( (fTreeVariableInvMassLambda < lUpperLimitLambda && fTreeVariableInvMassLambda > lLowerLimitLambda ) ||
710 (fTreeVariableInvMassAntiLambda < lUpperLimitLambda && fTreeVariableInvMassAntiLambda > lLowerLimitLambda ) ||
711 (fTreeVariableInvMassK0s < lUpperLimitK0Short && fTreeVariableInvMassK0s > lLowerLimitK0Short ) ){
712 //Pre-selection in case this is AA...
713 if( fkIsNuclear == kFALSE ) fTree->Fill();
714 if( fkIsNuclear == kTRUE){
715 //If this is a nuclear collision___________________
716 // ... pre-filter with TPC, daughter eta selection
717 if( (fTreeVariableInvMassLambda < lUpperLimitLambda && fTreeVariableInvMassLambda > lLowerLimitLambda
718 && TMath::Abs(fTreeVariableNSigmasPosProton) < 6.0 && TMath::Abs(fTreeVariableNSigmasNegPion) < 6.0 ) ||
719 (fTreeVariableInvMassAntiLambda < lUpperLimitLambda && fTreeVariableInvMassAntiLambda > lLowerLimitLambda
720 && TMath::Abs(fTreeVariableNSigmasNegProton) < 6.0 && TMath::Abs(fTreeVariableNSigmasPosPion) < 6.0 ) ||
721 (fTreeVariableInvMassK0s < lUpperLimitK0Short && fTreeVariableInvMassK0s > lLowerLimitK0Short
722 && TMath::Abs(fTreeVariableNSigmasNegPion) < 6.0 && TMath::Abs(fTreeVariableNSigmasPosPion) < 6.0 ) ){
723 //insane test
724 if ( TMath::Abs(fTreeVariableNegEta)<0.8 && TMath::Abs(fTreeVariablePosEta)<0.8 ) fTree->Fill();
725 }
726 }//end nuclear_____________________________________
727 }
728 }
729
730//------------------------------------------------
731// Fill tree over.
732//------------------------------------------------
733
734 }// This is the end of the V0 loop
735
736 // Post output data.
737 PostData(1, fListHistV0);
738 PostData(2, fTree);
739
740
741}
742
743//________________________________________________________________________
744void AliAnalysisTaskExtractV0AOD::Terminate(Option_t *)
745{
746 // Draw result to the screen
747 // Called once at the end of the query
748 // This will draw the V0 candidate multiplicity, whose
749 // number of entries corresponds to the number of triggered events.
750 TList *cRetrievedList = 0x0;
751 cRetrievedList = (TList*)GetOutputData(1);
752 if(!cRetrievedList){
753 Printf("ERROR - AliAnalysisTaskExtractV0AOD : ouput data container list not available\n");
754 return;
755 }
756 fHistV0MultiplicityForTrigEvt = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistV0MultiplicityForTrigEvt") );
757 if (!fHistV0MultiplicityForTrigEvt) {
758 Printf("ERROR - AliAnalysisTaskExtractV0AOD : fHistV0MultiplicityForTrigEvt not available");
759 return;
760 }
761 TCanvas *canCheck = new TCanvas("AliAnalysisTaskExtractV0AOD","V0 Multiplicity",10,10,510,510);
762 canCheck->cd(1)->SetLogy();
763 fHistV0MultiplicityForTrigEvt->SetMarkerStyle(22);
764 fHistV0MultiplicityForTrigEvt->DrawCopy("E");
765}
766
767//________________________________________________________________________
768void AliAnalysisTaskExtractV0AOD::CheckChargeV0(AliESDv0 *v0)
769{
770 // This function checks charge of negative and positive daughter tracks.
771 // If incorrectly defined (onfly vertexer), swaps out.
772 if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ){
773 //V0 daughter track swapping is required! Note: everything is swapped here... P->N, N->P
774 Long_t lCorrectNidx = v0->GetPindex();
775 Long_t lCorrectPidx = v0->GetNindex();
776 Double32_t lCorrectNmom[3];
777 Double32_t lCorrectPmom[3];
778 v0->GetPPxPyPz( lCorrectNmom[0], lCorrectNmom[1], lCorrectNmom[2] );
779 v0->GetNPxPyPz( lCorrectPmom[0], lCorrectPmom[1], lCorrectPmom[2] );
780
781 AliExternalTrackParam lCorrectParamN(
782 v0->GetParamP()->GetX() ,
783 v0->GetParamP()->GetAlpha() ,
784 v0->GetParamP()->GetParameter() ,
785 v0->GetParamP()->GetCovariance()
786 );
787 AliExternalTrackParam lCorrectParamP(
788 v0->GetParamN()->GetX() ,
789 v0->GetParamN()->GetAlpha() ,
790 v0->GetParamN()->GetParameter() ,
791 v0->GetParamN()->GetCovariance()
792 );
793 lCorrectParamN.SetMostProbablePt( v0->GetParamP()->GetMostProbablePt() );
794 lCorrectParamP.SetMostProbablePt( v0->GetParamN()->GetMostProbablePt() );
795
796 //Get Variables___________________________________________________
797 Double_t lDcaV0Daughters = v0 -> GetDcaV0Daughters();
798 Double_t lCosPALocal = v0 -> GetV0CosineOfPointingAngle();
799 Bool_t lOnFlyStatusLocal = v0 -> GetOnFlyStatus();
800
801 //Create Replacement Object_______________________________________
802 AliESDv0 *v0correct = new AliESDv0(lCorrectParamN,lCorrectNidx,lCorrectParamP,lCorrectPidx);
803 v0correct->SetDcaV0Daughters ( lDcaV0Daughters );
804 v0correct->SetV0CosineOfPointingAngle ( lCosPALocal );
805 v0correct->ChangeMassHypothesis ( kK0Short );
806 v0correct->SetOnFlyStatus ( lOnFlyStatusLocal );
807
808 //Reverse Cluster info..._________________________________________
809 v0correct->SetClusters( v0->GetClusters( 1 ), v0->GetClusters ( 0 ) );
810
811 *v0 = *v0correct;
812 //Proper cleanup..._______________________________________________
813 v0correct->Delete();
814 v0correct = 0x0;
815
816 //Just another cross-check and output_____________________________
817 if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ) {
818 AliWarning("Found Swapped Charges, tried to correct but something FAILED!");
819 }else{
820 //AliWarning("Found Swapped Charges and fixed.");
821 }
822 //________________________________________________________________
823 }else{
824 //Don't touch it! ---
825 //Printf("Ah, nice. Charges are already ordered...");
826 }
827 return;
828}