First commit of V0 QA Task. A long way from perfect but probably ok to catch early...
[u/mrichter/AliRoot.git] / PWGLF / QATasks / AliAnalysisTaskQAV0.cxx
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 // QA Task designed to investigate V0 characteristics in any data
19 // 
20 // Stripped down and adapted from AliAnalysisTaskExtractV0: 
21 // --- smaller, more convenient output for debugging
22 // --- optional TH3 list of histos to enable a light-weight analysis
23 // --- This code is under development, so...
24 //
25 //  Please Report Any Bugs! 
26 //
27 //   --- David Dobrigkeit Chinellato
28 //        (david.chinellato@gmail.com)
29 //
30 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
31
32 class TTree;
33 class TParticle;
34 class TVector3;
35
36 //class AliMCEventHandler;
37 //class AliMCEvent;
38 //class AliStack;
39
40 class AliESDVertex;
41 class AliAODVertex;
42 class AliESDv0;
43 class AliAODv0;
44
45 #include <Riostream.h>
46 #include "TList.h"
47 #include "TH1.h"
48 #include "TH2.h"
49 #include "TH3.h"
50 #include "THnSparse.h"
51 #include "TVector3.h"
52 #include "TCanvas.h"
53 #include "TMath.h"
54 #include "TLegend.h"
55 #include "AliLog.h"
56 #include "AliCentrality.h"
57 #include "AliESDEvent.h"
58 #include "AliAODEvent.h"
59 #include "AliV0vertexer.h"
60 #include "AliCascadeVertexer.h"
61 #include "AliESDpid.h"
62 #include "AliESDtrack.h"
63 #include "AliESDtrackCuts.h"
64 #include "AliInputEventHandler.h"
65 #include "AliAnalysisManager.h"
66 #include "AliMCEventHandler.h"
67
68 #include "AliCFContainer.h"
69 #include "AliMultiplicity.h"
70
71 #include "AliESDcascade.h"
72 #include "AliAODcascade.h"
73 #include "AliESDUtils.h"
74 #include "AliESDHeader.h"
75
76 #include "AliAnalysisUtils.h"
77 #include "AliAnalysisTaskQAV0.h"
78
79 //debugging purposes
80 #include "TObjectTable.h"
81
82 ClassImp(AliAnalysisTaskQAV0)
83
84 AliAnalysisTaskQAV0::AliAnalysisTaskQAV0() 
85   : AliAnalysisTaskSE(), 
86   //Output lists
87   fOutput(0), 
88
89   //Histos
90   fHistEvent(0),
91   fHistTopDCANegToPV(0),
92   fHistTopDCAPosToPV(0),
93   fHistTopDCAV0Daughters(0),
94   fHistTopCosinePA(0),
95   fHistTopV0Radius(0),
96   fHistSelectedTopDCANegToPV(0),
97   fHistSelectedTopDCAPosToPV(0),
98   fHistSelectedTopDCAV0Daughters(0),
99   fHistSelectedTopCosinePA(0),
100   fHistSelectedTopV0Radius(0),
101
102   f2dHistInvMassK0Short(0),
103   f2dHistInvMassLambda(0),
104   f2dHistInvMassAntiLambda(0),
105
106   f2dHistInvMassWithdEdxK0Short(0),
107   f2dHistInvMassWithdEdxLambda(0),
108   f2dHistInvMassWithdEdxAntiLambda(0),
109
110   f2dHistResponseNegativeAsPion(0),
111   f2dHistResponseNegativeAsProton(0),
112   f2dHistResponsePositiveAsPion(0),
113   f2dHistResponsePositiveAsProton(0),
114
115   f2dHistdEdxSignalPionFromLambda(0),
116   f2dHistdEdxSignalProtonFromLambda(0),
117   f2dHistResponsePionFromLambda(0),
118   f2dHistResponseProtonFromLambda(0),
119
120
121   //Task Control / Utils
122   fPIDResponse(0),
123   fkRunV0Vertexer ( kFALSE )
124 {
125   // Dummy Constructor
126   for(Int_t iV0selIdx   = 0; iV0selIdx   < 7; iV0selIdx++   ) { fV0Sels          [iV0selIdx   ] = -1.; }
127 }
128
129 AliAnalysisTaskQAV0::AliAnalysisTaskQAV0(const char *name) 
130   : AliAnalysisTaskSE(name), 
131   //Output lists
132   fOutput(0), 
133
134   //Histos
135   fHistEvent(0),
136   fHistTopDCANegToPV(0),
137   fHistTopDCAPosToPV(0),
138   fHistTopDCAV0Daughters(0),
139   fHistTopCosinePA(0),
140   fHistTopV0Radius(0),
141   fHistSelectedTopDCANegToPV(0),
142   fHistSelectedTopDCAPosToPV(0),
143   fHistSelectedTopDCAV0Daughters(0),
144   fHistSelectedTopCosinePA(0),
145   fHistSelectedTopV0Radius(0),
146
147   f2dHistInvMassK0Short(0),
148   f2dHistInvMassLambda(0),
149   f2dHistInvMassAntiLambda(0),
150
151   f2dHistInvMassWithdEdxK0Short(0),
152   f2dHistInvMassWithdEdxLambda(0),
153   f2dHistInvMassWithdEdxAntiLambda(0),
154
155   f2dHistResponseNegativeAsPion(0),
156   f2dHistResponseNegativeAsProton(0),
157   f2dHistResponsePositiveAsPion(0),
158   f2dHistResponsePositiveAsProton(0),
159
160   f2dHistdEdxSignalPionFromLambda(0),
161   f2dHistdEdxSignalProtonFromLambda(0),
162   f2dHistResponsePionFromLambda(0),
163   f2dHistResponseProtonFromLambda(0),
164
165   //Task Control / Utils
166   fPIDResponse(0),
167   fkRunV0Vertexer ( kFALSE )
168 {
169   // Constructor
170   //VERTEXER CUTS
171   // REALLY LOOSE? Be careful when attempting to run over PbPb if fkRunV0Vertexer is set! 
172   fV0VertexerSels[0] =  33.  ;  // max allowed chi2
173   fV0VertexerSels[1] =   0.02;  // min allowed impact parameter for the 1st daughter (LHC09a4 : 0.05)
174   fV0VertexerSels[2] =   0.02;  // min allowed impact parameter for the 2nd daughter (LHC09a4 : 0.05)
175   fV0VertexerSels[3] =   2.0 ;  // max allowed DCA between the daughter tracks       (LHC09a4 : 0.5)
176   fV0VertexerSels[4] =   0.95;  // min allowed cosine of V0's pointing angle         (LHC09a4 : 0.99)
177   fV0VertexerSels[5] =   0.5 ;  // min radius of the fiducial volume                 (LHC09a4 : 0.2)
178   fV0VertexerSels[6] = 200.  ;  // max radius of the fiducial volume                 (LHC09a4 : 100.0)
179
180   //SELECTION CUTS
181   // REALLY LOOSE? Be careful when attempting to run over PbPb if fkRunV0Vertexer is set! 
182   fV0Sels[0] =  33.  ;  // max allowed chi2
183   fV0Sels[1] =   0.02;  // min allowed impact parameter for the 1st daughter (LHC09a4 : 0.05)
184   fV0Sels[2] =   0.02;  // min allowed impact parameter for the 2nd daughter (LHC09a4 : 0.05)
185   fV0Sels[3] =   2.0 ;  // max allowed DCA between the daughter tracks       (LHC09a4 : 0.5)
186   fV0Sels[4] =   0.95;  // min allowed cosine of V0's pointing angle         (LHC09a4 : 0.99)
187   fV0Sels[5] =   0.5 ;  // min radius of the fiducial volume                 (LHC09a4 : 0.2)
188   fV0Sels[6] = 200.  ;  // max radius of the fiducial volume                 (LHC09a4 : 100.0)
189
190   // Output slot #0 writes into a TList container (Lambda Histos and fTree)
191   DefineOutput(1, TList::Class());
192 }
193
194
195 AliAnalysisTaskQAV0::~AliAnalysisTaskQAV0()
196 {
197 //------------------------------------------------
198 // DESTRUCTOR
199 //------------------------------------------------
200
201    if (fOutput){
202       delete fOutput;
203       fOutput = 0x0;
204    }
205 }
206
207 //________________________________________________________________________
208 void AliAnalysisTaskQAV0::UserCreateOutputObjects()
209 {
210   //Define Output Lists
211   fOutput = new TList();
212   fOutput->SetOwner(); 
213
214   //Histogram Output: Event-by-Event
215   fHistEvent = new TH1D( "fHistEvent", ";Evt. Sel. Step;Count",4,0,4); 
216   fHistEvent->GetXaxis()->SetBinLabel(1, "Processed");
217   fHistEvent->GetXaxis()->SetBinLabel(2, "Phys-Sel");  
218   fHistEvent->GetXaxis()->SetBinLabel(3, "Has Vtx");  
219   fHistEvent->GetXaxis()->SetBinLabel(4, "Vtx |z|<10cm");  
220   fOutput->Add(fHistEvent); 
221
222   //Topological Selection Histograms, 1D 
223   fHistTopDCANegToPV        = new TH1D( "fHistTopDCANegToPV",";DCA Neg. Daughter to PV (cm);Counts",200,0,1); 
224   fHistTopDCAPosToPV        = new TH1D( "fHistTopDCAPosToPV",";DCA Pos. Daughter to PV (cm);Counts",200,0,1); 
225   fHistTopDCAV0Daughters    = new TH1D( "fHistTopDCAV0Daughters",";DCA V0 Daughters (#sigma);Counts",200,0,2); 
226   fHistTopCosinePA          = new TH1D( "fHistTopCosinePA",";Cosine of PA;Counts",200,-1,1); 
227   fHistTopV0Radius          = new TH1D( "fHistTopV0Radius",";Decay Radius (cm);Counts",200,0.,10);
228
229   fOutput->Add( fHistTopDCANegToPV     );
230   fOutput->Add( fHistTopDCAPosToPV     );
231   fOutput->Add( fHistTopDCAV0Daughters );
232   fOutput->Add( fHistTopCosinePA       );
233   fOutput->Add( fHistTopV0Radius       );
234
235   //Zoomed In 
236   fHistSelectedTopDCANegToPV        = new TH1D( "fHistSelectedTopDCANegToPV",";DCA Neg. Daughter to PV (cm);Counts",200,fV0Sels[1],1); 
237   fHistSelectedTopDCAPosToPV        = new TH1D( "fHistSelectedTopDCAPosToPV",";DCA Pos. Daughter to PV (cm);Counts",200,fV0Sels[2],1); 
238   fHistSelectedTopDCAV0Daughters    = new TH1D( "fHistSelectedTopDCAV0Daughters",";DCA V0 Daughters (#sigma);Counts",200,0,fV0Sels[3]); 
239   fHistSelectedTopCosinePA          = new TH1D( "fHistSelectedTopCosinePA",";Cosine of PA;Counts",200,fV0Sels[4],1); 
240   fHistSelectedTopV0Radius          = new TH1D( "fHistSelectedTopV0Radius",";Decay Radius (cm);Counts",200,fV0Sels[5],10);
241
242   fOutput->Add( fHistSelectedTopDCANegToPV     );
243   fOutput->Add( fHistSelectedTopDCAPosToPV     );
244   fOutput->Add( fHistSelectedTopDCAV0Daughters );
245   fOutput->Add( fHistSelectedTopCosinePA       );
246   fOutput->Add( fHistSelectedTopV0Radius       );
247
248   //Invariant Mass Plots
249   f2dHistInvMassK0Short     = new TH2D( "f2dHistInvMassK0Short"   , ";p_{T};M(#pi^{+},#pi^{-})"   ,250,0,25,500,0.25,0.75);
250   f2dHistInvMassLambda      = new TH2D( "f2dHistInvMassLambda"    , ";p_{T};M(p,#pi^{-})"         ,250,0,25,300,1.07,1.115+0.255);
251   f2dHistInvMassAntiLambda  = new TH2D( "f2dHistInvMassAntiLambda", ";p_{T};M(#pi^{+},#bar{p})"   ,250,0,25,300,1.07,1.115+0.255);
252
253   fOutput->Add( f2dHistInvMassK0Short     );
254   fOutput->Add( f2dHistInvMassLambda      );
255   fOutput->Add( f2dHistInvMassAntiLambda  );
256
257   //Invariant Mass Plots, with dEdx
258   f2dHistInvMassWithdEdxK0Short     = new TH2D( "f2dHistInvMassWithdEdxK0Short"   , ";p_{T};M(#pi^{+},#pi^{-})"   ,250,0,25,500,0.25,0.75);
259   f2dHistInvMassWithdEdxLambda      = new TH2D( "f2dHistInvMassWithdEdxLambda"    , ";p_{T};M(p,#pi^{-})"         ,250,0,25,300,1.07,1.115+0.255);
260   f2dHistInvMassWithdEdxAntiLambda  = new TH2D( "f2dHistInvMassWithdEdxAntiLambda", ";p_{T};M(#pi^{+},#bar{p})"   ,250,0,25,300,1.07,1.115+0.255);
261
262   fOutput->Add( f2dHistInvMassWithdEdxK0Short     );
263   fOutput->Add( f2dHistInvMassWithdEdxLambda      );
264   fOutput->Add( f2dHistInvMassWithdEdxAntiLambda  );
265
266   //dE/dx QA for main analysis and for calibration check 
267   f2dHistResponseNegativeAsPion    = new TH2D( "f2dHistResponseNegativeAsPion", ";p_{T}^{V0};N#sigma",500,0,5,400,-20,20);
268   f2dHistResponseNegativeAsProton  = new TH2D( "f2dHistResponseNegativeAsProton", ";p_{T}^{V0};N#sigma",500,0,5,400,-20,20);
269   f2dHistResponsePositiveAsPion    = new TH2D( "f2dHistResponsePositiveAsPion", ";p_{T}^{V0};N#sigma",500,0,5,400,-20,20);
270   f2dHistResponsePositiveAsProton  = new TH2D( "f2dHistResponsePositiveAsProton", ";p_{T}^{V0};N#sigma",500,0,5,400,-20,20);
271
272   //Clean Signal Check from Lambdas: stricter cuts, raw signal check 
273   f2dHistdEdxSignalPionFromLambda    = new TH2D( "f2dHistdEdxSignalPionFromLambda", ";p_{T}^{V0};TPC Signal",500,0,5,8000,0,800);
274   f2dHistdEdxSignalProtonFromLambda  = new TH2D( "f2dHistdEdxSignalProtonFromLambda", ";p_{T}^{V0};TPC Signal",500,0,5,8000,0,800);
275   f2dHistResponsePionFromLambda      = new TH2D( "f2dHistResponsePionFromLambda", ";p_{T}^{V0};N#sigma",500,0,5,400,-20,20);
276   f2dHistResponseProtonFromLambda    = new TH2D( "f2dHistResponseProtonFromLambda", ";p_{T}^{V0};N#sigma",500,0,5,400,-20,20);
277
278
279   fOutput->Add( f2dHistResponseNegativeAsPion        );
280   fOutput->Add( f2dHistResponseNegativeAsProton      );
281   fOutput->Add( f2dHistResponsePositiveAsPion        );
282   fOutput->Add( f2dHistResponsePositiveAsProton      );
283
284   fOutput->Add( f2dHistdEdxSignalPionFromLambda        );
285   fOutput->Add( f2dHistdEdxSignalProtonFromLambda      );
286   fOutput->Add( f2dHistResponsePionFromLambda          );
287   fOutput->Add( f2dHistResponseProtonFromLambda        );
288
289 //------------------------------------------------
290 // Particle Identification Setup
291 //------------------------------------------------
292
293   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
294   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
295   fPIDResponse = inputHandler->GetPIDResponse();
296   
297    //Regular output: Histograms
298    PostData(1, fOutput);
299 }// end UserCreateOutputObjects
300
301
302 //________________________________________________________________________
303 void AliAnalysisTaskQAV0::UserExec(Option_t *) 
304 {
305    // Main loop
306    // Called for each event
307    //gObjectTable->Print();
308    AliESDEvent *lESDevent = 0x0;
309
310    //AliAODEvent *lAODevent = 0x0;
311    Int_t    nV0s                        = -1;
312
313    Double_t lTrkgPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
314    Double_t lBestPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
315    Double_t lMagneticField                 = -10.;
316         
317    // Connect to the InputEvent 
318    // After these lines, we should have an ESD/AOD event + the number of cascades in it.
319
320    // Appropriate for ESD analysis! 
321    lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
322    if (!lESDevent) {
323       AliWarning("ERROR: lESDevent not available \n");
324       return;
325    }
326   
327   //------------------------------------------------
328   // Rerun V0 vertexer, if asked for
329   // --- WARNING: Be careful when using in PbPb
330   //------------------------------------------------
331   if( fkRunV0Vertexer ){
332     lESDevent->ResetV0s();
333     AliV0vertexer lV0vtxer;
334     lV0vtxer.SetDefaultCuts(fV0VertexerSels);
335     lV0vtxer.Tracks2V0vertices(lESDevent);
336   }
337
338   fHistEvent->Fill(0.5); 
339
340 //------------------------------------------------
341 // Physics Selection
342 //------------------------------------------------
343   
344   // new method
345   UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
346   Bool_t isSelected = 0;
347   isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
348   
349   //pA triggering: CINT7
350   //if( fkSwitchINT7 ) isSelected = (maskIsSelected & AliVEvent::kINT7) == AliVEvent::kINT7;
351
352   //Standard Min-Bias Selection
353   if ( !isSelected ) {
354     PostData(1, fOutput);
355     return;
356   }
357   fHistEvent->Fill(1.5); 
358   
359 //------------------------------------------------
360 // After Trigger Selection
361 //------------------------------------------------
362
363
364         nV0s = lESDevent->GetNumberOfV0s();
365   
366 //------------------------------------------------
367 // Getting: Primary Vertex + MagField Info
368 //------------------------------------------------
369
370         const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();
371         // get the vtx stored in ESD found with tracks
372         lPrimaryTrackingESDVtx->GetXYZ( lTrkgPrimaryVtxPos );
373         
374         const AliESDVertex *lPrimaryBestESDVtx = lESDevent->GetPrimaryVertex(); 
375   // get the best primary vertex available for the event
376         // As done in AliCascadeVertexer, we keep the one which is the best one available.
377         // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
378         // This one will be used for next calculations (DCA essentially)
379    lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
380
381    Double_t tPrimaryVtxPosition[3];
382    const AliVVertex *primaryVtx = lESDevent->GetPrimaryVertex();
383    tPrimaryVtxPosition[0] = primaryVtx->GetX();
384    tPrimaryVtxPosition[1] = primaryVtx->GetY();
385    tPrimaryVtxPosition[2] = primaryVtx->GetZ();
386
387   //------------------------------------------------
388   // Primary Vertex Requirements Section:
389   //  ---> pp and PbPb: Only requires |z|<10cm
390   //  ---> pPb: all requirements checked at this stage
391   //------------------------------------------------
392   
393   //Roberto's PV selection criteria, implemented 17th April 2013
394   
395   // vertex selection 
396   Bool_t fHasVertex = kFALSE;
397   const AliESDVertex *vertex = lESDevent->GetPrimaryVertexTracks();
398   if (vertex->GetNContributors() < 1) {
399     vertex = lESDevent->GetPrimaryVertexSPD();
400     if (vertex->GetNContributors() < 1) fHasVertex = kFALSE;
401     else fHasVertex = kTRUE;
402     TString vtxTyp = vertex->GetTitle();
403     Double_t cov[6]={0};
404     vertex->GetCovarianceMatrix(cov);
405     Double_t zRes = TMath::Sqrt(cov[5]);
406     if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) fHasVertex = kFALSE;
407   }
408   else fHasVertex = kTRUE;
409   
410   //Is First event in chunk rejection: Still present!
411   if(fHasVertex == kFALSE) {
412     AliWarning("Pb / | PV does not satisfy selection criteria!");
413     PostData(1, fOutput);
414     return;
415   }
416   
417   fHistEvent->Fill(2.5); 
418
419   //17 April Fix: Always do primary vertex Z selection, after pA vertex selection from Roberto
420   if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0) {
421     AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
422     PostData(1, fOutput);
423     return;
424   }
425
426   fHistEvent->Fill(3.5); 
427   
428   lMagneticField = lESDevent->GetMagneticField( );
429
430 //------------------------------------------------
431 // Only look at events with well-established PV
432 //------------------------------------------------
433         
434    //const AliESDVertex *lPrimaryTrackingESDVtxCheck = lESDevent->GetPrimaryVertexTracks();
435    const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
436    //if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtxCheck->GetStatus() ){
437    //   AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
438    //   PostData(1, fOutput);
439    //   return;
440    //}
441   
442
443
444 //------------------------------------------------
445 // MAIN LAMBDA LOOP STARTS HERE
446 //------------------------------------------------
447
448 //Variable definition
449   Int_t    lOnFlyStatus = 0;// nv0sOn = 0, nv0sOff = 0;
450   Double_t lChi2V0 = 0;
451   Double_t lDcaV0Daughters = 0, lDcaV0ToPrimVertex = 0;
452   Double_t lDcaPosToPrimVertex = 0, lDcaNegToPrimVertex = 0;
453   Double_t lV0CosineOfPointingAngle = 0;
454   Double_t lV0Radius = 0, lPt = 0;
455   Double_t lRapK0Short = 0, lRapLambda = 0;
456   Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
457   Double_t lAlphaV0 = 0, lPtArmV0 = 0;
458   Double_t lNegEta = -100, lPosEta = -100;
459   Double_t fMinV0Pt = 0; 
460   Double_t fMaxV0Pt = 100; 
461
462    Int_t nv0s = 0;
463    nv0s = lESDevent->GetNumberOfV0s();
464
465    //for (Int_t iV0 = 0; iV0 < nv0s; iV0++) 
466    for (Int_t iV0 = 0; iV0 < nv0s; iV0++) //extra-crazy test
467    {// This is the begining of the V0 loop
468       AliESDv0 *v0 = ((AliESDEvent*)lESDevent)->GetV0(iV0);
469       if (!v0) continue;
470
471       //Only use Offline Candidates for QA       
472       lOnFlyStatus = v0->GetOnFlyStatus();
473       if( lOnFlyStatus == kTRUE ) continue; 
474
475       Double_t tDecayVertexV0[3]; v0->GetXYZ(tDecayVertexV0[0],tDecayVertexV0[1],tDecayVertexV0[2]); 
476
477       Double_t tV0mom[3];
478       v0->GetPxPyPz( tV0mom[0],tV0mom[1],tV0mom[2] ); 
479       Double_t lV0TotalMomentum = TMath::Sqrt(
480       tV0mom[0]*tV0mom[0]+tV0mom[1]*tV0mom[1]+tV0mom[2]*tV0mom[2] );
481
482       lV0Radius = TMath::Sqrt(tDecayVertexV0[0]*tDecayVertexV0[0]+tDecayVertexV0[1]*tDecayVertexV0[1]);
483
484       lPt = v0->Pt();
485       lRapK0Short = v0->RapK0Short();
486       lRapLambda  = v0->RapLambda();
487       if ((lPt<fMinV0Pt)||(fMaxV0Pt<lPt)) continue;
488
489       UInt_t lKeyPos = (UInt_t)TMath::Abs(v0->GetPindex());
490       UInt_t lKeyNeg = (UInt_t)TMath::Abs(v0->GetNindex());
491
492       Double_t lMomPos[3]; v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
493       Double_t lMomNeg[3]; v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
494
495       AliESDtrack *pTrack=((AliESDEvent*)lESDevent)->GetTrack(lKeyPos);
496       AliESDtrack *nTrack=((AliESDEvent*)lESDevent)->GetTrack(lKeyNeg);
497       if (!pTrack || !nTrack) {
498          Printf("ERROR: Could not retreive one of the daughter track");
499          continue;
500       }
501
502       //Daughter Eta for Eta selection, afterwards
503       lNegEta = nTrack->Eta();
504       lPosEta = pTrack->Eta();
505
506       // Filter like-sign V0 (next: add counter and distribution)
507       if ( pTrack->GetSign() == nTrack->GetSign()){
508          continue;
509       } 
510
511       //________________________________________________________________________
512       // Track quality cuts 
513       Float_t lPosTrackCrossedRows = pTrack->GetTPCClusterInfo(2,1);
514       Float_t lNegTrackCrossedRows = nTrack->GetTPCClusterInfo(2,1);
515       Int_t lLeastNbrCrossedRows = (Int_t) lPosTrackCrossedRows;
516       if( lNegTrackCrossedRows < lLeastNbrCrossedRows )
517          lLeastNbrCrossedRows = (Int_t) lNegTrackCrossedRows;
518
519       // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)
520       if( !(pTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
521       if( !(nTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
522      
523       //Get status flags
524       //lPosTrackStatus = pTrack->GetStatus();
525       //lNegTrackStatus = nTrack->GetStatus();
526      
527       if ( ( ( pTrack->GetTPCClusterInfo(2,1) ) < 70 ) || ( ( nTrack->GetTPCClusterInfo(2,1) ) < 70 ) ) continue;
528         
529       //GetKinkIndex condition
530       if( pTrack->GetKinkIndex(0)>0 || nTrack->GetKinkIndex(0)>0 ) continue;
531
532       //Findable clusters > 0 condition
533       if( pTrack->GetTPCNclsF()<=0 || nTrack->GetTPCNclsF()<=0 ) continue;
534
535       //Compute ratio Crossed Rows / Findable clusters
536       //Note: above test avoids division by zero! 
537       Float_t lPosTrackCrossedRowsOverFindable = lPosTrackCrossedRows / ((double)(pTrack->GetTPCNclsF())); 
538       Float_t lNegTrackCrossedRowsOverFindable = lNegTrackCrossedRows / ((double)(nTrack->GetTPCNclsF())); 
539
540       Float_t lLeastRatioCrossedRowsOverFindable = lPosTrackCrossedRowsOverFindable;
541       if( lNegTrackCrossedRowsOverFindable < lLeastRatioCrossedRowsOverFindable )
542          lLeastRatioCrossedRowsOverFindable = lNegTrackCrossedRowsOverFindable;
543
544       //Lowest Cut Level for Ratio Crossed Rows / Findable = 0.8, set here
545       if ( lLeastRatioCrossedRowsOverFindable < 0.8 ) continue;
546
547       //End track Quality Cuts
548       //________________________________________________________________________
549
550       lDcaPosToPrimVertex = TMath::Abs(pTrack->GetD(tPrimaryVtxPosition[0],
551                                                         tPrimaryVtxPosition[1],
552                                                         lMagneticField) );
553
554       lDcaNegToPrimVertex = TMath::Abs(nTrack->GetD(tPrimaryVtxPosition[0],
555                                                         tPrimaryVtxPosition[1],
556                                                         lMagneticField) );
557
558
559       lChi2V0 = v0->GetChi2V0();
560       lDcaV0Daughters = v0->GetDcaV0Daughters();
561       lDcaV0ToPrimVertex = v0->GetD(tPrimaryVtxPosition[0],tPrimaryVtxPosition[1],tPrimaryVtxPosition[2]);
562       lV0CosineOfPointingAngle = v0->GetV0CosineOfPointingAngle(tPrimaryVtxPosition[0],tPrimaryVtxPosition[1],tPrimaryVtxPosition[2]);
563       //fTreeVariableV0CosineOfPointingAngle=lV0CosineOfPointingAngle;
564
565       // Getting invariant mass infos directly from ESD
566       v0->ChangeMassHypothesis(310);
567       lInvMassK0s = v0->GetEffMass();
568       v0->ChangeMassHypothesis(3122);
569       lInvMassLambda = v0->GetEffMass();
570       v0->ChangeMassHypothesis(-3122);
571       lInvMassAntiLambda = v0->GetEffMass();
572       lAlphaV0 = v0->AlphaV0();
573       lPtArmV0 = v0->PtArmV0();
574       
575       //Official means of acquiring N-sigmas 
576       Float_t lNSigmasPosProton = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kProton );
577       Float_t lNSigmasPosPion   = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kPion );
578       Float_t lNSigmasNegProton = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kProton );
579       Float_t lNSigmasNegPion   = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kPion );
580
581 //This requires an Invariant Mass Hypothesis afterwards
582       Float_t lDistOverTotMom = TMath::Sqrt(
583                                                 TMath::Power( tDecayVertexV0[0] - lBestPrimaryVtxPos[0] , 2) +
584                                                 TMath::Power( tDecayVertexV0[1] - lBestPrimaryVtxPos[1] , 2) +
585                                                 TMath::Power( tDecayVertexV0[2] - lBestPrimaryVtxPos[2] , 2)
586                                         );
587       lDistOverTotMom /= (lV0TotalMomentum+1e-10); //avoid division by zero, to be sure
588
589 //------------------------------------------------
590 // Fill Main Output Histograms
591 //------------------------------------------------
592
593   //Topological Variable Checks, One-Dimensional      
594   fHistTopDCANegToPV      -> Fill( lDcaNegToPrimVertex      ) ; 
595   fHistTopDCAPosToPV      -> Fill( lDcaPosToPrimVertex      ) ; 
596   fHistTopDCAV0Daughters  -> Fill( lDcaV0Daughters          ) ; 
597   fHistTopCosinePA        -> Fill( lV0CosineOfPointingAngle ) ; 
598   fHistTopV0Radius        -> Fill( lV0Radius                ) ; 
599
600
601   if( lDcaNegToPrimVertex > fV0Sels[1] && lDcaPosToPrimVertex > fV0Sels[2]      && 
602       lDcaV0Daughters     < fV0Sels[3] && lV0CosineOfPointingAngle > fV0Sels[4] &&
603       lV0Radius           > fV0Sels[5] && lV0Radius < fV0Sels [6] ){ 
604
605     //Topological Variables zoomed in at selection level (whatever that may be)
606     //May be slightly redundant if no specific extra configuration was done 
607     fHistSelectedTopDCANegToPV      -> Fill( lDcaNegToPrimVertex      ) ; 
608     fHistSelectedTopDCAPosToPV      -> Fill( lDcaPosToPrimVertex      ) ; 
609     fHistSelectedTopDCAV0Daughters  -> Fill( lDcaV0Daughters          ) ; 
610     fHistSelectedTopCosinePA        -> Fill( lV0CosineOfPointingAngle ) ; 
611     fHistSelectedTopV0Radius        -> Fill( lV0Radius                ) ; 
612
613     //Specific fV0Sel selection level, but no dEdx applied 
614     f2dHistInvMassK0Short     -> Fill ( lPt , lInvMassK0s        )   ; 
615     f2dHistInvMassLambda      -> Fill ( lPt , lInvMassLambda     )   ;
616     f2dHistInvMassAntiLambda  -> Fill ( lPt , lInvMassAntiLambda )   ;
617
618     //General dE/dx QA 
619     f2dHistResponseNegativeAsPion     -> Fill( lPt, lNSigmasNegPion      ); 
620     f2dHistResponseNegativeAsProton   -> Fill( lPt, lNSigmasNegProton    ); 
621     f2dHistResponsePositiveAsPion     -> Fill( lPt, lNSigmasPosPion      ); 
622     f2dHistResponsePositiveAsProton   -> Fill( lPt, lNSigmasPosProton    ); 
623
624     //Clean Sample From Lambdas
625     //Very strict cuts to ensure dealing with good Lambdas 
626     if ( lDcaV0Daughters < 1.0 && lV0CosineOfPointingAngle > 0.999 && TMath::Abs( lInvMassK0s - 0.497614 ) > 0.012 
627           && TMath::Abs( lInvMassAntiLambda - 1.115683) > 0.08 && TMath::Abs( lInvMassLambda - 1.115683) < 0.002 ) { 
628       
629       f2dHistdEdxSignalPionFromLambda     -> Fill( lPt, nTrack-> GetTPCsignal() );
630       f2dHistdEdxSignalProtonFromLambda   -> Fill( lPt, pTrack-> GetTPCsignal() );
631       f2dHistResponsePionFromLambda     -> Fill( lPt, lNSigmasNegPion   );
632       f2dHistResponseProtonFromLambda   -> Fill( lPt, lNSigmasPosProton );      
633     }
634
635     //Specific fV0Sel selection level, dE/dx applied
636     if ( TMath::Abs(lNSigmasPosPion)   < 5 && TMath::Abs(lNSigmasNegPion)   < 5 )    f2dHistInvMassWithdEdxK0Short       -> Fill ( lPt , lInvMassK0s        )   ;
637     if ( TMath::Abs(lNSigmasPosProton) < 5 && TMath::Abs(lNSigmasNegPion)   < 5 )    f2dHistInvMassWithdEdxLambda        -> Fill ( lPt , lInvMassLambda     )   ;
638     if ( TMath::Abs(lNSigmasPosPion)   < 5 && TMath::Abs(lNSigmasNegProton) < 5 )    f2dHistInvMassWithdEdxAntiLambda    -> Fill ( lPt , lInvMassAntiLambda )   ;
639
640     
641
642
643   }
644
645
646   //Topological Variable Checks, Two-Dimensional 
647
648 //------------------------------------------------
649 // End Filling of main Histograms
650 //------------------------------------------------
651
652   }// This is the end of the V0 loop
653
654
655   // Post output data.
656   PostData(1, fOutput);
657
658 }
659
660 //________________________________________________________________________
661 void AliAnalysisTaskQAV0::Terminate(Option_t *)
662 {
663    // Draw result to the screen
664    // Called once at the end of the query
665    // This will draw the V0 candidate multiplicity, whose 
666    // number of entries corresponds to the number of triggered events. 
667    TList *cRetrievedList = 0x0;
668    cRetrievedList = (TList*)GetOutputData(1);
669    if(!cRetrievedList){
670       Printf("ERROR - AliAnalysisTaskQAV0 : ouput data container list not available\n");
671       return;
672    }            
673    fHistEvent = dynamic_cast<TH1D*> (  cRetrievedList->FindObject("fHistEvent")  );
674    if (!fHistEvent) {
675       Printf("ERROR - AliAnalysisTaskQAV0 : fHistEvent not available");
676       return;
677    }
678    TCanvas *canCheck = new TCanvas("AliAnalysisTaskQAV0","V0 Multiplicity",10,10,510,510);
679    canCheck->cd(1)->SetLogy();
680    fHistEvent->SetMarkerStyle(22);
681    fHistEvent->DrawCopy("E");
682 }