]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelor.cxx
649303c44fc5d7c44927a874ab15a6900c075779
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSELc2V0bachelor.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, 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  * appeuear 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 /* $Id$ */
17
18 //
19 //
20 //               Base class for Lc2V0 Analysis
21 //
22 //
23 //  The Lc spectra study is done 2D histograms:
24 //   cascade_invariantMass VS cascade_pT
25 //
26 //  Cuts have been centralized in AliRDHFCutsLctoV0 class
27 //
28 //-------------------------------------------------------------------------
29 //
30 //                 Authors: A.De Caro(a,b), P. Pagano(b)
31 //  (a) Centro 'E.Fermi' - Roma
32 //  (b) INFN and University of Salerno
33 //
34 //  Contatcs: decaro@sa.infn.it
35 //            paola.pagano@sa.infn.it
36 //-------------------------------------------------------------------------
37
38 #include <TSystem.h>
39 #include <TParticle.h>
40 #include <TParticlePDG.h>
41 #include <TH1F.h>
42 #include <TH1F.h>
43 #include <TH2F.h>
44 #include <TTree.h>
45 #include "TROOT.h"
46 #include <TDatabasePDG.h>
47 #include <AliAnalysisDataSlot.h>
48 #include <AliAnalysisDataContainer.h>
49 #include "AliStack.h"
50 #include "AliMCEvent.h"
51 #include "AliAnalysisManager.h"
52 #include "AliAODMCHeader.h"
53 #include "AliAODHandler.h"
54 #include "AliLog.h"
55 #include "AliExternalTrackParam.h"
56 #include "AliAODVertex.h"
57 #include "AliAODRecoDecay.h"
58 #include "AliAODRecoDecayHF.h"
59 #include "AliAODRecoCascadeHF.h"
60 #include "AliAnalysisVertexingHF.h"
61 #include "AliESDtrack.h"
62 #include "AliAODTrack.h"
63 #include "AliAODv0.h"
64 #include "AliAODMCParticle.h"
65 #include "AliAnalysisTaskSE.h"
66 #include "AliAnalysisTaskSELc2V0bachelor.h"
67 #include "AliNormalizationCounter.h"
68 #include "AliAODPidHF.h"
69 #include "AliInputEventHandler.h"
70 #include "AliESDtrackCuts.h"
71 #include "AliNeutralTrackParam.h"
72
73 using std::cout;
74 using std::endl;
75
76 ClassImp(AliAnalysisTaskSELc2V0bachelor)
77
78 //__________________________________________________________________________
79 AliAnalysisTaskSELc2V0bachelor::AliAnalysisTaskSELc2V0bachelor() : AliAnalysisTaskSE(),
80   fUseMCInfo(kFALSE),
81   fOutput(0),
82   fOutputAll(0),
83   fOutputPIDBach(0),
84   fCEvents(0),
85   fIsK0SAnalysis(kFALSE),
86   fCounter(0),
87   fAnalCuts(0),
88   fUseOnTheFlyV0(kFALSE),
89   fIsEventSelected(kFALSE),
90   fWriteVariableTree(kFALSE),
91   fVariablesTree(0),
92   fCandidateVariables(),
93   fVtx1(0),
94   fBzkG(0),
95   fAdditionalChecks(kFALSE)
96 {
97   //
98   // Default ctor
99   //
100 }
101 //___________________________________________________________________________
102 AliAnalysisTaskSELc2V0bachelor::AliAnalysisTaskSELc2V0bachelor(const Char_t* name,
103                                                                AliRDHFCutsLctoV0* analCuts, Bool_t useOnTheFly,
104                                                                Bool_t writeVariableTree, Bool_t additionalChecks) :
105   AliAnalysisTaskSE(name),
106   fUseMCInfo(kFALSE),
107   fOutput(0),
108   fOutputAll(0),
109   fOutputPIDBach(0),
110   fCEvents(0),
111   fIsK0SAnalysis(kFALSE),
112   fCounter(0),
113   fAnalCuts(analCuts),
114   fUseOnTheFlyV0(useOnTheFly),
115   fIsEventSelected(kFALSE),
116   fWriteVariableTree(writeVariableTree),
117   fVariablesTree(0),
118   fCandidateVariables(),
119   fVtx1(0),
120   fBzkG(0),
121   fAdditionalChecks(additionalChecks)
122 {
123   //
124   // Constructor. Initialization of Inputs and Outputs
125   //
126   Info("AliAnalysisTaskSELc2V0bachelor","Calling Constructor");
127
128   DefineOutput(1,TList::Class());  //conters
129   DefineOutput(2,AliNormalizationCounter::Class());
130   DefineOutput(3,AliRDHFCutsLctoV0::Class());
131   if (!writeVariableTree) {
132     DefineOutput(4,TList::Class());  //All Entries output
133     DefineOutput(5,TList::Class());  //3sigma PID output
134   } else {
135     // Output slot #4 keeps a tree of the candidate variables after track selection
136     DefineOutput(4,TTree::Class());  //My private output
137   }
138
139 }
140
141 //___________________________________________________________________________
142 AliAnalysisTaskSELc2V0bachelor::~AliAnalysisTaskSELc2V0bachelor() {
143   //
144   // destructor
145   //
146   Info("~AliAnalysisTaskSELc2V0bachelor","Calling Destructor");
147   
148   if (fOutput) {
149     delete fOutput;
150     fOutput = 0;
151   }
152
153   if (fOutputAll) {
154     delete fOutputAll;
155     fOutputAll = 0;
156   }
157
158   if (fOutputPIDBach) {
159     delete fOutputPIDBach;
160     fOutputPIDBach = 0;
161   }
162
163   if (fCounter) {
164     delete fCounter;
165     fCounter = 0;
166   }
167
168   if (fAnalCuts) {
169     delete fAnalCuts;
170     fAnalCuts = 0;
171   }
172
173   if (fVariablesTree) {
174     delete fVariablesTree;
175     fVariablesTree = 0;
176   }
177
178 }
179 //_________________________________________________
180 void AliAnalysisTaskSELc2V0bachelor::Init() {
181   //
182   // Initialization
183   //
184
185   fIsEventSelected=kFALSE;
186
187   if (fDebug > 1) AliInfo("Init");
188
189   PostData(3,fAnalCuts);
190
191   return;
192 }
193
194 //_________________________________________________
195 void AliAnalysisTaskSELc2V0bachelor::UserExec(Option_t *)
196 {
197   // user exec
198   if (!fInputEvent) {
199     AliError("NO EVENT FOUND!");
200     return;
201   }
202
203   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
204   TClonesArray *arrayLctopKos=0;
205
206   if (!aodEvent && AODEvent() && IsStandardAOD()) {
207     // In case there is an AOD handler writing a standard AOD, use the AOD 
208     // event in memory rather than the input (ESD) event.    
209     aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
210     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
211     // have to taken from the AOD event hold by the AliAODExtension
212     AliAODHandler* aodHandler = (AliAODHandler*) 
213       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
214
215     if (aodHandler->GetExtensions()) {
216       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
217       AliAODEvent *aodFromExt = ext->GetAOD();
218       arrayLctopKos=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
219     }
220   } else {
221     arrayLctopKos=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
222   }
223
224   fCEvents->Fill(1);
225
226   if (fUseMCInfo)
227     fAnalCuts->SetTriggerClass("");
228
229   // AOD primary vertex
230   fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
231   if (!fVtx1) return;
232
233   fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent); // better to initialize before CheckEventSelection call
234
235   CheckEventSelection(aodEvent);
236
237
238   // fix for temporary bug in ESDfilter 
239   fBzkG = (Double_t)aodEvent->GetMagneticField(); 
240   if (TMath::Abs(fBzkG)<0.001) return;
241   fCEvents->Fill(2);
242
243   if (!arrayLctopKos) {
244     AliInfo("Could not find array of HF cascades, skipping the event");
245     return;
246   } else {
247     if (arrayLctopKos->GetEntriesFast()) {
248       AliInfo(Form("Found %d cascades",arrayLctopKos->GetEntriesFast()));
249     }
250   }
251   fCEvents->Fill(3);
252
253   // mc analysis 
254   TClonesArray *mcArray = 0;
255   AliAODMCHeader *mcHeader=0;
256
257   if (fUseMCInfo) {
258     // MC array need for maching
259     mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
260     if (!mcArray) {
261       AliError("Could not find Monte-Carlo in AOD");
262       return;
263     }
264     fCEvents->Fill(4); // in case of MC events
265
266     // load MC header
267     mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
268     if (!mcHeader) {
269       AliError("AliAnalysisTaskSELc2V0bachelor::UserExec: MC header branch not found!\n");
270       return;
271     }
272     fCEvents->Fill(5); // in case of MC events
273
274     Double_t zMCVertex = mcHeader->GetVtxZ();
275     if (TMath::Abs(zMCVertex) > fAnalCuts->GetMaxVtxZ()) {
276       AliDebug(2,Form("Event rejected: abs(zVtxMC)=%f > fAnalCuts->GetMaxVtxZ()=%f",zMCVertex,fAnalCuts->GetMaxVtxZ()));
277       return;
278     } else {
279       fCEvents->Fill(17); // in case of MC events
280     }
281   }
282
283   fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo); // it is very important that it stays BEFORE any other event selection
284
285   if (fVtx1->GetNContributors()>0) // this check is done in IsEventSelected
286     fCEvents->Fill(6);
287
288   if ( !fIsEventSelected ) return; // don't take into account not selected events 
289   fCEvents->Fill(7);
290
291   Int_t nSelectedAnal = 0;
292   if (fIsK0SAnalysis) {
293     MakeAnalysisForLc2prK0S(arrayLctopKos,mcArray, nSelectedAnal, fAnalCuts);
294
295     if (nSelectedAnal)
296       CheckEventSelectionWithCandidates(aodEvent);
297
298   }
299
300   fCounter->StoreCandidates(aodEvent,nSelectedAnal,kTRUE);
301   fCounter->StoreCandidates(aodEvent,nSelectedAnal,kFALSE);
302
303   PostData(1,fOutput);
304   PostData(2,fCounter);
305   if (!fWriteVariableTree) {
306     PostData(4,fOutputAll);
307     PostData(5,fOutputPIDBach);
308   } else {
309     PostData(4,fVariablesTree);
310   }
311
312   fIsEventSelected=kFALSE;
313
314   return;
315 }
316
317 //________________________________________ terminate ___________________________
318 void AliAnalysisTaskSELc2V0bachelor::Terminate(Option_t*)
319 {    
320   // The Terminate() function is the last function to be called during
321   // a query. It always runs on the client, it can be used to present
322   // the results graphically or save the results to file.
323   
324   //AliInfo("Terminate","");
325   AliAnalysisTaskSE::Terminate();
326   
327   fOutput = dynamic_cast<TList*> (GetOutputData(1));
328   if (!fOutput) {     
329     AliError("fOutput not available");
330     return;
331   }
332   
333   //fCEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
334   if (!fWriteVariableTree) {
335     fOutputAll = dynamic_cast<TList*> (GetOutputData(4));
336     if (!fOutputAll) {
337       AliError("fOutputAll not available");
338       return;
339     }
340
341     fOutputPIDBach = dynamic_cast<TList*> (GetOutputData(5));
342     if (!fOutputPIDBach) {
343       AliError("fOutputPIDBach not available");
344       return;
345     }
346   }
347
348   return;
349 }
350 //___________________________________________________________________________
351 void AliAnalysisTaskSELc2V0bachelor::UserCreateOutputObjects() { 
352   // output
353   AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
354
355   fOutput = new TList();
356   fOutput->SetOwner();
357   fOutput->SetName("chist0");
358   DefineGeneralHistograms(); // define general histograms
359   PostData(1,fOutput);
360
361   fCounter = new AliNormalizationCounter("NormalizationCounter");
362   fCounter->Init();
363   PostData(2,fCounter);
364
365   if (!fWriteVariableTree) {
366
367     fOutputAll = new TList();
368     fOutputAll->SetOwner();
369     fOutputAll->SetName("listAll");
370
371     fOutputPIDBach = new TList();
372     fOutputPIDBach->SetOwner();
373     fOutputPIDBach->SetName("listPIDBach");
374
375     DefineAnalysisHistograms(); // define analysis histograms
376   
377     PostData(4,fOutputAll);
378     PostData(5,fOutputPIDBach);
379   }
380   else {
381     DefineTreeVariables();
382     PostData(4,fVariablesTree);
383   }
384
385   return;
386 }
387
388 //-------------------------------------------------------------------------------
389 void AliAnalysisTaskSELc2V0bachelor::MakeAnalysisForLc2prK0S(TClonesArray *arrayLctopKos,
390                                                              TClonesArray *mcArray,
391                                                              Int_t &nSelectedAnal,
392                                                              AliRDHFCutsLctoV0 *cutsAnal)
393 {
394
395   // make the analysis
396
397   Int_t pdgCand = 4122;
398   Int_t pdgDgLctoV0bachelor[2]={2212,310}; // always 1st bachelor, 2nd V0
399   Int_t pdgDgV0toDaughters[2]={211,211};
400
401   // loop over cascades to search for candidates Lc->p+K0S
402   Int_t nCascades= arrayLctopKos->GetEntriesFast();
403   if (nCascades==0) {
404     AliInfo("Could not find cascades, skipping the event");
405     return;
406   }
407
408   for (Int_t iLctopK0S = 0; iLctopK0S<nCascades; iLctopK0S++) {
409
410     ((TH1F*)(fOutput->FindObject("hCandidateSelection")))->Fill(0);
411
412     // Lc candidates and K0S from Lc
413     AliAODRecoCascadeHF* lcK0Spr = dynamic_cast<AliAODRecoCascadeHF*>(arrayLctopKos->At(iLctopK0S));
414     if (!lcK0Spr) {
415       AliDebug(2,Form("Cascade %d doens't exist, skipping",iLctopK0S));
416       continue;
417     }
418
419     Bool_t unsetvtx=kFALSE;
420     if (!lcK0Spr->GetOwnPrimaryVtx()) {
421       lcK0Spr->SetOwnPrimaryVtx(fVtx1);
422       unsetvtx=kTRUE;
423     }
424
425     if (!lcK0Spr->GetSecondaryVtx()) {
426       AliInfo("No secondary vertex"); // it will be done in AliRDHFCutsLctoV0::IsSelected
427       continue;
428     }
429
430     if (lcK0Spr->GetNDaughters()!=2) {
431       AliDebug(2,Form("Cascade %d has not 2 daughters (nDaughters=%d)",iLctopK0S,lcK0Spr->GetNDaughters())); // it will be done in AliRDHFCutsLctoV0::IsSelected
432       continue;
433     }
434
435     AliAODv0 * v0part = dynamic_cast<AliAODv0*>(lcK0Spr->Getv0());
436     AliAODTrack * bachPart = dynamic_cast<AliAODTrack*>(lcK0Spr->GetBachelor());
437     if (!v0part || !bachPart) {
438       AliDebug(2,Form("Cascade %d has no V0 or no bachelor object",iLctopK0S)); // it will be done in AliRDHFCutsLctoV0::IsSelected
439       continue;
440     }
441
442     if (!v0part->GetSecondaryVtx()) {
443       AliDebug(2,Form("No secondary vertex for V0 by cascade %d",iLctopK0S)); // it will be done in AliRDHFCutsLctoV0::IsSelected
444       continue;
445     }
446
447     if (v0part->GetNDaughters()!=2) {
448       AliDebug(2,Form("current V0 has not 2 daughters (onTheFly=%d, nDaughters=%d)",v0part->GetOnFlyStatus(),v0part->GetNDaughters())); // it will be done in AliRDHFCutsLctoV0::IsSelected
449       continue;
450     }
451
452     AliAODTrack * v0Pos = dynamic_cast<AliAODTrack*>(lcK0Spr->Getv0PositiveTrack());
453     AliAODTrack * v0Neg = dynamic_cast<AliAODTrack*>(lcK0Spr->Getv0NegativeTrack());
454     if (!v0Neg || !v0Neg) {
455       AliDebug(2,Form("V0 by cascade %d has no V0positive of V0negative object",iLctopK0S)); // it will be done in AliRDHFCutsLctoV0::IsSelected
456       continue;
457     }
458
459     ((TH1F*)(fOutput->FindObject("hCandidateSelection")))->Fill(1);
460
461     if (v0Pos->Charge() == v0Neg->Charge()) continue;
462
463     ((TH1F*)(fOutput->FindObject("hCandidateSelection")))->Fill(2);
464
465     Int_t isLc = 0;
466
467     if (fUseMCInfo) {
468
469       Int_t pdgCode=-2;
470
471       // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
472       Int_t mcLabel = lcK0Spr->MatchToMC(pdgCand,pdgDgLctoV0bachelor[1],pdgDgLctoV0bachelor,pdgDgV0toDaughters,mcArray,kTRUE);
473       if (mcLabel!=-1) {
474         AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0S,nCascades));
475
476         AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
477         if (partLc) {
478           pdgCode = partLc->GetPdgCode();
479           if (pdgCode<0) AliDebug(2,Form(" Â¯Â¯Â¯Â¯Â¯Â¯Â¯Â¯Â¯Â¯Â¯Â¯Â¯Â¯Â¯Â¯Â¯Â¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
480           pdgCode = TMath::Abs(pdgCode);
481           isLc = 1;
482         }
483       } else {
484         AliDebug(2,Form("No MC candidate (cascade number %d -total cascade number = %d -)", iLctopK0S,nCascades));
485         pdgCode=-1;
486       }
487     }
488
489     FillLc2pK0Sspectrum(lcK0Spr, isLc,
490                         nSelectedAnal, cutsAnal,
491                         mcArray);
492
493     if (unsetvtx) lcK0Spr->UnsetOwnPrimaryVtx();
494
495   }
496
497   AliDebug(2, Form("Found %d Reco particles that are Lc!!", nSelectedAnal));
498
499   return;
500 }
501
502 //________________________________________________________________________
503 void AliAnalysisTaskSELc2V0bachelor::FillLc2pK0Sspectrum(AliAODRecoCascadeHF *part,
504                                                          Int_t isLc,
505                                                          Int_t &nSelectedAnal,
506                                                          AliRDHFCutsLctoV0 *cutsAnal,
507                                                          TClonesArray *mcArray)
508 {
509   //
510   // Fill histos for Lc -> K0S+proton
511   //
512
513   TString fillthis="";
514
515   AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
516   Double_t momBach  = bachelor->P();
517
518   AliAODv0 * v0part = (AliAODv0*)part->Getv0();
519   Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
520
521   Bool_t areCutsUsingPID = cutsAnal->GetIsUsePID();
522   cutsAnal->SetUsePID(kFALSE);
523   Bool_t isInCascadeWindow = (((cutsAnal->IsSelectedSingleCut(part,AliRDHFCuts::kCandidate,0))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on Lc->p+K0S invMass
524   cutsAnal->SetUsePID(areCutsUsingPID);
525
526   //if ( !( !onFlyV0 || (onFlyV0 && fUseOnTheFlyV0) ) ) return;
527   if ( onFlyV0 && !fUseOnTheFlyV0 ) return;
528
529   if (fAdditionalChecks) CheckCandidatesAtDifferentLevels(part,cutsAnal);
530
531   if ( !(cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) return;
532
533   if ( !( ( (cutsAnal->IsSelected(part,AliRDHFCuts::kTracks))&(AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) ) return;
534
535   if ( ( ( (cutsAnal->IsSelected(part,AliRDHFCuts::kAll))&(AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) ) nSelectedAnal++;
536
537   // Fill candidate variable Tree (track selection, V0 invMass selection)
538   if ( fWriteVariableTree ) {
539     Double_t invmassK0S = v0part->MassK0Short();
540     Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
541     if ( !onFlyV0 && isInCascadeWindow && part->CosV0PointingAngle()>0.99 && TMath::Abs(invmassK0S-mk0sPDG)<=0.05)
542       FillTheTree(part,cutsAnal,mcArray,isLc);
543     return;
544   }
545
546   cutsAnal->SetUsePID(kFALSE);
547   Bool_t isCandidateSelectedCuts = (((cutsAnal->IsSelected(part,AliRDHFCuts::kCandidate))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // kinematic/topological cuts
548   cutsAnal->SetUsePID(areCutsUsingPID);
549   Bool_t isBachelorID = (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // ID x bachelor
550
551   //if (((cutsAnal->IsSelected(part,AliRDHFCuts::kCandidate))&AliRDHFCutsLctoV0::kLcToK0Spr)==AliRDHFCutsLctoV0::kLcToK0Spr) {
552   if (((cutsAnal->IsSelected(part,AliRDHFCuts::kAll))&AliRDHFCutsLctoV0::kLcToK0Spr)==AliRDHFCutsLctoV0::kLcToK0Spr) {
553     if (fUseMCInfo && isLc && !fWriteVariableTree) {
554       Int_t pdgCand1 = 4122;
555       Int_t pdgDgLctoV0bachelor1[2]={2212,310};
556       Int_t pdgDgV0toDaughters1[2]={211,211};
557       Int_t mcLabel1=part->MatchToMC(pdgCand1,pdgDgLctoV0bachelor1[1],pdgDgLctoV0bachelor1,pdgDgV0toDaughters1,mcArray,kTRUE);
558       AliDebug(2,Form(" Found true MC candidate: Lc->pK0S(%d) - onTheFly=%1d",mcLabel1,onFlyV0));
559     }
560   }
561
562   Double_t nSigmaTPCpr=-999.;
563   cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,4,nSigmaTPCpr);
564   Double_t nSigmaTOFpr=-999.;
565   cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,4,nSigmaTOFpr);
566
567   Double_t nSigmaTPCpi=-999.;
568   cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,2,nSigmaTPCpi);
569   Double_t nSigmaTOFpi=-999.;
570   cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,2,nSigmaTOFpi);
571
572   Double_t nSigmaTPCka=-999.;
573   cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,3,nSigmaTPCka);
574   Double_t nSigmaTOFka=-999.;
575   cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,3,nSigmaTOFka);
576
577   //if (onFlyV0 && fUseOnTheFlyV0) {  
578   if (onFlyV0) {  
579
580     if (isCandidateSelectedCuts)
581       FillAnalysisHistograms(part,isBachelorID,"");
582
583   }
584   //else if (!onFlyV0) {
585   else {
586
587     if (isCandidateSelectedCuts) {
588
589       fillthis="histoprotonBachSigmaVspTOF";
590       ((TH2F*)(fOutput->FindObject(fillthis)))->Fill(momBach,nSigmaTOFpr);
591       fillthis="histoprotonBachSigmaVspTPC";
592       ((TH2F*)(fOutput->FindObject(fillthis)))->Fill(momBach,nSigmaTPCpr);
593
594       FillAnalysisHistograms(part,isBachelorID,"Offline");
595
596     }
597
598   }
599
600
601   if (fUseMCInfo) {
602     if (isLc==1) {
603       //if (onFlyV0 && fUseOnTheFlyV0) {
604       if (onFlyV0) {
605
606         if (isCandidateSelectedCuts)
607           FillAnalysisHistograms(part,isBachelorID,"Sgn");
608     
609       }     
610       //else if (!onFlyV0) {  
611       else {  
612
613         if (isCandidateSelectedCuts) {
614
615           fillthis="histoprotonBachSigmaVspTOFsgn";
616           ((TH2F*)(fOutput->FindObject(fillthis)))->Fill(momBach,nSigmaTOFpr);
617           fillthis="histoprotonBachSigmaVspTPCsgn";
618           ((TH2F*)(fOutput->FindObject(fillthis)))->Fill(momBach,nSigmaTPCpr);
619
620           FillAnalysisHistograms(part,isBachelorID,"OfflineSgn");
621
622         }
623
624       }
625
626     }// sgn
627     else { // bkg
628       //if (onFlyV0 && fUseOnTheFlyV0) {
629       if (onFlyV0) {
630
631         if (isCandidateSelectedCuts)
632           FillAnalysisHistograms(part,isBachelorID,"Bkg");
633
634       }
635       //else if (!onFlyV0) {
636       else {
637
638         if (isCandidateSelectedCuts) {
639
640           fillthis="histoprotonBachSigmaVspTOFbkg";
641           ((TH2F*)(fOutput->FindObject(fillthis)))->Fill(momBach,nSigmaTOFpr);
642           fillthis="histoprotonBachSigmaVspTPCbkg";
643           ((TH2F*)(fOutput->FindObject(fillthis)))->Fill(momBach,nSigmaTPCpr);
644
645           FillAnalysisHistograms(part,isBachelorID,"OfflineBkg");
646         }
647
648       }
649
650     }
651   } // if fUseMCInfo
652  
653   return;
654 }
655
656 //----------------------------------------------------
657 void AliAnalysisTaskSELc2V0bachelor::DefineK0SHistos()
658
659
660   TString nameMass=" ", nameSgn=" ", nameBkg=" ";
661
662   Double_t mLcPDG  = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
663   Double_t mK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
664
665   if (fUseOnTheFlyV0) {
666
667     // V0 invariant masses (on-the-fly)
668     nameMass="histK0SMass";
669     TH2F* spectrumK0SMass = new TH2F(nameMass.Data(),"K^{0}_{S} invariant mass VS p_{T}; M(#pi^{+}#pi^{-}) [GeV/c^{2}]; p_{T}(K^{0}_{S}) [GeV/c]; Entries",
670                                     1000,mK0SPDG-0.050,mK0SPDG+0.050,175,0.,35.);
671
672     // Lc invariant masses (x K0S on-the-fly)
673     nameMass="histLcMassByK0S";
674     TH2F* spectrumLcMassByK0S = new TH2F(nameMass.Data(),"#Lambda_{c} invariant mass (by K^{0}_{S}) vs p_{T} ; m_{inv}(p-K^{0}_{S}) [GeV/c^{2}]; p_{T}(#Lambda_{c}) [GeV/c]",
675                                          1000,mLcPDG-0.250,mLcPDG+0.250,175,0.,35.);
676
677     nameMass="histpK0Svsp";
678     TH2F* momentumDistributionK0Svsp = new TH2F(nameMass.Data(),"#Lambda_{c}: p(K^{0}_{S}) vs p(p);  p_{p}; p_{K^{0}_{S}}  ",
679                                                 175,0.,35.,175,0.,35.);
680
681     nameMass="histArmPodK0S";
682     TH2F* armenterosPodK0S = new TH2F(nameMass.Data(),"K^{0}_{S}: Armenteros-Podolanski distribution; #frac{p_{L}^{+}-p_{L}^{-}}{p_{L}^{+}+p_{L}^{-}}; p_{T}^{+} [GeV/c]",
683                                       200,-1.,1.,300,0.,0.3);
684  
685     nameMass="histDCAtoPVvspK0S";
686     TH2F *dcatoPVvspK0S = new TH2F(nameMass.Data(),"K^{0}_{S}: DCA to Primary Vertex vs K^{0}_{S} momentum ; p(K^{0}_{S}) [GeV/c]; DCA to Primary Vertex [n#sigmas]; Entries",
687                                    175,0.,35.,50,0.,5.);
688
689     nameMass="histK0ScosPAwrtPVvspK0S";
690     TH2F *cosPAwrtPVvspK0S = new TH2F(nameMass.Data(),"K^{0}_{S}: cosine of pointing angle wrt primary vertex vs K^{0}_{S} momentum ; p(K^{0}_{S}) [GeV/c]; cosine; Entries",
691                                       175,0.,35.,100,0.99,1.);
692
693     TH2F* allspectrumK0SMass = (TH2F*)spectrumK0SMass->Clone(); 
694     TH2F* allspectrumLcMassByK0S    = (TH2F*)spectrumLcMassByK0S->Clone(); 
695     TH2F* allmomentumDistributionK0Svsp= (TH2F*)momentumDistributionK0Svsp->Clone(); 
696     TH2F* alldcatoPVvspK0S=(TH2F*)dcatoPVvspK0S->Clone(); 
697     TH2F* allcosV0PAwrtPVvspK0S=(TH2F*)cosPAwrtPVvspK0S->Clone(); 
698
699     TH2F* pidBachspectrumK0SMass = (TH2F*)spectrumK0SMass->Clone(); 
700     TH2F* pidBachspectrumLcMassByK0S    = (TH2F*)spectrumLcMassByK0S->Clone(); 
701     TH2F* pidBachmomentumDistributionK0Svsp= (TH2F*)momentumDistributionK0Svsp->Clone(); 
702     TH2F* pidBachdcatoPVvspK0S=(TH2F*)dcatoPVvspK0S->Clone(); 
703     TH2F* pidBachcosV0PAwrtPVvspK0S=(TH2F*)cosPAwrtPVvspK0S->Clone(); 
704
705     TH2F* allArmenterosPodK0S = (TH2F*)armenterosPodK0S->Clone();
706     TH2F* pidBachArmenterosPodK0S = (TH2F*)armenterosPodK0S->Clone();
707
708     fOutputAll->Add(allspectrumK0SMass);
709     fOutputAll->Add(allspectrumLcMassByK0S);
710     fOutputAll->Add(allmomentumDistributionK0Svsp); 
711     fOutputAll->Add(allArmenterosPodK0S);
712     fOutputAll->Add(alldcatoPVvspK0S);
713     fOutputAll->Add(allcosV0PAwrtPVvspK0S);
714
715     fOutputPIDBach->Add(pidBachspectrumK0SMass);
716     fOutputPIDBach->Add(pidBachspectrumLcMassByK0S);
717     fOutputPIDBach->Add(pidBachmomentumDistributionK0Svsp); 
718     fOutputPIDBach->Add(pidBachArmenterosPodK0S);
719     fOutputPIDBach->Add(pidBachdcatoPVvspK0S);
720     fOutputPIDBach->Add(pidBachcosV0PAwrtPVvspK0S);
721  
722   }
723
724   // V0 invariant masses (offline)
725   nameMass="histK0SMassOffline";
726   TH2F* spectrumK0SMassOffline = new TH2F(nameMass.Data(),"K^{0}_{S} invariant mass VS p_{T}; M(#pi^{+}#pi^{-}) [GeV/c^{2}]; p_{T}(K^{0}_{S}) [GeV/c]; Entries",
727                                          1000,mK0SPDG-0.050,mK0SPDG+0.050,175,0.,35.);
728
729   // Lc invariant masses (x K0S offline)
730   nameMass="histLcMassByK0SOffline";
731   TH2F* spectrumLcMassOfflineByK0S = new TH2F(nameMass.Data(),"#Lambda_{c} invariant mass (by K^{0}_{S}) vs p_{T}; M(K^{0}_{S}p) [GeV/c^{2}]; p_{T} [GeV/c]",
732                                               1000,mLcPDG-0.250,mLcPDG+0.250,175,0.,35.);
733
734   nameMass="histpK0SvspOffline";
735   TH2F* momentumDistributionK0SvspOffline = new TH2F(nameMass.Data(),"#Lambda_{c}: p(K^{0}_{S}) vs p(p) - Offline ;  p_{p} [GeV/c]; p_{K^{0}_{S}} [GeV/c]",
736                                                      175,0.,35.,175,0.,35.);
737
738   nameMass="histArmPodK0SOffline";
739   TH2F* armenterosPodK0SOff = new TH2F(nameMass.Data(),"K^{0}_{S}  Armenteros-Podolanski distribution - Offline; #frac{p_{L}^{+}-p_{L}^{-}}{p_{L}^{+}+p_{L}^{-}}; p_{T}^{+} [GeV/c]",
740                                        200,-1.,1.,300,0.,0.3);
741
742   nameMass="histDCAtoPVvspK0SOffline";
743   TH2F *dcatoPVvspK0SOffline = new TH2F(nameMass.Data(),"K^{0}_{S}: DCA to Primary Vertex vs  K^{0}_{S} invariant mass - Offline; p(K^{0}_{S}) [GeV/c]; DCA to Primary Vertex [n#sigmas]; Entries",
744                                         175,0.,35.,50,0.,5.);
745
746   nameMass="histK0ScosPAwrtPVvspK0SOffline";
747   TH2F *cosPAwrtPVvspK0SOffline = new TH2F(nameMass.Data(),"K^{0}_{S}: cosine of pointing angle wrt primary vertex vs K^{0}_{S} momentum - Offline; p(K^{0}_{S}) [GeV/c]; cosine; Entries",
748                                            175,0.,35.,100,0.99,1.);
749
750
751
752   TH2F* allspectrumK0SMassOffline = (TH2F*)spectrumK0SMassOffline->Clone(); 
753   TH2F* allspectrumLcMassOfflineByK0S    = (TH2F*)spectrumLcMassOfflineByK0S->Clone(); 
754   TH2F* allmomentumDistributionK0SvspOffline= (TH2F*)momentumDistributionK0SvspOffline->Clone(); 
755   TH2F* alldcatoPVvspK0SOffline=(TH2F*)dcatoPVvspK0SOffline->Clone(); 
756   TH2F* allcosPAwrtPVvspK0SOffline=(TH2F*)cosPAwrtPVvspK0SOffline->Clone(); 
757   TH2F* allArmenterosPodK0SOff = (TH2F*)armenterosPodK0SOff->Clone();
758
759   TH2F* pidBachspectrumK0SMassOffline = (TH2F*)spectrumK0SMassOffline->Clone(); 
760   TH2F* pidBachspectrumLcMassOfflineByK0S    = (TH2F*)spectrumLcMassOfflineByK0S->Clone(); 
761   TH2F* pidBachmomentumDistributionK0SvspOffline= (TH2F*)momentumDistributionK0SvspOffline->Clone(); 
762   TH2F* pidBachdcatoPVvspK0SOffline=(TH2F*)dcatoPVvspK0SOffline->Clone(); 
763   TH2F* pidBachcosPAwrtPVvspK0SOffline=(TH2F*)cosPAwrtPVvspK0SOffline->Clone(); 
764   TH2F* pidBachArmenterosPodK0SOff = (TH2F*)armenterosPodK0SOff->Clone();
765
766
767   fOutputAll->Add(allspectrumK0SMassOffline);
768   fOutputAll->Add(allspectrumLcMassOfflineByK0S);
769   fOutputAll->Add(allmomentumDistributionK0SvspOffline); 
770   fOutputAll->Add(allArmenterosPodK0SOff);
771   fOutputAll->Add(alldcatoPVvspK0SOffline);
772   fOutputAll->Add(allcosPAwrtPVvspK0SOffline);
773
774   fOutputPIDBach->Add(pidBachspectrumK0SMassOffline);
775   fOutputPIDBach->Add(pidBachspectrumLcMassOfflineByK0S);
776   fOutputPIDBach->Add(pidBachmomentumDistributionK0SvspOffline); 
777   fOutputPIDBach->Add(pidBachArmenterosPodK0SOff);
778   fOutputPIDBach->Add(pidBachdcatoPVvspK0SOffline);
779   fOutputPIDBach->Add(pidBachcosPAwrtPVvspK0SOffline);
780
781   if (fUseMCInfo) {
782
783     if (fUseOnTheFlyV0) {
784
785       nameSgn="histK0SMassSgn";
786       nameBkg="histK0SMassBkg";
787       TH2F* spectrumK0SMassSgn = new TH2F(nameSgn.Data(), "K^{0}_{S} - sgn: invariant mass VS p_{T} - MC; m_{inv}(#pi^{+}#pi^{-}) [GeV/c^{2}]; p_{T}(K^{0}_{S}) [GeV/c]; Entries",
788                                           1000,mK0SPDG-0.050,mK0SPDG+0.050,175,0.,35.);
789       TH2F* spectrumK0SMassBkg = new TH2F(nameBkg.Data(), "K^{0}_{S} - bkg: invariant mass VS p_{T} - MC; m_{inv}(#pi^{+}#pi^{-}) [GeV/c^{2}]; p_{T}(K^{0}_{S}) [GeV/c]; Entries",
790                                           1000,mK0SPDG-0.050,mK0SPDG+0.050,175,0.,35.);
791
792       TH2F* allspectrumK0SMassSgn = (TH2F*)spectrumK0SMassSgn->Clone(); 
793       TH2F* allspectrumK0SMassBkg = (TH2F*) spectrumK0SMassBkg->Clone();  
794       TH2F* pidBachspectrumK0SMassSgn = (TH2F*)spectrumK0SMassSgn->Clone(); 
795       TH2F* pidBachspectrumK0SMassBkg = (TH2F*) spectrumK0SMassBkg->Clone();  
796       fOutputAll->Add(allspectrumK0SMassSgn);
797       fOutputAll->Add(allspectrumK0SMassBkg);
798       fOutputPIDBach->Add(pidBachspectrumK0SMassSgn);
799       fOutputPIDBach->Add(pidBachspectrumK0SMassBkg);
800
801
802       nameSgn="histLcMassByK0SSgn";
803       nameBkg="histLcMassByK0SBkg";
804       TH2F* spectrumLcMassByK0SSgn = new TH2F(nameSgn.Data(), "#Lambda_{c} - sgn: invariant mass (by K^{0}_{S}) vs p_{T}  - MC; m_{inv}(p-K^{0}_{S}) [GeV/c^{2}];  p_{T}",
805                                               1000,mLcPDG-0.250,mLcPDG+0.250,175,0.,35.);
806       TH2F* spectrumLcMassByK0SBkg = new TH2F(nameBkg.Data(), "#Lambda_{c} - bkg: invariant mass (by K^{0}_{S}) vs p_{T}  - MC; m_{inv}(p-K^{0}_{S}) [GeV/c^{2}]; p_{T}",
807                                               1000,mLcPDG-0.250,mLcPDG+0.250,175,0.,35.);
808
809       TH2F* allspectrumLcMassByK0SSgn = (TH2F*)spectrumLcMassByK0SSgn->Clone(); 
810       TH2F* allspectrumLcMassByK0SBkg = (TH2F*) spectrumLcMassByK0SBkg->Clone();  
811       TH2F* pidBachspectrumLcMassByK0SSgn = (TH2F*)spectrumLcMassByK0SSgn->Clone(); 
812       TH2F* pidBachspectrumLcMassByK0SBkg = (TH2F*) spectrumLcMassByK0SBkg->Clone();  
813       fOutputAll->Add(allspectrumLcMassByK0SSgn);
814       fOutputAll->Add(allspectrumLcMassByK0SBkg);
815       fOutputPIDBach->Add(pidBachspectrumLcMassByK0SSgn);
816       fOutputPIDBach->Add(pidBachspectrumLcMassByK0SBkg);
817
818       nameSgn="histpK0SvspSgn";
819       nameBkg="histpK0SvspBkg";
820       TH2F* momentumDistributionK0SvspSgn= new TH2F(nameSgn.Data(),"#Lambda_{c} - sgn: K^{0}_{S} vs p Total Momentum Distribution - MC; p_{p}; p_{K^{0}_{S}}",
821                                                     175,0.,35.,175,0.,35.);
822       TH2F* momentumDistributionK0SvspBkg= new TH2F(nameBkg.Data(),"#Lambda_{c} - bkg: K^{0}_{S} vs p Total Momentum Distribution - MC; p_{p}; p_{K^{0}_{S}}",
823                                                     175,0.,35.,175,0.,35.);
824
825       TH2F* allmomentumDistributionK0SvspSgn= (TH2F*)momentumDistributionK0SvspSgn->Clone(); 
826       TH2F* allmomentumDistributionK0SvspBkg= (TH2F*)momentumDistributionK0SvspBkg->Clone(); 
827       TH2F* pidBachmomentumDistributionK0SvspSgn= (TH2F*)momentumDistributionK0SvspSgn->Clone(); 
828       TH2F* pidBachmomentumDistributionK0SvspBkg= (TH2F*)momentumDistributionK0SvspBkg->Clone(); 
829       fOutputAll->Add(allmomentumDistributionK0SvspSgn); 
830       fOutputAll->Add(allmomentumDistributionK0SvspBkg); 
831       fOutputPIDBach->Add(pidBachmomentumDistributionK0SvspSgn); 
832       fOutputPIDBach->Add(pidBachmomentumDistributionK0SvspBkg); 
833
834
835       // armenteros-podolanski plots K0S
836       nameSgn="histArmPodK0SSgn";
837       nameBkg="histArmPodK0SBkg";
838       TH2F* armenterosPodK0SSgn = new TH2F(nameSgn.Data(),"K^{0}_{S}  Armenteros-Podolanski distribution (sgn); #frac{p_{L}^{+}-p_{L}^{-}}{p_{L}^{+}+p_{L}^{-}}; p_{T}^{+} [GeV/c]",
839                                            200,-1.,1.,300,0.,0.3);
840       TH2F* armenterosPodK0SBkg = new TH2F(nameBkg.Data(),"K^{0}_{S}  Armenteros-Podolanski distribution (bkg); #frac{p_{L}^{+}-p_{L}^{-}}{p_{L}^{+}+p_{L}^{-}}; p_{T}^{+} [GeV/c]",
841                                            200,-1.,1.,300,0.,0.3);
842
843       TH2F* allArmenterosPodK0SSgn = (TH2F*)armenterosPodK0SSgn->Clone();
844       TH2F* allArmenterosPodK0SBkg = (TH2F*)armenterosPodK0SBkg->Clone();
845       TH2F* pidBachArmenterosPodK0SSgn = (TH2F*)armenterosPodK0SSgn->Clone();
846       TH2F* pidBachArmenterosPodK0SBkg = (TH2F*)armenterosPodK0SBkg->Clone();
847
848       fOutputAll->Add(allArmenterosPodK0SSgn);
849       fOutputAll->Add(allArmenterosPodK0SBkg);
850       fOutputPIDBach->Add(pidBachArmenterosPodK0SSgn);
851       fOutputPIDBach->Add(pidBachArmenterosPodK0SBkg);
852
853
854       nameSgn="histDCAtoPVvspK0SSgn";
855       nameBkg="histDCAtoPVvspK0SBkg";
856       TH2F *dcatoPVvspK0SSgn=new TH2F(nameSgn.Data(),"K^{0}_{S} - sgn: DCA to Primary Vertex vs  K^{0}_{S} invariant mass (sgn); p(K^{0}_{S}) [GeV/c]; DCA to Primary Vertex [n#sigmas]; Entries",175,0.,35.,50,0.,5.);
857       TH2F *dcatoPVvspK0SBkg=new TH2F(nameBkg.Data(),"K^{0}_{S} - bkg: DCA to Primary Vertex vs  K^{0}_{S} invariant mass (bkg); p(K^{0}_{S}) [GeV/c]; DCA to Primary Vertex [n#sigmas]; Entries",175,0.,35.,50,0.,5.);
858
859       TH2F* alldcatoPVvspK0SSgn= (TH2F*)dcatoPVvspK0SSgn->Clone();
860       TH2F* alldcatoPVvspK0SBkg= (TH2F*)dcatoPVvspK0SBkg->Clone();
861       TH2F* pidBachdcatoPVvspK0SSgn= (TH2F*)dcatoPVvspK0SSgn->Clone();
862       TH2F* pidBachdcatoPVvspK0SBkg= (TH2F*)dcatoPVvspK0SBkg->Clone();
863
864       fOutputAll->Add(alldcatoPVvspK0SSgn);
865       fOutputPIDBach->Add(pidBachdcatoPVvspK0SSgn);
866       fOutputAll->Add(alldcatoPVvspK0SBkg);
867       fOutputPIDBach->Add(pidBachdcatoPVvspK0SBkg);
868
869     }
870
871
872     nameSgn="histK0SMassOfflineSgn";
873     nameBkg="histK0SMassOfflineBkg";
874     TH2F* spectrumK0SMassOfflineSgn = new TH2F(nameSgn.Data(), "K^{0}_{S} - sgn: invariant mass VS p_{T} - MC; m_{inv}(#pi^{+}-#pi^{-}) [GeV/c^{2}]; p_{T}(K^{0}_{S}) [GeV/c]; Entries",
875                                               1000,mK0SPDG-0.050,mK0SPDG+0.050,175,0.,35.);
876     TH2F* spectrumK0SMassOfflineBkg = new TH2F(nameBkg.Data(), "K^{0}_{S} - bkg: invariant mass VS p_{T} - MC; m_{inv}(#pi^{+}-#pi^{-}) [GeV/c^{2}]; p_{T}(K^{0}_{S}) [GeV/c]; Entries",
877                                               1000,mK0SPDG-0.050,mK0SPDG+0.050,175,0.,35.);
878
879     TH2F* allspectrumK0SMassOfflineSgn = (TH2F*)spectrumK0SMassOfflineSgn->Clone(); 
880     TH2F* allspectrumK0SMassOfflineBkg = (TH2F*) spectrumK0SMassOfflineBkg->Clone();  
881     fOutputAll->Add(allspectrumK0SMassOfflineSgn);
882     fOutputAll->Add(allspectrumK0SMassOfflineBkg);
883
884
885     TH2F* pidBachspectrumK0SMassOfflineSgn = (TH2F*)spectrumK0SMassOfflineSgn->Clone(); 
886     TH2F* pidBachspectrumK0SMassOfflineBkg = (TH2F*) spectrumK0SMassOfflineBkg->Clone();
887     fOutputPIDBach->Add(pidBachspectrumK0SMassOfflineSgn);
888     fOutputPIDBach->Add(pidBachspectrumK0SMassOfflineBkg);
889
890
891     nameSgn="histLcMassByK0SOfflineSgn";
892     nameBkg="histLcMassByK0SOfflineBkg";
893     TH2F* spectrumLcMassOfflineByK0SSgn = new TH2F(nameSgn.Data(), "#Lambda_{c} - sgn: invariant mass (by K^{0}_{S})  vs p_{T} - MC; M(#Lambda_{c}) [GeV/c^{2}]; p_{T}",
894                                                    1000,mLcPDG-0.250,mLcPDG+0.250,175,0.,35.);
895     TH2F* spectrumLcMassOfflineByK0SBkg = new TH2F(nameBkg.Data(), "#Lambda_{c} - bkg: invariant mass (by K^{0}_{S})  vs p_{T} - MC; M(#Lambda_{c}) [GeV/c^{2}]; p_{T}",
896                                                    1000,mLcPDG-0.250,mLcPDG+0.250,175,0.,35.);
897
898
899     TH2F* allspectrumLcMassOfflineByK0SSgn = (TH2F*)spectrumLcMassOfflineByK0SSgn->Clone(); 
900     TH2F* allspectrumLcMassOfflineByK0SBkg = (TH2F*) spectrumLcMassOfflineByK0SBkg->Clone();  
901     TH2F* pidBachspectrumLcMassOfflineByK0SSgn = (TH2F*)spectrumLcMassOfflineByK0SSgn->Clone(); 
902     TH2F* pidBachspectrumLcMassOfflineByK0SBkg = (TH2F*) spectrumLcMassOfflineByK0SBkg->Clone();  
903     fOutputAll->Add(allspectrumLcMassOfflineByK0SSgn);
904     fOutputAll->Add(allspectrumLcMassOfflineByK0SBkg);
905     fOutputPIDBach->Add(pidBachspectrumLcMassOfflineByK0SSgn);
906     fOutputPIDBach->Add(pidBachspectrumLcMassOfflineByK0SBkg);
907   
908  
909     nameSgn="histpK0SvspOfflineSgn";
910     nameBkg="histpK0SvspOfflineBkg";
911     TH2F* momentumDistributionK0SvspOfflineSgn= new TH2F(nameSgn.Data(),"#Lambda_{c} - sgn: K^{0}_{S} vs p Total Momentum Distribution - Offline  - MC; p_{p}; p_{K^{0}_{S}}",
912                                                          175,0.,35.,175,0.,35.);
913     TH2F* momentumDistributionK0SvspOfflineBkg= new TH2F(nameBkg.Data(),"#Lambda_{c} - bkg: K^{0}_{S} vs p Total Momentum Distribution - Offline  - MC; p_{p}; p_{K^{0}_{S}}",
914                                                          175,0.,35.,175,0.,35.);
915
916
917     TH2F* allmomentumDistributionK0SvspOfflineSgn= (TH2F*)momentumDistributionK0SvspOfflineSgn->Clone(); 
918     TH2F* allmomentumDistributionK0SvspOfflineBkg= (TH2F*)momentumDistributionK0SvspOfflineBkg->Clone(); 
919     TH2F* pidBachmomentumDistributionK0SvspOfflineSgn= (TH2F*)momentumDistributionK0SvspOfflineSgn->Clone(); 
920     TH2F* pidBachmomentumDistributionK0SvspOfflineBkg= (TH2F*)momentumDistributionK0SvspOfflineBkg->Clone(); 
921     fOutputAll->Add(allmomentumDistributionK0SvspOfflineSgn); 
922     fOutputAll->Add(allmomentumDistributionK0SvspOfflineBkg); 
923     fOutputPIDBach->Add(pidBachmomentumDistributionK0SvspOfflineSgn); 
924     fOutputPIDBach->Add(pidBachmomentumDistributionK0SvspOfflineBkg); 
925
926
927
928
929
930     // armenteros-podolanski plots K0S (offline)
931     nameSgn="histArmPodK0SOfflineSgn";
932     nameBkg="histArmPodK0SOfflineBkg";
933     TH2F* armenterosPodK0SOffSgn = new TH2F(nameSgn.Data(),"K^{0}_{S}  Armenteros-Podolanski distribution (sgn) -offline-; #frac{p_{L}^{+}-p_{L}^{-}}{p_{L}^{+}+p_{L}^{-}}; p_{T}^{+} [GeV/c]",
934                                             200,-1.,1.,300,0.,0.3);
935     TH2F* armenterosPodK0SOffBkg = new TH2F(nameBkg.Data(),"K^{0}_{S}  Armenteros-Podolanski distribution (bkg) -offline-; #frac{p_{L}^{+}-p_{L}^{-}}{p_{L}^{+}+p_{L}^{-}}; p_{T}^{+} [GeV/c]",
936                                             200,-1.,1.,300,0.,0.3);
937   
938
939     TH2F* allArmenterosPodK0SOffSgn = (TH2F*)armenterosPodK0SOffSgn->Clone();
940     TH2F* allArmenterosPodK0SOffBkg = (TH2F*)armenterosPodK0SOffBkg->Clone();
941     TH2F* pidBachArmenterosPodK0SOffSgn = (TH2F*)armenterosPodK0SOffSgn->Clone();
942     TH2F* pidBachArmenterosPodK0SOffBkg = (TH2F*)armenterosPodK0SOffBkg->Clone();
943
944
945     fOutputAll->Add(allArmenterosPodK0SOffSgn);
946     fOutputAll->Add(allArmenterosPodK0SOffBkg);
947     fOutputPIDBach->Add(pidBachArmenterosPodK0SOffSgn);
948     fOutputPIDBach->Add(pidBachArmenterosPodK0SOffBkg);
949
950
951     nameSgn="histDCAtoPVvspK0SOfflineSgn";
952     nameBkg="histDCAtoPVvspK0SOfflineBkg";
953     TH2F *dcatoPVvspK0SOfflineSgn=new TH2F(nameSgn.Data(),"K^{0}_{S} -offline - (sgn): DCA to Primary Vertex vs  K^{0}_{S} invariant mass; p(K^{0}_{S}) [GeV/c]; DCA to Primary Vertex [n#sigmas]; Entries",
954                                            175,0.,35.,50,0.,5.);
955     TH2F *dcatoPVvspK0SOfflineBkg=new TH2F(nameBkg.Data(),"K^{0}_{S} -offline - (bkg): DCA to Primary Vertex vs  K^{0}_{S} invariant mass; p(K^{0}_{S}) [GeV/c]; DCA to Primary Vertex [n#sigmas]; Entries",
956                                            175,0.,35.,50,0.,5.);
957     
958
959     TH2F* alldcatoPVvspK0SOfflineSgn= (TH2F*)dcatoPVvspK0SOfflineSgn->Clone();
960     TH2F* pidBachdcatoPVvspK0SOfflineSgn= (TH2F*)dcatoPVvspK0SOfflineSgn->Clone();
961     TH2F* alldcatoPVvspK0SOfflineBkg= (TH2F*)dcatoPVvspK0SOfflineBkg->Clone();
962     TH2F* pidBachdcatoPVvspK0SOfflineBkg= (TH2F*)dcatoPVvspK0SOfflineBkg->Clone();
963
964
965
966     fOutputAll->Add(alldcatoPVvspK0SOfflineSgn);
967     fOutputPIDBach->Add(pidBachdcatoPVvspK0SOfflineSgn);
968     fOutputAll->Add(alldcatoPVvspK0SOfflineBkg);
969     fOutputPIDBach->Add(pidBachdcatoPVvspK0SOfflineBkg);
970
971   }
972
973   return;
974 }
975
976 //---------------------------
977 void AliAnalysisTaskSELc2V0bachelor::CheckEventSelection(AliAODEvent *aodEvent) {
978   //
979   // To fill control histograms
980   //
981
982   TClonesArray *arrayLctopKos=0;
983   if (!aodEvent){
984     if(AODEvent() && IsStandardAOD()) {
985       // In case there is an AOD handler writing a standard AOD, use the AOD 
986       // event in memory rather than the input (ESD) event.    
987       aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
988       // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
989       // have to taken from the AOD event hold by the AliAODExtension
990       AliAODHandler* aodHandler = (AliAODHandler*) 
991         ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
992
993       if (aodHandler->GetExtensions()) {
994         AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
995         AliAODEvent *aodFromExt = ext->GetAOD();
996         arrayLctopKos=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
997       }
998     }
999   } else {
1000     arrayLctopKos=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
1001   }
1002
1003   Float_t zVertex = fVtx1->GetZ();
1004   TString titleVtx=fVtx1->GetTitle();
1005
1006   if (TMath::Abs(fBzkG)>=0.001) {
1007
1008     if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ2")))->Fill(zVertex);
1009
1010     if (arrayLctopKos) {
1011
1012       if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ3")))->Fill(zVertex);
1013
1014       // mc analysis 
1015       TClonesArray *mcArray = 0;
1016       AliAODMCHeader *mcHeader=0;
1017
1018       if (fUseMCInfo) {
1019         // MC array need for maching
1020         mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1021
1022         if (mcArray) {
1023           // load MC header
1024           if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ4")))->Fill(zVertex);
1025           mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1026
1027           if (mcHeader && fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ5")))->Fill(zVertex);
1028
1029           // check on MC Lc Daughter
1030           if (fAdditionalChecks) {
1031             for (Int_t iii=0; iii<mcArray->GetEntries(); iii++)
1032               SearchLcDaughter(mcArray,iii);
1033           }
1034
1035         }
1036
1037       }
1038
1039       if (fVtx1->GetNContributors()>0) {
1040         if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ6")))->Fill(zVertex);
1041
1042         TString firedTriggerClasses = aodEvent->GetFiredTriggerClasses(); // trigger class
1043         ULong64_t fTriggerMask=AliVEvent::kAnyINT;
1044         Bool_t check1 = kFALSE;
1045         if ( !fUseMCInfo && // don't do for MC...
1046              (aodEvent->GetRunNumber()<136851 || aodEvent->GetRunNumber()>139517) ) { // ...and for PbPb 2010 data
1047           if ( !(firedTriggerClasses.Contains("CINT1")) ) {
1048             AliInfo(Form(" ======================== firedTriggerClasses.Data() = %s",firedTriggerClasses.Data()));
1049             fCEvents->Fill(8);
1050             if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ8")))->Fill(zVertex);
1051             check1 = kTRUE;
1052           }
1053         }
1054
1055         Bool_t isSelectedAAA = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
1056         if (!isSelectedAAA) {
1057           fCEvents->Fill(9);
1058           if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ9")))->Fill(zVertex);
1059         }
1060
1061         if (!isSelectedAAA || check1) {
1062           fCEvents->Fill(16);
1063           if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ16")))->Fill(zVertex);
1064         }
1065
1066         fTriggerMask=AliVEvent::kAny;
1067         Bool_t isSelectedBBB = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
1068         if (!isSelectedBBB) {
1069           fCEvents->Fill(10);
1070           if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ10")))->Fill(zVertex);
1071         }
1072
1073         if (titleVtx.Contains("Z")) {
1074           fCEvents->Fill(11);
1075           if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ11")))->Fill(zVertex);
1076         }
1077         else if (titleVtx.Contains("3D")) {
1078           fCEvents->Fill(12);
1079           if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ12")))->Fill(zVertex);
1080         } else {
1081           fCEvents->Fill(13);
1082           if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ13")))->Fill(zVertex);
1083         }
1084
1085         if (TMath::Abs(zVertex)<=fAnalCuts->GetMaxVtxZ()) {
1086           fCEvents->Fill(14);
1087           if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ14")))->Fill(zVertex);
1088         }
1089
1090         if ( fIsEventSelected ) {
1091           if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ7")))->Fill(zVertex);
1092         } else {
1093           fCEvents->Fill(15);
1094           if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ15")))->Fill(zVertex);
1095         }
1096
1097       } // nContributors>=1
1098     } // analysisArray exists
1099   } // magnetic field exists
1100
1101   return;
1102 }
1103
1104 //-----------------
1105 void AliAnalysisTaskSELc2V0bachelor::CheckEventSelectionWithCandidates(AliAODEvent *aodEvent) {
1106   //
1107   // To fill control histograms
1108   //
1109
1110   Float_t zVertex = fVtx1->GetZ();
1111   TString titleVtx=fVtx1->GetTitle();
1112   TString firedTriggerClasses = aodEvent->GetFiredTriggerClasses(); // trigger class
1113   ULong64_t fTriggerMask=AliVEvent::kAnyINT;
1114
1115   ((TH1F*)(fOutput->FindObject("hEventsWithCandidates")))->Fill(6);
1116   if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ6a")))->Fill(zVertex);
1117
1118   Bool_t check1a = kFALSE;
1119   if ( !fUseMCInfo && // don't do for MC...
1120        (aodEvent->GetRunNumber()<136851 || aodEvent->GetRunNumber()>139517) ) { // ...and for PbPb 2010 data
1121     if ( !(firedTriggerClasses.Contains("CINT1")) ) {
1122       AliInfo(Form(" ======================== firedTriggerClasses.Data() = %s",firedTriggerClasses.Data()));
1123       ((TH1F*)(fOutput->FindObject("hEventsWithCandidates")))->Fill(8);
1124       if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ8a")))->Fill(zVertex);
1125       check1a = kTRUE;
1126     }
1127   }
1128
1129   Bool_t isSelectedAAAa = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
1130   if (!isSelectedAAAa) {
1131     ((TH1F*)(fOutput->FindObject("hEventsWithCandidates")))->Fill(9);
1132     if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ9a")))->Fill(zVertex);
1133   }
1134
1135   if (!isSelectedAAAa || check1a) {
1136     ((TH1F*)(fOutput->FindObject("hEventsWithCandidates")))->Fill(16);
1137     if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ16a")))->Fill(zVertex);
1138   }
1139
1140   fTriggerMask=AliVEvent::kAny;
1141   Bool_t isSelectedBBBa = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
1142   if (!isSelectedBBBa) {
1143     ((TH1F*)(fOutput->FindObject("hEventsWithCandidates")))->Fill(10);
1144     if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ10a")))->Fill(zVertex);
1145   }
1146
1147   if (titleVtx.Contains("Z")) {
1148     ((TH1F*)(fOutput->FindObject("hEventsWithCandidates")))->Fill(11);
1149     if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ11a")))->Fill(zVertex);
1150   }
1151   else if (titleVtx.Contains("3D")) {
1152     ((TH1F*)(fOutput->FindObject("hEventsWithCandidates")))->Fill(12);
1153     if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ12a")))->Fill(zVertex);
1154   } else {
1155     ((TH1F*)(fOutput->FindObject("hEventsWithCandidates")))->Fill(13);
1156     if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ13a")))->Fill(zVertex);
1157   }
1158
1159   if (TMath::Abs(zVertex)<=fAnalCuts->GetMaxVtxZ()) {
1160     ((TH1F*)(fOutput->FindObject("hEventsWithCandidates")))->Fill(14);
1161     if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ14a")))->Fill(zVertex);
1162   }
1163
1164   if ( fIsEventSelected ) {
1165     ((TH1F*)(fOutput->FindObject("hEventsWithCandidates")))->Fill(7);
1166     if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ7a")))->Fill(zVertex);
1167   } else {
1168     ((TH1F*)(fOutput->FindObject("hEventsWithCandidates")))->Fill(15);
1169     if (fAdditionalChecks) ((TH1F*)(fOutput->FindObject("hZ15a")))->Fill(zVertex);
1170   }
1171
1172   return;
1173 }
1174
1175 //-----------------------------
1176 Int_t AliAnalysisTaskSELc2V0bachelor::MatchToMC(AliAODRecoCascadeHF *lc2bacV0,
1177                                                 Int_t *pdgDgLc2bacV0, Int_t *pdgDgV0,
1178                                                 TClonesArray *mcArray) {
1179   //
1180   // This is now implemented in AliAODRecoCascadeHF
1181   //
1182
1183   // bachelor
1184   AliAODTrack *bachelor = (AliAODTrack*)lc2bacV0->GetBachelor();
1185   if (!bachelor) return -1;
1186   Int_t labBachelor = TMath::Abs(bachelor->GetLabel());
1187   if (labBachelor<0) return -1;
1188   AliAODMCParticle *partBachelor = (AliAODMCParticle*)mcArray->At(labBachelor);
1189   if (!partBachelor) return -1;
1190   if (TMath::Abs(partBachelor->GetPdgCode())!=pdgDgLc2bacV0[0]) return -1;
1191
1192   Int_t labBacMother = partBachelor->GetMother();
1193   if (labBacMother<0) return -1;
1194   AliAODMCParticle *partBacMother = (AliAODMCParticle*)mcArray->At(labBacMother);
1195   if (!partBacMother) return -1;
1196   if (TMath::Abs(partBacMother->GetPdgCode())!=4122) return -1;
1197
1198   // V0
1199   AliAODTrack *posV0Daugh = (AliAODTrack*)lc2bacV0->Getv0PositiveTrack();
1200   AliAODTrack *negV0Daugh = (AliAODTrack*)lc2bacV0->Getv0NegativeTrack();
1201   if (!posV0Daugh || !negV0Daugh) return -1;
1202
1203   Int_t labV0pos = TMath::Abs(posV0Daugh->GetLabel());
1204   Int_t labV0neg = TMath::Abs(negV0Daugh->GetLabel());
1205   if (labV0pos<0 || labV0neg<0) return -1;
1206
1207   AliAODMCParticle *partV0pos = (AliAODMCParticle*)mcArray->At(labV0neg);
1208   AliAODMCParticle *partV0neg = (AliAODMCParticle*)mcArray->At(labV0pos);
1209   if (!partV0pos || !partV0neg) return -1;
1210
1211   if ( ! ( (TMath::Abs(partV0pos->GetPdgCode())==pdgDgV0[0] &&
1212             TMath::Abs(partV0neg->GetPdgCode())==pdgDgV0[1]) ||
1213            (TMath::Abs(partV0pos->GetPdgCode())==pdgDgV0[1] &&
1214             TMath::Abs(partV0neg->GetPdgCode())==pdgDgV0[0]) ) ) return -1;
1215   Int_t labV0posMother = partV0pos->GetMother();
1216   Int_t labV0negMother = partV0neg->GetMother();
1217
1218   if (labV0posMother<0 || labV0negMother<0) return -1;
1219   if (labV0posMother!=labV0negMother) return -1;
1220
1221   AliAODMCParticle *motherV0 = (AliAODMCParticle*)mcArray->At(labV0posMother);
1222   if (!motherV0) return-1;
1223
1224   if (TMath::Abs(motherV0->GetPdgCode())!=pdgDgLc2bacV0[1]) return -1;
1225   Int_t labV0mother = motherV0->GetMother();
1226   if (labV0mother<0) return -1;
1227   AliAODMCParticle *gMotherV0 = (AliAODMCParticle*)mcArray->At(labV0mother);
1228   if (!gMotherV0) return-1;
1229
1230   if ( !(pdgDgLc2bacV0[1]==310 && TMath::Abs(gMotherV0->GetPdgCode())==311) &&
1231        !(pdgDgLc2bacV0[1]==3122 && TMath::Abs(motherV0->GetPdgCode())==3122) ) return -1;
1232
1233   if ( (pdgDgLc2bacV0[1]==310 && TMath::Abs(gMotherV0->GetPdgCode())==311) ) {
1234     Int_t labV0GMother = gMotherV0->GetMother();
1235     if (labV0GMother<0) return -1;
1236     AliAODMCParticle *ggMotherV0 = (AliAODMCParticle*)mcArray->At(labV0GMother);
1237     if (!ggMotherV0) return-1;
1238
1239     if (TMath::Abs(ggMotherV0->GetPdgCode())!=4122) return -1;
1240     gMotherV0 = (AliAODMCParticle*)ggMotherV0;
1241     labV0mother=labV0GMother;
1242   }
1243   else if (pdgDgLc2bacV0[1]==3122 && TMath::Abs(motherV0->GetPdgCode())==3122) {
1244     if (TMath::Abs(gMotherV0->GetPdgCode())!=4122) return -1;
1245   }
1246
1247   if (labBacMother!=labV0mother) {
1248     return -1;
1249   }
1250
1251   return labBacMother;
1252
1253 }
1254
1255 //________________________________________________________________
1256 Int_t AliAnalysisTaskSELc2V0bachelor::SearchLcDaughter(TClonesArray *arrayMC, Int_t iii) {
1257   //
1258   // This is to check Lc dinasty
1259   //
1260
1261   Int_t indexToBeReturned=-999;
1262
1263   Int_t pdgLc=4122;
1264   Int_t pdgLambda=3122;
1265   Int_t pdgV0=310;
1266   Int_t pdgK0=311;
1267   Int_t pdgBachelor=2212;
1268   Int_t pdgBachelorPi=211;
1269
1270   TString fillthis="";
1271   fillthis="histMcStatLc";
1272
1273   AliAODMCParticle *searchLc = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iii));
1274   if(!searchLc) return -999;
1275   if (TMath::Abs(searchLc->GetPdgCode()) != pdgLc) return -999;
1276
1277   ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(0);
1278   indexToBeReturned = 0;
1279
1280   ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*1);
1281   indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*1;
1282
1283   Int_t nDaughLc = searchLc->GetNDaughters();
1284   if (nDaughLc!=2) {
1285     ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*10);
1286     indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*10;
1287     return indexToBeReturned;
1288   }
1289
1290   Int_t index1=searchLc->GetDaughter(0);
1291   Int_t index2=searchLc->GetDaughter(1);
1292   if (index1<=0 || index2<=0) {
1293     return -999;
1294   }
1295
1296   AliAODMCParticle *daugh1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(index1));
1297   AliAODMCParticle *daugh2 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(index2));
1298   if (!daugh1 || !daugh2) return -999;
1299
1300   Int_t daughPdg1 = TMath::Abs(daugh1->GetPdgCode());
1301   Int_t daughPdg2 = TMath::Abs(daugh2->GetPdgCode());
1302   if ( !( (daughPdg1==pdgBachelor && daughPdg2==pdgK0) ||
1303           (daughPdg2==pdgBachelor && daughPdg1==pdgK0) ||
1304           (daughPdg1==pdgLambda && daughPdg2==pdgBachelorPi) ||
1305           (daughPdg2==pdgLambda && daughPdg1==pdgBachelorPi) ) ) {
1306     ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*10);
1307     indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*10;
1308     return indexToBeReturned;
1309   }
1310
1311   if (daughPdg1==pdgK0 || daughPdg1==pdgLambda) {
1312     index1=searchLc->GetDaughter(1);
1313     index2=searchLc->GetDaughter(0);
1314   }
1315   daugh1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(index1));
1316   daugh2 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(index2));
1317   if (!daugh1 || !daugh2) return -999;
1318
1319   daughPdg1=TMath::Abs(daugh1->GetPdgCode());
1320   daughPdg2=TMath::Abs(daugh2->GetPdgCode());
1321
1322   if ( daughPdg1==pdgBachelor && daughPdg2==pdgK0 ) { // Lc+ -> p K0bar
1323
1324     ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*2);
1325     indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*2;
1326
1327     Int_t nDaughK0 = daugh2->GetNDaughters();
1328     if (nDaughK0!=1) return -999;
1329
1330     Int_t indexK0daugh=daugh2->GetDaughter(0);
1331     if (indexK0daugh<=0) return -999;
1332
1333     AliAODMCParticle *daughK0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indexK0daugh));
1334     if (!daughK0) return -999;
1335
1336     Int_t daughK0Pdg=TMath::Abs(daughK0->GetPdgCode());
1337     if (daughK0Pdg!=pdgV0) {
1338       ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*4); // K0L
1339       indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*4;
1340       return indexToBeReturned;
1341     }
1342     ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*3); // K0S
1343     indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*3;
1344
1345     Int_t nDaughK0S = daughK0->GetNDaughters();
1346     if (nDaughK0S!=2) {
1347       ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*5); // other decays for K0S
1348       indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*5;
1349       return indexToBeReturned;
1350     }
1351
1352     index1=daughK0->GetDaughter(0);
1353     index2=daughK0->GetDaughter(1);
1354     if(index1<=0 || index2<=0) {
1355       return -999;
1356     }
1357
1358     AliAODMCParticle *daughK0S1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(index1));
1359     AliAODMCParticle *daughK0S2 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(index2));
1360     if (!daughK0S1 || !daughK0S2) return -999;
1361
1362     Int_t daughK0S1pdg=TMath::Abs(daughK0S1->GetPdgCode());
1363     Int_t daughK0S2pdg=TMath::Abs(daughK0S2->GetPdgCode());
1364
1365     if ( daughK0S1pdg==211 && daughK0S2pdg==211 ) {
1366       ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*6); // K0S -> pi+ pi-
1367       indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*6;
1368     } else {
1369       ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*5); // other decays for K0S
1370       indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*5;
1371     }
1372
1373   } //if (daughPdg1==pdgBachelor && daughPdg2==pdgK0)
1374   else if ( daughPdg1==pdgBachelorPi && daughPdg2==pdgLambda ) { // Lc+ -> pi+ Lambda
1375
1376     ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*7);
1377     indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*7;
1378
1379     Int_t nDaughL = daugh2->GetNDaughters();
1380     if (nDaughL!=2) {
1381       ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*8);
1382       indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*8;
1383       return indexToBeReturned;
1384     }
1385
1386     index1=daugh2->GetDaughter(0);
1387     index2=daugh2->GetDaughter(1);
1388     if(index1<=0 || index2<=0) {
1389       return -999;
1390     }
1391
1392     AliAODMCParticle *daughL1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(index1));
1393     AliAODMCParticle *daughL2 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(index2));
1394     if (!daughL1 || !daughL2) return -999;
1395
1396     Int_t daughL1pdg=TMath::Abs(daughL1->GetPdgCode());
1397     Int_t daughL2pdg=TMath::Abs(daughL2->GetPdgCode());
1398     if ( (daughL1pdg==211 && daughL2pdg==2212) ||
1399          (daughL2pdg==211 && daughL1pdg==2212) ) {
1400       ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*9);
1401       indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*9;
1402     } else {
1403       ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(TMath::Nint(searchLc->Charge()/3.)*8);
1404       indexToBeReturned = TMath::Nint(searchLc->Charge()/3.)*8;
1405     }
1406
1407   } //else if (daughPdg1==pdgBachelorPi && daughPdg2==pdgLambda)
1408
1409   return indexToBeReturned;
1410 }
1411
1412 //________________________________________________________________
1413 void AliAnalysisTaskSELc2V0bachelor::FillArmPodDistribution(AliAODv0 *vZero,
1414                                                             TString histoTitle,
1415                                                             TList *histoList) {
1416   //
1417   // This is to fill Armenteros Podolanski plots
1418   //
1419
1420   Double_t alpha = vZero->AlphaV0();
1421   Double_t qT    = vZero->PtArmV0();
1422
1423   ((TH2F*)(histoList->FindObject(histoTitle)))->Fill(alpha,qT);
1424
1425 }
1426
1427 //-------------------------------------------------------------------------------
1428 void AliAnalysisTaskSELc2V0bachelor::CheckCandidatesAtDifferentLevels(AliAODRecoCascadeHF *part, AliRDHFCutsLctoV0* cutsAnal) {
1429   //
1430   // This is to check candidates at different levels
1431   //
1432
1433   Bool_t areCutsUsingPID = cutsAnal->GetIsUsePID();
1434
1435   AliAODv0 * v0part = (AliAODv0*)part->Getv0();
1436   Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
1437
1438   AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
1439
1440   if ( !onFlyV0 )
1441     ((TH1F*)(fOutput->FindObject("hCandidateSelection")))->Fill(3); // it counts number of candidates coming from offline V0s
1442
1443   if ( cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122)) )
1444     ((TH1F*)(fOutput->FindObject("hCandidateSelection")))->Fill(4);
1445   if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kTracks))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) )
1446     ((TH1F*)(fOutput->FindObject("hCandidateSelection")))->Fill(5);
1447   cutsAnal->SetUsePID(kFALSE);
1448   if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kCandidate))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) )
1449     ((TH1F*)(fOutput->FindObject("hCandidateSelection")))->Fill(6);
1450   cutsAnal->SetUsePID(areCutsUsingPID);
1451   if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) )
1452     ((TH1F*)(fOutput->FindObject("hCandidateSelection")))->Fill(7);
1453   if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kCandidate))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) )
1454     ((TH1F*)(fOutput->FindObject("hCandidateSelection")))->Fill(8);
1455   if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kAll))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) )
1456     ((TH1F*)(fOutput->FindObject("hCandidateSelection")))->Fill(9);
1457
1458   if ( cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122)) ) {
1459
1460     if ( ( (cutsAnal->IsSelected(part,AliRDHFCuts::kTracks))&(AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) {
1461
1462       Int_t aaa = cutsAnal->IsSelected(part,AliRDHFCuts::kTracks);
1463       if ( (aaa&AliRDHFCutsLctoV0::kLcToK0Spr)==AliRDHFCutsLctoV0::kLcToK0Spr ) {
1464         if ( ( (aaa&AliRDHFCutsLctoV0::kLcToLpi)==AliRDHFCutsLctoV0::kLcToLpi && bachelor->Charge()==-1)  ||
1465              ( (aaa&AliRDHFCutsLctoV0::kLcToLBarpi)==AliRDHFCutsLctoV0::kLcToLBarpi && bachelor->Charge()==+1) )
1466           ((TH1F*)(fOutput->FindObject("hSwitchOnCandidates1")))->Fill( -aaa );
1467         else
1468           ((TH1F*)(fOutput->FindObject("hSwitchOnCandidates1")))->Fill( aaa );
1469       }
1470
1471       cutsAnal->SetUsePID(kFALSE);
1472       aaa = cutsAnal->IsSelected(part,AliRDHFCuts::kCandidate);
1473       if ((aaa&AliRDHFCutsLctoV0::kLcToK0Spr)==AliRDHFCutsLctoV0::kLcToK0Spr) {
1474         if ( ( (aaa&AliRDHFCutsLctoV0::kLcToLpi)==AliRDHFCutsLctoV0::kLcToLpi && bachelor->Charge()==-1) ||
1475              ( (aaa&AliRDHFCutsLctoV0::kLcToLBarpi)==AliRDHFCutsLctoV0::kLcToLBarpi && bachelor->Charge()==+1) )
1476           ((TH1F*)(fOutput->FindObject("hSwitchOnCandidates2")))->Fill( -aaa );
1477         else
1478           ((TH1F*)(fOutput->FindObject("hSwitchOnCandidates2")))->Fill( aaa );
1479       }
1480       cutsAnal->SetUsePID(areCutsUsingPID);
1481
1482       aaa = cutsAnal->IsSelected(part,AliRDHFCuts::kPID);
1483       if ((aaa&AliRDHFCutsLctoV0::kLcToK0Spr)==AliRDHFCutsLctoV0::kLcToK0Spr) {
1484         if ( ( (aaa&AliRDHFCutsLctoV0::kLcToLpi)==AliRDHFCutsLctoV0::kLcToLpi && bachelor->Charge()==-1) ||
1485              ( (aaa&AliRDHFCutsLctoV0::kLcToLBarpi)==AliRDHFCutsLctoV0::kLcToLBarpi && bachelor->Charge()==+1) )
1486           ((TH1F*)(fOutput->FindObject("hSwitchOnCandidates3")))->Fill( -aaa );
1487         else
1488           ((TH1F*)(fOutput->FindObject("hSwitchOnCandidates3")))->Fill( aaa );
1489       }
1490
1491       aaa = cutsAnal->IsSelected(part,AliRDHFCuts::kAll);
1492       if ((aaa&AliRDHFCutsLctoV0::kLcToK0Spr)==AliRDHFCutsLctoV0::kLcToK0Spr) {
1493         if ( ( (aaa&AliRDHFCutsLctoV0::kLcToLpi)==AliRDHFCutsLctoV0::kLcToLpi && bachelor->Charge()==-1) ||
1494              ( (aaa&AliRDHFCutsLctoV0::kLcToLBarpi)==AliRDHFCutsLctoV0::kLcToLBarpi && bachelor->Charge()==+1) )
1495           ((TH1F*)(fOutput->FindObject("hSwitchOnCandidates4")))->Fill( -aaa );
1496         else
1497           ((TH1F*)(fOutput->FindObject("hSwitchOnCandidates4")))->Fill( aaa );
1498       }
1499
1500     }
1501   }
1502
1503   return;
1504 }
1505
1506 //-------------------------------------------------------------------------------
1507 void AliAnalysisTaskSELc2V0bachelor::FillTheTree(AliAODRecoCascadeHF *part, AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *mcArray, Int_t isLc) {
1508   //
1509   // This is to fill tree
1510   //
1511
1512   Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
1513   Double_t mLPDG   = TDatabasePDG::Instance()->GetParticle(3122)->Mass();
1514
1515   Double_t invmassLc = part->InvMassLctoK0sP();
1516   Double_t invmassLc2Lpi = part->InvMassLctoLambdaPi();
1517
1518   AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
1519
1520   AliAODv0 * v0part = (AliAODv0*)part->Getv0();
1521   Double_t dcaV0ptp = v0part->GetDCA();
1522   Double_t invmassK0S = v0part->MassK0Short();
1523   Double_t invmassLambda = v0part->MassLambda();
1524   Double_t invmassLambdaBar = v0part->MassAntiLambda();
1525
1526   Int_t isLc2LBarpi=0, isLc2Lpi=0;
1527   Int_t mcLabel = -1;
1528   Int_t isDp2K0Spi=0, isDs2K0SK=0;
1529   Int_t mcLabel2 = -1;
1530   Int_t mcLabel3 = -1;
1531   Int_t isKstar12K0Spi=0, isKstar22K0Spi=0;
1532   Int_t mcLabel4 = -1;
1533   Int_t mcLabel5 = -1;
1534   Double_t ptCandByMC = 0.;//fmcPartCandidate->Pt();
1535   Double_t yCandByMC  = 0.;//fmcPartCandidate->Y() ;
1536   if (fUseMCInfo) {
1537     if (isLc) {
1538       Int_t pdgCand0 = 4122;
1539       Int_t pdgDgLctoV0bachelor0[2]={2212,310};
1540       Int_t pdgDgV0toDaughters0[2]={211,211};
1541       Int_t mcLabelLc2pK0S = part->MatchToMC(pdgCand0,pdgDgLctoV0bachelor0[1],pdgDgLctoV0bachelor0,pdgDgV0toDaughters0,mcArray,kTRUE);
1542       AliAODMCParticle *lambdaCpartMC = (AliAODMCParticle*)mcArray->At(mcLabelLc2pK0S);
1543       if (lambdaCpartMC) {
1544         ptCandByMC = lambdaCpartMC->Pt();
1545         yCandByMC  = lambdaCpartMC->Y() ;
1546       }
1547     }
1548
1549     Int_t pdgCand = 4122;
1550     Int_t pdgDgLctoV0bachelor[2]={211,3122};
1551     Int_t pdgDgV0toDaughters[2]={2212,211};
1552     mcLabel = part->MatchToMC(pdgCand,pdgDgLctoV0bachelor[1],pdgDgLctoV0bachelor,pdgDgV0toDaughters,mcArray,kTRUE);
1553     if (mcLabel!=-1) {
1554       if (bachelor->Charge()==-1) isLc2LBarpi=1;
1555       if (bachelor->Charge()==+1) isLc2Lpi=1;
1556       AliAODMCParticle *lambdaCpartMC = (AliAODMCParticle*)mcArray->At(mcLabel);
1557       if (lambdaCpartMC) {
1558         ptCandByMC = lambdaCpartMC->Pt();
1559         yCandByMC  = lambdaCpartMC->Y() ;
1560       }
1561     }
1562
1563     Int_t pdgCand2 = 411; // D+ -> pi+ K0S
1564     Int_t pdgCand3 = 431; // Ds+ -> K+ K0S
1565     Int_t pdgDgCand2[2]={211,310};
1566     Int_t pdgDgCand3[2]={321,310};
1567     pdgDgV0toDaughters[0]=211;
1568     pdgDgV0toDaughters[1]=211;
1569     mcLabel2 = part->MatchToMC(pdgCand2,pdgDgCand2[1],pdgDgCand2,pdgDgV0toDaughters,mcArray,kTRUE);
1570     mcLabel3 = part->MatchToMC(pdgCand3,pdgDgCand3[1],pdgDgCand3,pdgDgV0toDaughters,mcArray,kTRUE);
1571     if (mcLabel2!=-1) {
1572       isDp2K0Spi=1;
1573       AliAODMCParticle *lambdaCpartMC = (AliAODMCParticle*)mcArray->At(mcLabel2);
1574       if (lambdaCpartMC) {
1575         ptCandByMC = lambdaCpartMC->Pt();
1576         yCandByMC  = lambdaCpartMC->Y() ;
1577       }
1578     }
1579     if (mcLabel3!=-1) {
1580       isDs2K0SK=1;
1581       AliAODMCParticle *lambdaCpartMC = (AliAODMCParticle*)mcArray->At(mcLabel3);
1582       if (lambdaCpartMC) {
1583         ptCandByMC = lambdaCpartMC->Pt();
1584         yCandByMC  = lambdaCpartMC->Y() ;
1585       }
1586     }
1587
1588     Int_t pdgCand4 = 313; // K*(892)+ -> pi+ K0S
1589     Int_t pdgCand5 = 325; // K*(1430)+ -> pi+ K0S
1590     Int_t pdgDgCand4[2]={211,310};
1591     Int_t pdgDgCand5[2]={211,310};
1592     pdgDgV0toDaughters[0]=211;
1593     pdgDgV0toDaughters[1]=211;
1594     mcLabel4 = part->MatchToMC(pdgCand4,pdgDgCand4[1],pdgDgCand4,pdgDgV0toDaughters,mcArray,kTRUE);
1595     mcLabel5 = part->MatchToMC(pdgCand5,pdgDgCand5[1],pdgDgCand5,pdgDgV0toDaughters,mcArray,kTRUE);
1596     if (mcLabel4!=-1) {
1597       isKstar12K0Spi=1;
1598       AliAODMCParticle *lambdaCpartMC = (AliAODMCParticle*)mcArray->At(mcLabel4);
1599       if (lambdaCpartMC) {
1600         ptCandByMC = lambdaCpartMC->Pt();
1601         yCandByMC  = lambdaCpartMC->Y() ;
1602       }
1603     }
1604     if (mcLabel5!=-1) {
1605       isKstar22K0Spi=1;
1606       AliAODMCParticle *lambdaCpartMC = (AliAODMCParticle*)mcArray->At(mcLabel5);
1607       if (lambdaCpartMC) {
1608         ptCandByMC = lambdaCpartMC->Pt();
1609         yCandByMC  = lambdaCpartMC->Y() ;
1610       }
1611     }
1612   }
1613
1614   Int_t isLcByMC = isLc+isLc2LBarpi*2+isLc2Lpi*4+isDp2K0Spi*8+isDs2K0SK*16+isKstar12K0Spi*32+isKstar22K0Spi*64;
1615
1616   Bool_t isMCparticleInFiducialAcceptance = kTRUE;
1617   if (isLc || isLc2LBarpi || isLc2Lpi || isDp2K0Spi || isDs2K0SK || isKstar12K0Spi || isKstar22K0Spi) {
1618     isMCparticleInFiducialAcceptance = cutsAnal->IsInFiducialAcceptance(ptCandByMC,yCandByMC);
1619   }
1620
1621   Int_t isK0S = 0;
1622   Int_t isLambda = 0;
1623   Int_t isLambdaBar = 0;
1624   Int_t isGamma = 0;
1625   if (fUseMCInfo) {
1626     Int_t pdgDg2prong[2] = {211, 211};
1627     Int_t labelK0S = v0part->MatchToMC(310,mcArray,2,pdgDg2prong);
1628     if (labelK0S>=0) isK0S = 1;
1629
1630     pdgDg2prong[0] = 211;
1631     pdgDg2prong[1] = 2212;
1632     Int_t lambdaLabel = v0part->MatchToMC(3122,mcArray,2,pdgDg2prong);
1633     if (lambdaLabel>=0) {
1634       AliAODMCParticle *lambdaTrack = (AliAODMCParticle*)mcArray->At(lambdaLabel);
1635       if (lambdaTrack->GetPdgCode()==3122) isLambda = 1;
1636       else if (lambdaTrack->GetPdgCode()==-3122) isLambdaBar = 1;
1637     }
1638
1639     pdgDg2prong[0] = 11;
1640     pdgDg2prong[1] = 11;
1641     Int_t gammaLabel = v0part->MatchToMC(22,mcArray,2,pdgDg2prong);
1642     if (gammaLabel>=0) {
1643       AliAODMCParticle *gammaTrack = (AliAODMCParticle*)mcArray->At(gammaLabel);
1644       if (gammaTrack->GetPdgCode()==22) isGamma = 1;
1645     }
1646   }
1647
1648   Int_t isV0ByMC = isK0S+isLambdaBar*2+isLambda*4+isGamma*8;
1649
1650   Int_t isBachelorSelected = (bachelor->TestFilterMask(BIT(4)))*1 + (!(bachelor->TestFilterMask(BIT(4))))*2;
1651   isBachelorSelected += (bachelor->GetLabel()<0)*4 + (bachelor->GetLabel()>=0)*8;
1652   if ( ( !(bachelor->HasPointOnITSLayer(0)) && !(bachelor->HasPointOnITSLayer(1)) ) )
1653     isBachelorSelected += 16;
1654   else {
1655     if ( bachelor->HasPointOnITSLayer(0) && !(bachelor->HasPointOnITSLayer(1)) )
1656       isBachelorSelected += 32;
1657     else if ( !(bachelor->HasPointOnITSLayer(0)) && bachelor->HasPointOnITSLayer(1) )
1658       isBachelorSelected += 64;
1659     else
1660       isBachelorSelected += 128;
1661   }
1662
1663   AliAODTrack *v0pos = (AliAODTrack*)part->Getv0PositiveTrack();
1664   AliAODTrack *v0neg = (AliAODTrack*)part->Getv0NegativeTrack(); 
1665
1666   Int_t areV0daughtersSelected = (v0pos->TestFilterMask(BIT(4)))*1 + (!(v0pos->TestFilterMask(BIT(4))))*2;
1667   areV0daughtersSelected += (v0pos->GetLabel()<0)*4 + (v0pos->GetLabel()>=0)*8;
1668   areV0daughtersSelected += (v0neg->TestFilterMask(BIT(4)))*16 + (!(v0neg->TestFilterMask(BIT(4))))*32;
1669   areV0daughtersSelected += (v0neg->GetLabel()<0)*64 + (v0neg->GetLabel()>=0)*128;
1670
1671   Double_t nSigmaITSpr=-999.;
1672   cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,4,nSigmaITSpr);
1673   Double_t nSigmaTPCpr=-999.;
1674   cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,4,nSigmaTPCpr);
1675   Double_t nSigmaTOFpr=-999.;
1676   cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,4,nSigmaTOFpr);
1677
1678   Double_t nSigmaITSpi=-999.;
1679   cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,2,nSigmaITSpi);
1680   Double_t nSigmaTPCpi=-999.;
1681   cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,2,nSigmaTPCpi);
1682   Double_t nSigmaTOFpi=-999.;
1683   cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,2,nSigmaTOFpi);
1684
1685   Double_t nSigmaITSka=-999.;
1686   cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,3,nSigmaITSka);
1687   Double_t nSigmaTPCka=-999.;
1688   cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,3,nSigmaTPCka);
1689   Double_t nSigmaTOFka=-999.;
1690   cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,3,nSigmaTOFka);
1691
1692
1693   Int_t flagToCheckCandidate = 1*(TMath::Abs(invmassK0S-mk0sPDG)<=0.050);
1694   flagToCheckCandidate+=2*((TMath::Abs(invmassLambdaBar-mLPDG)<=0.050) && (bachelor->Charge()==-1));
1695   flagToCheckCandidate+=4*((TMath::Abs(invmassLambda-mLPDG)<=0.050) && (bachelor->Charge()==+1));
1696   flagToCheckCandidate+=8*((TMath::Abs(invmassLambdaBar-mLPDG)<=0.050) && (bachelor->Charge()==+1));
1697   flagToCheckCandidate+=16*((TMath::Abs(invmassLambda-mLPDG)<=0.050) && (bachelor->Charge()==-1));
1698
1699   /*
1700   Bool_t areCutsUsingPID = cutsAnal->GetIsUsePID();
1701   cutsAnal->SetUsePID(kFALSE);
1702   Int_t aaa = cutsAnal->IsSelected(part,AliRDHFCuts::kCandidate);
1703   if ( (aaa&AliRDHFCutsLctoV0::kLcToK0Spr)==AliRDHFCutsLctoV0::kLcToK0Spr ) {
1704     if ( aaa==AliRDHFCutsLctoV0::kLcToK0Spr ) {
1705       flagToCheckCandidate = aaa; // Lc->K0S+p OK
1706     } else {
1707       if ( (aaa&AliRDHFCutsLctoV0::kLcToLpi)==AliRDHFCutsLctoV0::kLcToLpi ) {
1708         if (bachelor->Charge()==+1)
1709           flagToCheckCandidate = aaa; // Lc->Lambda+pi+
1710         else if (bachelor->Charge()==-1)
1711           flagToCheckCandidate =-aaa; // Lambda+pi- AS Lc->K0S+p candidate
1712       }
1713       if ( (aaa&AliRDHFCutsLctoV0::kLcToLBarpi)==AliRDHFCutsLctoV0::kLcToLBarpi ) {
1714         if (bachelor->Charge()==-1)
1715           flagToCheckCandidate = aaa; // Lc->LambdaBar+pi-
1716         else if (bachelor->Charge()==+1)
1717           flagToCheckCandidate =-aaa; // LambdaBar+pi+ AS Lc->K0S+p candidate
1718       }
1719     }
1720   } else {
1721     //if ( aaa==AliRDHFCutsLctoV0::kLcToK0Spr ) {
1722     //flagToCheckCandidate = -10-(AliRDHFCutsLctoV0::kLcToK0Spr); // NEVER
1723     //} else {
1724       if ( (aaa&AliRDHFCutsLctoV0::kLcToLpi)==AliRDHFCutsLctoV0::kLcToLpi ) {
1725         if (bachelor->Charge()==+1)
1726           flagToCheckCandidate = aaa; // Lc->Lambda+pi+ OK
1727         else if (bachelor->Charge()==-1)
1728           flagToCheckCandidate =-aaa; // Lambda+pi- AS Lc->Lambda+pi+ candidate
1729       }
1730       if ( (aaa&AliRDHFCutsLctoV0::kLcToLBarpi)==AliRDHFCutsLctoV0::kLcToLBarpi ) {
1731         if (bachelor->Charge()==-1)
1732           flagToCheckCandidate = aaa; // Lc->LambdaBar+pi- OK
1733         else if (bachelor->Charge()==+1)
1734           flagToCheckCandidate =-aaa; // LambdaBar+pi+ AS Lc->LambdaBar+pi- candidate
1735       }
1736       //}
1737   }
1738   cutsAnal->SetUsePID(areCutsUsingPID);
1739   */
1740
1741   fCandidateVariables[ 0] = fUseMCInfo+isLcByMC; // 0: real data; 1: bkg; 2: Lc->K0S+p; 3: Lc->LambdaBar+pbar; 5: Lc->Lambda+p
1742   fCandidateVariables[ 1] = fUseMCInfo+isV0ByMC; // 0: real data; 1: bkg; 2: K0S->pi+pi; 3: LambdaBar->pbar+pi+; 5: Lambda->p+pi-
1743   fCandidateVariables[ 2] = isBachelorSelected;
1744   fCandidateVariables[ 3] = areV0daughtersSelected;
1745   fCandidateVariables[ 4] = flagToCheckCandidate;
1746   fCandidateVariables[ 5] = invmassLc;
1747   fCandidateVariables[ 6] = invmassLc2Lpi;
1748   fCandidateVariables[ 7] = part->InvMass2Prongs(0,1,211,310); // D+ -> pi+ K0S
1749   fCandidateVariables[ 8] = part->InvMass2Prongs(0,1,321,310); // D+S -> K+ K0S
1750   fCandidateVariables[ 9] = invmassK0S;
1751   fCandidateVariables[10] = invmassLambda;
1752   fCandidateVariables[11] = invmassLambdaBar;
1753   fCandidateVariables[12] = v0part->InvMass2Prongs(0,1,11,11);
1754   fCandidateVariables[13] = part->GetDCA();
1755   fCandidateVariables[14] = dcaV0ptp;
1756   fCandidateVariables[15] = part->Getd0Prong(0);
1757   fCandidateVariables[16] = part->Getd0Prong(1);
1758   fCandidateVariables[17] = v0part->Getd0Prong(0);
1759   fCandidateVariables[18] = v0part->Getd0Prong(1);
1760   fCandidateVariables[19] = part->CosPointingAngle();
1761   fCandidateVariables[20] = part->CosV0PointingAngle();
1762   fCandidateVariables[21] = v0part->RadiusSecVtx();
1763   fCandidateVariables[22] = nSigmaITSpr;
1764   fCandidateVariables[23] = nSigmaITSpi;
1765   fCandidateVariables[24] = nSigmaITSka;
1766   fCandidateVariables[25] = nSigmaTPCpr;
1767   fCandidateVariables[26] = nSigmaTPCpi;
1768   fCandidateVariables[27] = nSigmaTPCka;
1769   fCandidateVariables[28] = nSigmaTOFpr;
1770   fCandidateVariables[29] = nSigmaTOFpi;
1771   fCandidateVariables[30] = nSigmaTOFka;
1772   fCandidateVariables[31] = part->Y(4122);
1773   fCandidateVariables[32] = bachelor->Eta();
1774   fCandidateVariables[33] = v0pos->Eta();
1775   fCandidateVariables[34] = v0neg->Eta();
1776   fCandidateVariables[35] = part->P();
1777   fCandidateVariables[36] = part->Pt();
1778   fCandidateVariables[37] = v0part->P();
1779   fCandidateVariables[38] = v0part->Pt();
1780   fCandidateVariables[39] = bachelor->P();
1781   fCandidateVariables[40] = bachelor->Pt();
1782   fCandidateVariables[41] = v0pos->P();
1783   fCandidateVariables[42] = v0pos->Pt();
1784   fCandidateVariables[43] = v0neg->P();
1785   fCandidateVariables[44] = v0neg->Pt();
1786   fCandidateVariables[45] = part->DecayLength();
1787   fCandidateVariables[46] = part->DecayLengthV0();
1788   fCandidateVariables[47] = part->CosPointingAngleXY();
1789   fCandidateVariables[48] = part->CosV0PointingAngleXY();
1790   fCandidateVariables[49] = part->DecayLengthXY();
1791   fCandidateVariables[50] = part->DecayLengthXYV0();
1792   fCandidateVariables[51] = part->NormalizedDecayLength();
1793   fCandidateVariables[52] = part->NormalizedV0DecayLength();
1794   fCandidateVariables[53] = part->NormalizedDecayLengthXY();
1795   fCandidateVariables[54] = part->NormalizedV0DecayLengthXY();
1796   Double_t xVtxLc=0, yVtxLc=0, zVtxLc=0;
1797   Double_t xLcMC=0,yLcMC=0,zLcMC=0;
1798   Double_t pxVtxBachelor=0, pyVtxBachelor=0, pzVtxBachelor=0;
1799   Double_t dcaForLc = PropagateToDCA(v0part,bachelor,fBzkG, xVtxLc, yVtxLc, zVtxLc, pxVtxBachelor, pyVtxBachelor, pzVtxBachelor);
1800   if (isLc) {
1801     Int_t pdgCand0 = 4122;
1802     Int_t pdgDgLctoV0bachelor0[2]={2212,310};
1803     Int_t pdgDgV0toDaughters0[2]={211,211};
1804     Int_t mcLabel0 = part->MatchToMC(pdgCand0,pdgDgLctoV0bachelor0[1],pdgDgLctoV0bachelor0,pdgDgV0toDaughters0,mcArray,kTRUE);
1805     AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel0));
1806     if(partLc){
1807       AliAODMCParticle *partLcDaug0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(partLc->GetDaughter(0)));
1808       if(partLcDaug0){
1809         xLcMC=partLcDaug0->Xv(), yLcMC=partLcDaug0->Yv(), zLcMC=partLcDaug0->Zv();
1810       }
1811     }
1812   } else if (isLc2LBarpi || isLc2Lpi) {
1813     AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1814     AliAODMCParticle *partLcDaug0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(partLc->GetDaughter(0)));
1815     xLcMC=partLcDaug0->Xv(), yLcMC=partLcDaug0->Yv(), zLcMC=partLcDaug0->Zv();
1816   } else if (isDp2K0Spi) {
1817     AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel2));
1818     AliAODMCParticle *partLcDaug0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(partLc->GetDaughter(0)));
1819     xLcMC=partLcDaug0->Xv(), yLcMC=partLcDaug0->Yv(), zLcMC=partLcDaug0->Zv();
1820   } else if (isDs2K0SK) {
1821     AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel3));
1822     AliAODMCParticle *partLcDaug0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(partLc->GetDaughter(0)));
1823     xLcMC=partLcDaug0->Xv(), yLcMC=partLcDaug0->Yv(), zLcMC=partLcDaug0->Zv();
1824   } else if (isKstar12K0Spi) {
1825     AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel4));
1826     AliAODMCParticle *partLcDaug0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(partLc->GetDaughter(0)));
1827     xLcMC=partLcDaug0->Xv(), yLcMC=partLcDaug0->Yv(), zLcMC=partLcDaug0->Zv();
1828   } else if (isKstar22K0Spi) {
1829     AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel5));
1830     AliAODMCParticle *partLcDaug0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(partLc->GetDaughter(0)));
1831     xLcMC=partLcDaug0->Xv(), yLcMC=partLcDaug0->Yv(), zLcMC=partLcDaug0->Zv();
1832   }
1833   fCandidateVariables[55]=dcaForLc;
1834   fCandidateVariables[56]=part->GetSecVtxX();
1835   fCandidateVariables[57]=part->GetSecVtxY();
1836   fCandidateVariables[58]=part->GetSecVtxZ();
1837   fCandidateVariables[59]=xVtxLc;
1838   fCandidateVariables[60]=yVtxLc;
1839   fCandidateVariables[61]=zVtxLc;
1840   fCandidateVariables[62]=xLcMC;
1841   fCandidateVariables[63]=yLcMC;
1842   fCandidateVariables[64]=zLcMC;
1843   fCandidateVariables[65]=bachelor->Px();
1844   fCandidateVariables[66]=bachelor->Py();
1845   fCandidateVariables[67]=pxVtxBachelor;
1846   fCandidateVariables[68]=pyVtxBachelor;
1847   fCandidateVariables[69]=v0part->Px();
1848   fCandidateVariables[70]=v0part->Py();
1849   fCandidateVariables[71]=fVtx1->GetX();
1850   fCandidateVariables[72]=fVtx1->GetY();
1851   fCandidateVariables[73]=fVtx1->GetZ();
1852   fCandidateVariables[74]=part->CosThetaStar(0,4122,2212,310);
1853   fCandidateVariables[75]=part->CosThetaStar(1,4122,2212,310);
1854   fCandidateVariables[76]=v0part->Eta();
1855   fCandidateVariables[77]=v0part->Y(310);
1856   fCandidateVariables[78]=pzVtxBachelor;
1857   fCandidateVariables[79]=v0part->Pz();
1858   fCandidateVariables[80]=bachelor->Charge();
1859   fCandidateVariables[81]=isMCparticleInFiducialAcceptance;
1860   if (fUseMCInfo) {
1861     fCandidateVariables[82]=0;
1862     if (bachelor->GetLabel()!=-1) {
1863       AliAODMCParticle *partBachelor = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(bachelor->GetLabel())));
1864       if(partBachelor) fCandidateVariables[82]=partBachelor->GetPdgCode();
1865     }
1866   } else {
1867     fCandidateVariables[82]=-1;
1868   }
1869   fCandidateVariables[83] = part->InvMass2Prongs(0,1,211,310); // Kstar( 892)+ -> pi+K0S
1870   fCandidateVariables[84] = part->InvMass2Prongs(0,1,321,310); // Kstar(1430)+ -> pi+K0S
1871
1872   fCandidateVariables[85]=0;
1873   fCandidateVariables[86]=0;
1874   fCandidateVariables[87]=0;
1875   if (fUseMCInfo) {
1876     if (bachelor->GetLabel()!=-1 &&
1877         v0pos->GetLabel()!=-1 &&
1878         v0neg->GetLabel()!=-1) {
1879       const Int_t ndg=3;
1880       Int_t dgLabels[ndg]={TMath::Abs(bachelor->GetLabel()),
1881                            TMath::Abs(v0pos->GetLabel()),
1882                            TMath::Abs(v0neg->GetLabel())};
1883       Int_t ndgCk=0;
1884       Int_t *pdgDg=0;
1885       Int_t absLabelMother=-1;
1886       fCandidateVariables[85]=SearchForCommonMother(mcArray,
1887                                                     dgLabels,ndg,ndgCk,pdgDg,absLabelMother);
1888       AliAODMCParticle *part1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(v0pos->GetLabel())));
1889       AliAODMCParticle *part2 = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(v0neg->GetLabel())));
1890       fCandidateVariables[86]=part1->GetPdgCode();
1891       fCandidateVariables[87]=part2->GetPdgCode();
1892     }
1893   }
1894
1895   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1896   AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
1897   AliPIDResponse *pidResponse=inputHandler->GetPIDResponse();
1898   fCandidateVariables[88]=pidResponse->GetTOFResponse().GetStartTimeMask(bachelor->P());
1899
1900   //fCandidateVariables[65] = bachelor->Px();
1901   //fCandidateVariables[66] = bachelor->Py();
1902   //fCandidateVariables[67] = bachelor->Pz();
1903   //fCandidateVariables[68] = v0pos->Px();
1904   //fCandidateVariables[69] = v0pos->Py();
1905   //fCandidateVariables[70] = v0pos->Pz();
1906   //fCandidateVariables[71] = v0neg->Px();
1907   //fCandidateVariables[72] = v0neg->Py();
1908   //fCandidateVariables[73] = v0neg->Pz();
1909   //fCandidateVariables[74] = part->PxProng(0);
1910   //fCandidateVariables[75] = part->PyProng(0);
1911   //fCandidateVariables[76] = part->PzProng(0);
1912   //fCandidateVariables[77] = part->PxProng(1);
1913   //fCandidateVariables[78] = part->PyProng(1);
1914   //fCandidateVariables[79] = part->PzProng(1);
1915   //fCandidateVariables[80] = v0part->PxProng(0);
1916   //fCandidateVariables[81] = v0part->PyProng(0);
1917   //fCandidateVariables[82] = v0part->PzProng(0);
1918   //fCandidateVariables[83] = v0part->PxProng(1);
1919   //fCandidateVariables[84] = v0part->PyProng(1);
1920   //fCandidateVariables[85] = v0part->PzProng(1);
1921   //fCandidateVariables[86] = part->QtProng(0);
1922   //fCandidateVariables[87] = part->Alpha();
1923
1924   fVariablesTree->Fill();
1925
1926   return;
1927 }
1928
1929 //-------------------------------------------------------------------------------
1930 void AliAnalysisTaskSELc2V0bachelor::DefineTreeVariables() {
1931   //
1932   // This is to define tree variables
1933   //
1934
1935   const char* nameoutput = GetOutputSlot(4)->GetContainer()->GetName();
1936   fVariablesTree = new TTree(nameoutput,"Candidates variables tree");
1937   Int_t nVar = 89;
1938   fCandidateVariables = new Float_t [nVar];
1939   TString * fCandidateVariableNames = new TString[nVar];
1940   fCandidateVariableNames[ 0]="isLcByMC";
1941   fCandidateVariableNames[ 1]="isV0ByMC";
1942   fCandidateVariableNames[ 2]="flagToCheckBachelor";
1943   fCandidateVariableNames[ 3]="flagToCheckV0daughters";
1944   fCandidateVariableNames[ 4]="flagToCheckCandidate";
1945   fCandidateVariableNames[ 5]="massLc2K0Sp";
1946   fCandidateVariableNames[ 6]="massLc2Lambdapi";
1947   fCandidateVariableNames[ 7]="massD2K0Spi"; // D+ -> pi+ K0S
1948   fCandidateVariableNames[ 8]="massDS2K0SK"; // D+S -> K+ K0S
1949   fCandidateVariableNames[ 9]="massK0S";
1950   fCandidateVariableNames[10]="massLambda";
1951   fCandidateVariableNames[11]="massLambdaBar";
1952   fCandidateVariableNames[12]="massGamma";
1953   fCandidateVariableNames[13]="dcaLcptp"; // DCA Lc prong-to-prong
1954   fCandidateVariableNames[14]="dcaV0ptp";
1955   fCandidateVariableNames[15]="tImpParBach";
1956   fCandidateVariableNames[16]="tImpParV0";
1957   fCandidateVariableNames[17]="dcaV0postoPV";
1958   fCandidateVariableNames[18]="dcaV0negtoPV";
1959   fCandidateVariableNames[19]="cosPALc";
1960   fCandidateVariableNames[20]="cosPAK0S";
1961   fCandidateVariableNames[21]="rhoV0";
1962   fCandidateVariableNames[22]="nSigmaITSpr";
1963   fCandidateVariableNames[23]="nSigmaITSpi";
1964   fCandidateVariableNames[24]="nSigmaITSka";
1965   fCandidateVariableNames[25]="nSigmaTPCpr";
1966   fCandidateVariableNames[26]="nSigmaTPCpi";
1967   fCandidateVariableNames[27]="nSigmaTPCka";
1968   fCandidateVariableNames[28]="nSigmaTOFpr";
1969   fCandidateVariableNames[29]="nSigmaTOFpi";
1970   fCandidateVariableNames[30]="nSigmaTOFka";
1971   fCandidateVariableNames[31]="yLc";
1972   fCandidateVariableNames[32]="etaBach"; // etaBachelor
1973   fCandidateVariableNames[33]="etaV0pos"; // etaV0pos
1974   fCandidateVariableNames[34]="etaV0neg"; // etaV0neg
1975   fCandidateVariableNames[35]="LcP"; // @ DCA
1976   fCandidateVariableNames[36]="LcPt"; // @ DCA
1977   fCandidateVariableNames[37]="v0P"; // @ V0 DCA
1978   fCandidateVariableNames[38]="v0Pt"; // @ V0 DCA
1979   fCandidateVariableNames[39]="bachelorP"; // @ prim vtx
1980   fCandidateVariableNames[40]="bachelorPt"; // @ prim vtx
1981   fCandidateVariableNames[41]="V0positiveP"; // @ prim vtx
1982   fCandidateVariableNames[42]="V0positivePt"; // @ prim vtx
1983   fCandidateVariableNames[43]="V0negativeP"; // @ prim vtx
1984   fCandidateVariableNames[44]="V0negativePt"; // @ prim vtx
1985   fCandidateVariableNames[45]="decayLengthLc";
1986   fCandidateVariableNames[46]="decayLengthV0";
1987   fCandidateVariableNames[47]="cosPALcXY"; // cosPA XY x Lc
1988   fCandidateVariableNames[48]="cosPAV0XY"; // cosPA XY x V0
1989   fCandidateVariableNames[49]="decayLengthLcXY"; // decay length XY x Lc
1990   fCandidateVariableNames[50]="decayLengthV0XY"; // decay length XY x V0
1991   fCandidateVariableNames[51]="normalizedDecayLengthLc"; // normalized decay length x Lc
1992   fCandidateVariableNames[52]="normalizedDecayLengthV0"; // normalized decay length x V0
1993   fCandidateVariableNames[53]="normalizedDecayLengthXYLc"; // normalized decay length XY x Lc
1994   fCandidateVariableNames[54]="normalizedDecayLengthXYV0"; // normalized decay length XY x V0
1995   fCandidateVariableNames[55]="newLcDCA";
1996   fCandidateVariableNames[56]="xVtxLcBad";
1997   fCandidateVariableNames[57]="yVtxLcBad";
1998   fCandidateVariableNames[58]="zVtxLcBad";
1999   fCandidateVariableNames[59]="xVtxLcGood";
2000   fCandidateVariableNames[60]="yVtxLcGood";
2001   fCandidateVariableNames[61]="zVtxLcGood";
2002   fCandidateVariableNames[62]="xVtxLcMC";
2003   fCandidateVariableNames[63]="yVtxLcMC";
2004   fCandidateVariableNames[64]="zVtxLcMC";
2005   fCandidateVariableNames[65]="pxVtxBachelorBad";
2006   fCandidateVariableNames[66]="pyVtxBachelorBad";
2007   fCandidateVariableNames[67]="pxVtxBachelorGood";
2008   fCandidateVariableNames[68]="pyVtxBachelorGood";
2009   fCandidateVariableNames[69]="pxVtxV0";
2010   fCandidateVariableNames[70]="pyVtxV0";
2011   fCandidateVariableNames[71]="xPvtx";
2012   fCandidateVariableNames[72]="yPvtx";
2013   fCandidateVariableNames[73]="zPvtx";
2014   fCandidateVariableNames[74]="cosThetaStarBachelor";
2015   fCandidateVariableNames[75]="cosThetaStarV0";
2016   fCandidateVariableNames[76]="etaV0";
2017   fCandidateVariableNames[77]="yV0";
2018   fCandidateVariableNames[78]="pzVtxBachelorGood";
2019   fCandidateVariableNames[79]="pzVtxV0";
2020   fCandidateVariableNames[80]="bachelorCharge";
2021   fCandidateVariableNames[81]="isMCparticleInFiducialAcceptance";
2022   fCandidateVariableNames[82]="pdgBachelor"; // pdg MC bachelor
2023   fCandidateVariableNames[83]="massKstar12K0Spi"; // Kstar( 892)+ -> pi+ K0S
2024   fCandidateVariableNames[84]="massKstar22K0Spi"; // Kstar(1430)+ -> pi+ K0S
2025   fCandidateVariableNames[85]="pdgCandidate"; // pdg MC candidate recovered via new method
2026   fCandidateVariableNames[86]="pdgV0pos"; // pdg MC V0 positive
2027   fCandidateVariableNames[87]="pdgV0neg"; // pdg MC V0 negative
2028   fCandidateVariableNames[88]="startTimeMask"; // start time mask
2029
2030   //fCandidateVariableNames[65]="bachelorPx";
2031   //fCandidateVariableNames[66]="bachelorPy";
2032   //fCandidateVariableNames[67]="bachelorPz";
2033   //fCandidateVariableNames[68]="V0positivePx";
2034   //fCandidateVariableNames[69]="V0positivePy";
2035   //fCandidateVariableNames[70]="V0positivePz";
2036   //fCandidateVariableNames[71]="V0negativePx";
2037   //fCandidateVariableNames[72]="V0negativePy";
2038   //fCandidateVariableNames[73]="V0negativePz";
2039   //fCandidateVariableNames[74]="bachelorPxDCA";
2040   //fCandidateVariableNames[75]="bachelorPyDCA";
2041   //fCandidateVariableNames[76]="bachelorPzDCA";
2042   //fCandidateVariableNames[77]="v0PxDCA";
2043   //fCandidateVariableNames[78]="v0PyDCA";
2044   //fCandidateVariableNames[79]="v0PzDCA";
2045   //fCandidateVariableNames[80]="V0positivePxDCA";
2046   //fCandidateVariableNames[81]="V0positivePyDCA";
2047   //fCandidateVariableNames[82]="V0positivePzDCA";
2048   //fCandidateVariableNames[83]="V0negativePxDCA";
2049   //fCandidateVariableNames[84]="V0negativePyDCA";
2050   //fCandidateVariableNames[85]="V0negativePzDCA";
2051   //fCandidateVariableNames[86]="qtLc";
2052   //fCandidateVariableNames[87]="alphaLc";
2053
2054   for (Int_t ivar=0; ivar<nVar; ivar++) {
2055     fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
2056   }
2057
2058   return;
2059 }
2060
2061 //__________________________________________________________________________
2062 void  AliAnalysisTaskSELc2V0bachelor::DefineGeneralHistograms() {
2063   //
2064   // This is to define general histograms
2065   //
2066
2067   fCEvents = new TH1F("fCEvents","conter",18,0,18);
2068   fCEvents->SetStats(kTRUE);
2069   fCEvents->GetXaxis()->SetBinLabel(1,"X1");
2070   fCEvents->GetXaxis()->SetBinLabel(2,"Analyzed events");
2071   fCEvents->GetXaxis()->SetBinLabel(3,"AliAODVertex exists");
2072   fCEvents->GetXaxis()->SetBinLabel(4,"CascadesHF exists");
2073   fCEvents->GetXaxis()->SetBinLabel(5,"MCarray exists");
2074   fCEvents->GetXaxis()->SetBinLabel(6,"MCheader exists");
2075   fCEvents->GetXaxis()->SetBinLabel(7,"GetNContributors()>0");
2076   fCEvents->GetXaxis()->SetBinLabel(8,"IsEventSelected");
2077   fCEvents->GetXaxis()->SetBinLabel(9,"triggerClass!=CINT1");
2078   fCEvents->GetXaxis()->SetBinLabel(10,"triggerMask!=kAnyINT");
2079   fCEvents->GetXaxis()->SetBinLabel(11,"triggerMask!=kAny");
2080   fCEvents->GetXaxis()->SetBinLabel(12,"vtxTitle.Contains(Z)");
2081   fCEvents->GetXaxis()->SetBinLabel(13,"vtxTitle.Contains(3D)");
2082   fCEvents->GetXaxis()->SetBinLabel(14,"vtxTitle.Doesn'tContain(Z-3D)");
2083   fCEvents->GetXaxis()->SetBinLabel(15,Form("zVtx<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
2084   fCEvents->GetXaxis()->SetBinLabel(16,"!IsEventSelected");
2085   fCEvents->GetXaxis()->SetBinLabel(17,"triggerMask!=kAnyINT || triggerClass!=CINT1");
2086   fCEvents->GetXaxis()->SetBinLabel(18,Form("zVtxMC<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
2087   //fCEvents->GetXaxis()->SetTitle("");
2088   fCEvents->GetYaxis()->SetTitle("counts");
2089
2090   fOutput->Add(fCEvents);
2091   TString fillthis="";
2092
2093   if (fUseMCInfo && fAdditionalChecks) {
2094     fillthis="histMcStatLc";
2095     TH1F* mcStatisticLc = new TH1F(fillthis.Data(),"#Lambda_{c} generated and their decays",21,-10.5,10.5);
2096     fOutput->Add(mcStatisticLc);
2097   }
2098
2099   //fillthis="histopionV0SigmaVspTOF";
2100   //TH2F *hpionV0SigmaVspTOF=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2101   fillthis="histoprotonBachSigmaVspTOF";
2102   TH2F *hprotonBachSigmaVspTOF=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2103
2104   //fOutput->Add(hpionV0SigmaVspTOF);
2105   fOutput->Add(hprotonBachSigmaVspTOF);
2106
2107   //fillthis="histopionV0SigmaVspTPC";
2108   //TH2F *hpionV0SigmaVspTPC=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2109   fillthis="histoprotonBachSigmaVspTPC";
2110   TH2F *hprotonBachSigmaVspTPC=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2111
2112   //fOutput->Add(hpionV0SigmaVspTPC);
2113   fOutput->Add(hprotonBachSigmaVspTPC);
2114
2115   if (fUseMCInfo) {
2116
2117     //fillthis="histopionV0SigmaVspTOFsgn";
2118     //TH2F *hpionV0SigmaVspTOFsgn=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2119     fillthis="histoprotonBachSigmaVspTOFsgn";
2120     TH2F *hprotonBachSigmaVspTOFsgn=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2121
2122     //fOutput->Add(hpionV0SigmaVspTOFsgn);
2123     fOutput->Add(hprotonBachSigmaVspTOFsgn);
2124
2125     //fillthis="histopionV0SigmaVspTPCsgn";
2126     //TH2F *hpionV0SigmaVspTPCsgn=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2127     fillthis="histoprotonBachSigmaVspTPCsgn";
2128     TH2F *hprotonBachSigmaVspTPCsgn=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2129
2130     //fOutput->Add(hpionV0SigmaVspTPCsgn);
2131     fOutput->Add(hprotonBachSigmaVspTPCsgn);
2132
2133
2134     //fillthis="histopionV0SigmaVspTOFbkg";
2135     //TH2F *hpionV0SigmaVspTOFbkg=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2136     fillthis="histoprotonBachSigmaVspTOFbkg";
2137     TH2F *hprotonBachSigmaVspTOFbkg=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2138
2139     //fOutput->Add(hpionV0SigmaVspTOFbkg);
2140     fOutput->Add(hprotonBachSigmaVspTOFbkg);
2141
2142     //fillthis="histopionV0SigmaVspTPCbkg";
2143     //TH2F *hpionV0SigmaVspTPCbkg=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2144     fillthis="histoprotonBachSigmaVspTPCbkg";
2145     TH2F *hprotonBachSigmaVspTPCbkg=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2146
2147     //fOutput->Add(hpionV0SigmaVspTPCbkg);
2148     fOutput->Add(hprotonBachSigmaVspTPCbkg);
2149
2150   }
2151
2152   if (fAdditionalChecks) {
2153
2154     TH1F *hZ2 = new TH1F("hZ2","",100,-50.,50.);
2155     fOutput->Add(hZ2);
2156     TH1F *hZ3 = new TH1F("hZ3","",100,-50.,50.);
2157     fOutput->Add(hZ3);
2158     TH1F *hZ4 = new TH1F("hZ4","",100,-50.,50.);
2159     fOutput->Add(hZ4);
2160     TH1F *hZ5 = new TH1F("hZ5","",100,-50.,50.);
2161     fOutput->Add(hZ5);
2162     TH1F *hZ6 = new TH1F("hZ6","",100,-50.,50.);
2163     fOutput->Add(hZ6);
2164     TH1F *hZ7 = new TH1F("hZ7","",100,-50.,50.);
2165     fOutput->Add(hZ7);
2166     TH1F *hZ8 = new TH1F("hZ8","",100,-50.,50.);
2167     fOutput->Add(hZ8);
2168     TH1F *hZ9 = new TH1F("hZ9","",100,-50.,50.);
2169     fOutput->Add(hZ9);
2170     TH1F *hZ10 = new TH1F("hZ10","",100,-50.,50.);
2171     fOutput->Add(hZ10);
2172     TH1F *hZ11 = new TH1F("hZ11","",100,-50.,50.);
2173     fOutput->Add(hZ11);
2174     TH1F *hZ12 = new TH1F("hZ12","",100,-50.,50.);
2175     fOutput->Add(hZ12);
2176     TH1F *hZ13 = new TH1F("hZ13","",100,-50.,50.);
2177     fOutput->Add(hZ13);
2178     TH1F *hZ14 = new TH1F("hZ14","",100,-50.,50.);
2179     fOutput->Add(hZ14);
2180     TH1F *hZ15 = new TH1F("hZ15","",100,-50.,50.);
2181     fOutput->Add(hZ15);
2182     TH1F *hZ16 = new TH1F("hZ16","",100,-50.,50.);
2183     fOutput->Add(hZ16);
2184   }
2185
2186   TH1F *hCandidateSelection = new TH1F("hCandidateSelection","",10,-0.5,9.5);
2187   hCandidateSelection->GetXaxis()->SetBinLabel(1,"IsEventSelected");
2188   hCandidateSelection->GetXaxis()->SetBinLabel(2,"IsSecondaryVtx");
2189   hCandidateSelection->GetXaxis()->SetBinLabel(3,"V0toPosNeg");
2190   hCandidateSelection->GetXaxis()->SetBinLabel(4,"offlineV0");
2191   hCandidateSelection->GetXaxis()->SetBinLabel(5,"isInFiducialAcceptance");
2192   hCandidateSelection->GetXaxis()->SetBinLabel(6,"analCuts::kTracks");
2193   hCandidateSelection->GetXaxis()->SetBinLabel(7,"analCuts::kCandidateNoPID");
2194   hCandidateSelection->GetXaxis()->SetBinLabel(8,"analCuts::kPID");
2195   hCandidateSelection->GetXaxis()->SetBinLabel(9,"analCuts::kCandidateWithPID");
2196   hCandidateSelection->GetXaxis()->SetBinLabel(10,"analCuts::kAll");
2197   fOutput->Add(hCandidateSelection);
2198
2199   TH1F *hEventsWithCandidates = new TH1F("hEventsWithCandidates","conter",11,5.5,16.5);
2200   hEventsWithCandidates->GetXaxis()->SetBinLabel(1,"GetNContributors()>0");
2201   hEventsWithCandidates->GetXaxis()->SetBinLabel(2,"IsEventSelected");
2202   hEventsWithCandidates->GetXaxis()->SetBinLabel(3,"triggerClass!=CINT1");
2203   hEventsWithCandidates->GetXaxis()->SetBinLabel(4,"triggerMask!=kAnyINT");
2204   hEventsWithCandidates->GetXaxis()->SetBinLabel(5,"triggerMask!=kAny");
2205   hEventsWithCandidates->GetXaxis()->SetBinLabel(6,"vtxTitle.Contains(Z)");
2206   hEventsWithCandidates->GetXaxis()->SetBinLabel(7,"vtxTitle.Contains(3D)");
2207   hEventsWithCandidates->GetXaxis()->SetBinLabel(8,"vtxTitle.Doesn'tContain(Z-3D)");
2208   hEventsWithCandidates->GetXaxis()->SetBinLabel(9,Form("zVtx<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
2209   hEventsWithCandidates->GetXaxis()->SetBinLabel(10,"!IsEventSelected");
2210   hEventsWithCandidates->GetXaxis()->SetBinLabel(11,"triggerMask!=kAnyINT || triggerClass!=CINT1");
2211   fOutput->Add(hEventsWithCandidates);
2212
2213   if (fAdditionalChecks) {
2214
2215     TH1F *hZ6a = new TH1F("hZ6a","",100,-50.,50.);
2216     fOutput->Add(hZ6a);
2217     TH1F *hZ7a = new TH1F("hZ7a","",100,-50.,50.);
2218     fOutput->Add(hZ7a);
2219     TH1F *hZ8a = new TH1F("hZ8a","",100,-50.,50.);
2220     fOutput->Add(hZ8a);
2221     TH1F *hZ9a = new TH1F("hZ9a","",100,-50.,50.);
2222     fOutput->Add(hZ9a);
2223     TH1F *hZ10a = new TH1F("hZ10a","",100,-50.,50.);
2224     fOutput->Add(hZ10a);
2225     TH1F *hZ11a = new TH1F("hZ11a","",100,-50.,50.);
2226     fOutput->Add(hZ11a);
2227     TH1F *hZ12a = new TH1F("hZ12a","",100,-50.,50.);
2228     fOutput->Add(hZ12a);
2229     TH1F *hZ13a = new TH1F("hZ13a","",100,-50.,50.);
2230     fOutput->Add(hZ13a);
2231     TH1F *hZ14a = new TH1F("hZ14a","",100,-50.,50.);
2232     fOutput->Add(hZ14a);
2233     TH1F *hZ15a = new TH1F("hZ15a","",100,-50.,50.);
2234     fOutput->Add(hZ15a);
2235     TH1F *hZ16a = new TH1F("hZ16a","",100,-50.,50.);
2236     fOutput->Add(hZ16a);
2237   }
2238
2239   TH1F *hSwitchOnCandidates1 = new TH1F("hSwitchOnCandidates1","",15,-7.5,7.5);
2240   fOutput->Add(hSwitchOnCandidates1);
2241   TH1F *hSwitchOnCandidates2 = new TH1F("hSwitchOnCandidates2","",15,-7.5,7.5);
2242   fOutput->Add(hSwitchOnCandidates2);
2243   TH1F *hSwitchOnCandidates3 = new TH1F("hSwitchOnCandidates3","",15,-7.5,7.5);
2244   fOutput->Add(hSwitchOnCandidates3);
2245   TH1F *hSwitchOnCandidates4 = new TH1F("hSwitchOnCandidates4","",15,-7.5,7.5);
2246   fOutput->Add(hSwitchOnCandidates4);
2247
2248   return;
2249 }
2250
2251 //________________________________________________________________________
2252 void  AliAnalysisTaskSELc2V0bachelor::DefineAnalysisHistograms() {
2253   //
2254   // This is to define analysis histograms
2255   //
2256
2257   if (fIsK0SAnalysis) DefineK0SHistos();// hK0S histos declarations
2258
2259   return;
2260 }
2261
2262 //________________________________________________________________________
2263 void  AliAnalysisTaskSELc2V0bachelor::FillAnalysisHistograms(AliAODRecoCascadeHF *part, Bool_t isBachelorID, TString appendthis) {
2264   //
2265   // This is to fill analysis histograms
2266   //
2267
2268   TString fillthis="";
2269
2270   Double_t invmassLc = part->InvMassLctoK0sP();
2271   Double_t lambdacpt = part->Pt();
2272
2273   AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
2274   Double_t momBach  = bachelor->P();
2275
2276   AliAODv0 * v0part = (AliAODv0*)part->Getv0();
2277   Double_t momK0S = v0part->P();
2278   Double_t ptK0S = v0part->Pt();
2279   Double_t dcaV0ptp = v0part->GetDCA();
2280   Double_t invmassK0S = v0part->MassK0Short();
2281
2282     fillthis="histK0SMass"+appendthis;
2283     //    cout << fillthis << endl;
2284     ((TH2F*)(fOutputAll->FindObject(fillthis)))->Fill(invmassK0S,ptK0S);
2285     if (isBachelorID)  ((TH2F*)(fOutputPIDBach->FindObject(fillthis)))->Fill(invmassK0S,ptK0S);
2286
2287     fillthis="histpK0Svsp"+appendthis;
2288     //    cout << fillthis << endl;
2289     ((TH2F*)(fOutputAll->FindObject(fillthis)))->Fill(momBach,momK0S);
2290     if (isBachelorID)  ((TH2F*)(fOutputPIDBach->FindObject(fillthis)))->Fill(momBach,momK0S);
2291
2292     fillthis="histDCAtoPVvspK0S"+appendthis;
2293     //    cout << fillthis << endl;
2294     ((TH2F*)(fOutputAll->FindObject(fillthis)))->Fill(momK0S,dcaV0ptp);
2295     if (isBachelorID)  ((TH2F*)(fOutputPIDBach->FindObject(fillthis)))->Fill(momK0S,dcaV0ptp);
2296
2297     fillthis="histArmPodK0S"+appendthis;
2298     //    cout << fillthis << endl;
2299     FillArmPodDistribution(v0part,fillthis,fOutputAll);
2300     if (isBachelorID) FillArmPodDistribution(v0part,fillthis,fOutputPIDBach);
2301
2302     fillthis="histLcMassByK0S"+appendthis;
2303     //    cout << fillthis << endl;
2304     ((TH2F*)(fOutputAll->FindObject(fillthis)))->Fill(invmassLc,lambdacpt);
2305     if (isBachelorID)((TH2F*)(fOutputPIDBach->FindObject(fillthis)))->Fill(invmassLc,lambdacpt);
2306
2307     return;
2308 }
2309 //---------------------------
2310 Double_t AliAnalysisTaskSELc2V0bachelor::PropagateToDCA(AliAODv0 *v, AliAODTrack *bachelor, Double_t b,
2311                                                         Double_t &xVtxLc, Double_t &yVtxLc, Double_t &zVtxLc,
2312                                                         Double_t &pxVtxBachelor, Double_t &pyVtxBachelor, Double_t &pzVtxBachelor) {
2313   //--------------------------------------------------------------------
2314   // This function returns the DCA between the V0 and the track
2315   // This is a copy of AliCascadeVertexer::PropagateToDCA(...) method
2316   //--------------------------------------------------------------------
2317
2318   // Get AliExternalTrackParam out of the AliAODTracks                                                      
2319   Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign;
2320   bachelor->PxPyPz(pxpypz);
2321   bachelor->XvYvZv(xyz);
2322   bachelor->GetCovarianceXYZPxPyPz(cv);
2323   sign=bachelor->Charge();
2324   AliExternalTrackParam *t = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
2325
2326   Double_t alpha=t->GetAlpha(), cs1=TMath::Cos(alpha), sn1=TMath::Sin(alpha);
2327   //Double_t alpha = GetAlpha(xyz,pxpypz), cs1=TMath::Cos(alpha), sn1=TMath::Sin(alpha);
2328
2329   // position and momentum of bachelor
2330   Double_t x1=xyz[0], y1=xyz[1], z1=xyz[2];
2331   Double_t px1=pxpypz[0], py1=pxpypz[1], pz1=pxpypz[2];
2332
2333   // position and momentum of V0
2334   Double_t x2=v->DecayVertexV0X(),
2335     y2=v->DecayVertexV0Y(),
2336     z2=v->DecayVertexV0Z();
2337   Double_t px2=v->Px(),
2338     py2=v->Py(),
2339     pz2=v->Pz();
2340
2341   /*
2342   AliAODTrack *trackP = (AliAODTrack*) v->GetDaughter(0);
2343   //Double_t pxpypzP[3]; trackP->PxPyPz(pxpypzP);
2344   //Double_t xyzP[3]; trackP->XvYvZv(xyzP);
2345   Double_t cvP[21]; trackP->GetCovarianceXYZPxPyPz(cvP);
2346   //Short_t signP=trackP->Charge();
2347   //AliExternalTrackParam *tP = new AliExternalTrackParam(xyzP,pxpypzP,cvP,signP);
2348
2349   // Get AliExternalTrackParam out of the AliAODTrack
2350   AliAODTrack *trackN = (AliAODTrack*) v->GetDaughter(1);
2351   //Double_t pxpypzN[3]; trackN->PxPyPz(pxpypzN);
2352   //Double_t xyzN[3]; trackN->XvYvZv(xyzN);
2353   Double_t cvN[21]; trackN->GetCovarianceXYZPxPyPz(cvN);
2354   //Short_t signN=trackN->Charge();
2355   //AliExternalTrackParam *tN = new AliExternalTrackParam(xyzN,pxpypzN,cvN,signN);
2356
2357   Double_t xyzV0[3]={x2,y2,z2};
2358   Double_t pxpypzV0[3]={px2,py2,pz2};
2359   Double_t cvV0[21]; for (Int_t ii=0; ii<21; ii++) cvV0[ii]=cvP[ii]+cvN[ii];
2360   AliNeutralTrackParam *trackV0 = new AliNeutralTrackParam(xyzV0,pxpypzV0,cvV0,0);
2361   */
2362  
2363   // calculation dca
2364   Double_t dd= Det(x2-x1,y2-y1,z2-z1,px1,py1,pz1,px2,py2,pz2);
2365   Double_t ax= Det(py1,pz1,py2,pz2);
2366   Double_t ay=-Det(px1,pz1,px2,pz2);
2367   Double_t az= Det(px1,py1,px2,py2);
2368
2369   Double_t dca=TMath::Abs(dd)/TMath::Sqrt(ax*ax + ay*ay + az*az);
2370
2371   // bachelor point @ the DCA
2372   Double_t t1 = Det(x2-x1,y2-y1,z2-z1,px2,py2,pz2,ax,ay,az)/
2373     Det(px1,py1,pz1,px2,py2,pz2,ax,ay,az);
2374   x1 += px1*t1; y1 += py1*t1; z1 += pz1*t1;
2375
2376   //propagate track to the point of DCA
2377   Double_t rho1=x1*cs1 + y1*sn1;
2378   if (!t->PropagateTo(rho1,b)) {
2379     Error("PropagateToDCA","Propagation failed !");
2380     delete t; t=NULL;
2381     return 1.e+33;
2382   }
2383
2384   Double_t pBachelorDCA[3]; t->GetPxPyPz(pBachelorDCA);
2385   pxVtxBachelor=pBachelorDCA[0], pyVtxBachelor=pBachelorDCA[1], pzVtxBachelor=pBachelorDCA[2];
2386
2387   delete t; t=NULL;
2388
2389   // V0 point @ the DCA
2390   Double_t t2 = Det(x1-x2,y1-y2,z1-z2,px1,py1,pz1,ax,ay,az)/
2391     Det(px2,py2,pz2,px1,py1,pz1,ax,ay,az);
2392   x2 += px2*t2; y2 += py2*t2; z2 += pz2*t2;
2393
2394
2395   // Lc decay vtx
2396   xVtxLc = 0.5*(x1+x2);
2397   yVtxLc = 0.5*(y1+y2);
2398   zVtxLc = 0.5*(z1+z2);
2399   
2400   return dca;
2401
2402 }
2403
2404 //---------------------------
2405 Double_t AliAnalysisTaskSELc2V0bachelor::GetAlpha(Double_t xyz[3],Double_t pxpypz[3])
2406 {
2407   //
2408   // To estimate alpha according to what done in the AliExternalTrackParam::Set(...) method
2409   //
2410
2411   Double_t alpha = 0.;
2412
2413   const double kSafe = 1e-5;
2414   Double_t radPos2 = xyz[0]*xyz[0]+xyz[1]*xyz[1];
2415   Double_t radMax  = 45.; // approximately ITS outer radius
2416   if (radPos2 < radMax*radMax) { // inside the ITS
2417     alpha = TMath::ATan2(pxpypz[1],pxpypz[0]);
2418   } else { // outside the ITS
2419     Float_t phiPos = TMath::Pi()+TMath::ATan2(-xyz[1], -xyz[0]);
2420      alpha =
2421        TMath::DegToRad()*(20*((((Int_t)(phiPos*TMath::RadToDeg()))/20))+10);
2422   }
2423
2424   Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
2425   // protection: avoid alpha being too close to 0 or +-pi/2
2426   if (TMath::Abs(sn)<2*kSafe) {
2427     if (alpha>0) alpha += alpha< TMath::Pi()/2. ?  2*kSafe : -2*kSafe;
2428     else         alpha += alpha>-TMath::Pi()/2. ? -2*kSafe :  2*kSafe;
2429     cs=TMath::Cos(alpha);
2430     sn=TMath::Sin(alpha);
2431   }
2432   else if (TMath::Abs(cs)<2*kSafe) {
2433     if (alpha>0) alpha += alpha> TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
2434     else         alpha += alpha>-TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
2435     cs=TMath::Cos(alpha);
2436     sn=TMath::Sin(alpha);
2437   }
2438
2439
2440   return alpha;
2441 }
2442
2443 //---------------------------
2444 Double_t AliAnalysisTaskSELc2V0bachelor::Det(Double_t a00, Double_t a01,
2445                                              Double_t a10, Double_t a11) const {
2446   //--------------------------------------------------------------------
2447   // This function calculates locally a 2x2 determinant.
2448   // This is a copy of the AliCascadeVertexer::Det(...) method
2449   //--------------------------------------------------------------------
2450   return a00*a11 - a01*a10;
2451 }
2452
2453 //---------------------------
2454 Double_t AliAnalysisTaskSELc2V0bachelor::Det(Double_t a00,Double_t a01,Double_t a02,
2455                                              Double_t a10,Double_t a11,Double_t a12,
2456                                              Double_t a20,Double_t a21,Double_t a22) const {
2457   //--------------------------------------------------------------------
2458   // This function calculates locally a 3x3 determinant
2459   // This is a copy of the AliCascadeVertexer::Det(...) method
2460   //--------------------------------------------------------------------
2461   return  a00*Det(a11,a12,a21,a22)-a01*Det(a10,a12,a20,a22)+a02*Det(a10,a11,a20,a21);
2462 }
2463
2464 //----------------------------------------------------------------------------
2465 Int_t AliAnalysisTaskSELc2V0bachelor::MatchToMClabelC(AliAODRecoCascadeHF *candidate,
2466                                                       TClonesArray *mcArray)
2467 {
2468   //
2469   // Check if this candidate is matched to a MC signal  Lc -> p K0S + X
2470   // If no, return -1
2471   // If yes, return label (>=0) of the AliAODMCParticle
2472   // 
2473
2474   AliAODv0 *theV0 = dynamic_cast<AliAODv0*>(candidate->Getv0()); // the V0
2475   AliVTrack *trk = dynamic_cast<AliVTrack*>(candidate->GetBachelor()); // the bachelor
2476   if (!trk || !theV0) return -1;
2477
2478   if (trk->GetLabel()==-1) return -1;
2479   Int_t bachLabels = TMath::Abs(trk->GetLabel());
2480   AliAODMCParticle*bachelorMC = dynamic_cast<AliAODMCParticle*>(mcArray->At(bachLabels));
2481   if (!bachelorMC) return -1;
2482   if (TMath::Abs(bachelorMC->GetPdgCode())!=2212) return -1;
2483   Int_t indexMotherBach = bachelorMC->GetMother();
2484   if (indexMotherBach==-1) return -1;
2485
2486   Int_t pdgDg2prong[2] = {211,211};
2487   Int_t lab2Prong = theV0->MatchToMC(310,mcArray,2,pdgDg2prong); // the V0
2488   if(lab2Prong<0) return -1;
2489   AliAODMCParticle*partK0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(lab2Prong));
2490   if (!partK0S) return -1;
2491   Int_t indexMotherK0S = partK0S->GetMother();
2492   if (indexMotherK0S==-1) return -1;
2493   AliAODMCParticle*partK0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(indexMotherK0S));
2494   if (!partK0) return -1;
2495   Int_t indexMotherK0 = partK0->GetMother();
2496   if (indexMotherK0==-1) return -1;
2497
2498   if (indexMotherBach!=indexMotherK0) return -1; // p e K0S sono fratelli
2499
2500   AliAODMCParticle*partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(indexMotherK0));
2501   if (!partLc) return -1;
2502   Int_t ndg2 = partLc->GetDaughter(1)-partLc->GetDaughter(0)+1;
2503   if (ndg2==2) return -1;
2504
2505   TString stringaCheck = Form(">>>>>>>> %d -> ",partLc->GetPdgCode());
2506   for(Int_t ii=0; ii<ndg2; ii++) {
2507     AliAODMCParticle* partDau=(AliAODMCParticle*)(mcArray->At(partLc->GetDaughter(0)+ii));
2508     stringaCheck.Append(Form("  %d",partDau->GetPdgCode()));
2509   }
2510   printf("%s \n",stringaCheck.Data());
2511
2512   return indexMotherBach;
2513
2514 }
2515 //-----------------------------------------------------------------------------
2516 Int_t AliAnalysisTaskSELc2V0bachelor::SearchForCommonMother(TClonesArray *mcArray,
2517                                                             Int_t dgLabels[10],Int_t ndg,
2518                                                             Int_t &ndgCk, Int_t *pdgDg, Int_t &absLabelMother) const
2519 {
2520   //
2521   // Check if this candidate is matched to a MC signal
2522   // If no, return 0
2523   // If yes, return pdgCode of particle
2524   // 
2525
2526   Int_t lab=-1,labMother=-1,pdgMother=0;
2527   AliAODMCParticle *part=0;
2528   AliAODMCParticle *mother=0;
2529
2530   // loop on daughter labels
2531   TArrayI **labelMother = new TArrayI*[ndg];
2532   for(Int_t i=0; i<ndg; i++) labelMother[i] = new TArrayI(0);
2533   for(Int_t i=0; i<ndg; i++) {
2534     lab = TMath::Abs(dgLabels[i]);
2535     if(lab<0) {
2536       printf("daughter with negative label %d\n",lab);
2537       return 0;
2538     }
2539     part = (AliAODMCParticle*)mcArray->At(lab);
2540     if(!part) { 
2541       printf("no MC particle\n");
2542       return 0;
2543     }
2544
2545     mother = part;
2546     while(mother->GetMother()>=0) {
2547       labMother=mother->GetMother();
2548       mother = (AliAODMCParticle*)mcArray->At(labMother);
2549       if(!mother) {
2550         printf("no MC mother particle\n");
2551         break;
2552       }
2553       pdgMother = TMath::Abs(mother->GetPdgCode());
2554       if (pdgMother<10 || (pdgMother>18 && pdgMother<111)) {
2555         break;
2556       }
2557       labelMother[i]->Set(labelMother[i]->GetSize()+1);
2558       labelMother[i]->AddAt(labMother,labelMother[i]->GetSize()-1);
2559     }
2560
2561   } // end loop on daughters
2562
2563   for(Int_t i=0; i<ndg; i++) {
2564     AliAODMCParticle*part0 = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels[i]));
2565     printf("part[%d]->GetLabel()=%d(%d) | ",i,dgLabels[i],part0->GetPdgCode());
2566     printf("labelMother[%d] = ",i);
2567     for (Int_t jj=0;jj<labelMother[i]->GetSize(); jj++)
2568       printf("%d, ",labelMother[i]->At(jj));
2569     printf("\n");
2570   }
2571
2572   Int_t pdgToBeReturned=0;
2573   Bool_t found=kFALSE;
2574   for (Int_t ii=0;ii<labelMother[0]->GetSize(); ii++) {
2575     for (Int_t jj=0;jj<labelMother[1]->GetSize(); jj++) {
2576       for (Int_t kk=0;kk<labelMother[2]->GetSize(); kk++) {
2577         if (labelMother[0]->At(ii)==labelMother[1]->At(jj) &&
2578             labelMother[1]->At(jj)==labelMother[2]->At(kk) &&
2579             labelMother[0]->At(ii)!=0 && labelMother[0]->At(ii)!=1 && !found) {
2580           mother = (AliAODMCParticle*)mcArray->At(labelMother[0]->At(ii));
2581           pdgToBeReturned=mother->GetPdgCode();
2582           absLabelMother=labelMother[0]->At(ii);
2583           printf("FOUND label for the mother of this candidate: %d (PDG=%d)\n",labelMother[0]->At(ii),pdgToBeReturned);
2584           mother->Print();
2585           found = kTRUE;
2586           ndgCk=3;
2587           pdgDg = new Int_t[ndgCk];
2588           for (Int_t aa=0; aa<ndgCk; aa++) {
2589             AliAODMCParticle *partMC = (AliAODMCParticle*)mcArray->At(dgLabels[aa]);
2590             pdgDg[aa]=partMC->GetPdgCode();
2591           }
2592           break;
2593         }
2594       }
2595     }
2596   }
2597
2598   delete [] labelMother;
2599
2600   return pdgToBeReturned;
2601
2602 }