Updates in LC->Kos+proton code (Annalisa)
[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   Double_t ptCandByMC = 0.;//fmcPartCandidate->Pt();
1532   Double_t yCandByMC  = 0.;//fmcPartCandidate->Y() ;
1533   if (fUseMCInfo) {
1534     if (isLc) {
1535       Int_t pdgCand0 = 4122;
1536       Int_t pdgDgLctoV0bachelor0[2]={2212,310};
1537       Int_t pdgDgV0toDaughters0[2]={211,211};
1538       Int_t mcLabelLc2pK0S = part->MatchToMC(pdgCand0,pdgDgLctoV0bachelor0[1],pdgDgLctoV0bachelor0,pdgDgV0toDaughters0,mcArray,kTRUE);
1539       AliAODMCParticle *lambdaCpartMC = (AliAODMCParticle*)mcArray->At(mcLabelLc2pK0S);
1540       if (lambdaCpartMC) {
1541         ptCandByMC = lambdaCpartMC->Pt();
1542         yCandByMC  = lambdaCpartMC->Y() ;
1543       }
1544     }
1545
1546     Int_t pdgCand = 4122;
1547     Int_t pdgDgLctoV0bachelor[2]={211,3122};
1548     Int_t pdgDgV0toDaughters[2]={2212,211};
1549     mcLabel = part->MatchToMC(pdgCand,pdgDgLctoV0bachelor[1],pdgDgLctoV0bachelor,pdgDgV0toDaughters,mcArray,kTRUE);
1550     if (mcLabel!=-1) {
1551       if (bachelor->Charge()==-1) isLc2LBarpi=1;
1552       if (bachelor->Charge()==+1) isLc2Lpi=1;
1553       AliAODMCParticle *lambdaCpartMC = (AliAODMCParticle*)mcArray->At(mcLabel);
1554       if (lambdaCpartMC) {
1555         ptCandByMC = lambdaCpartMC->Pt();
1556         yCandByMC  = lambdaCpartMC->Y() ;
1557       }
1558     }
1559
1560     Int_t pdgCand2 = 411; // D+ -> pi+ K0S
1561     Int_t pdgCand3 = 431; // Ds+ -> K+ K0S
1562     Int_t pdgDgCand2[2]={211,310};
1563     Int_t pdgDgCand3[2]={321,310};
1564     pdgDgV0toDaughters[0]=211;
1565     pdgDgV0toDaughters[1]=211;
1566     mcLabel2 = part->MatchToMC(pdgCand2,pdgDgCand2[1],pdgDgCand2,pdgDgV0toDaughters,mcArray,kTRUE);
1567     mcLabel3 = part->MatchToMC(pdgCand3,pdgDgCand3[1],pdgDgCand3,pdgDgV0toDaughters,mcArray,kTRUE);
1568     if (mcLabel2!=-1) {
1569       isDp2K0Spi=1;
1570       AliAODMCParticle *lambdaCpartMC = (AliAODMCParticle*)mcArray->At(mcLabel2);
1571       if (lambdaCpartMC) {
1572         ptCandByMC = lambdaCpartMC->Pt();
1573         yCandByMC  = lambdaCpartMC->Y() ;
1574       }
1575     }
1576     if (mcLabel3!=-1) {
1577       isDs2K0SK=1;
1578       AliAODMCParticle *lambdaCpartMC = (AliAODMCParticle*)mcArray->At(mcLabel3);
1579       if (lambdaCpartMC) {
1580         ptCandByMC = lambdaCpartMC->Pt();
1581         yCandByMC  = lambdaCpartMC->Y() ;
1582       }
1583     }
1584   }
1585
1586   Int_t isLcByMC = isLc+isLc2LBarpi*2+isLc2Lpi*4+isDp2K0Spi*8+isDs2K0SK*16;
1587
1588   Bool_t isMCparticleInFiducialAcceptance = kTRUE;
1589   if (isLc || isLc2LBarpi || isLc2Lpi || isDp2K0Spi || isDs2K0SK) {
1590     isMCparticleInFiducialAcceptance = cutsAnal->IsInFiducialAcceptance(ptCandByMC,yCandByMC);
1591   }
1592
1593   Int_t isK0S = 0;
1594   Int_t isLambda = 0;
1595   Int_t isLambdaBar = 0;
1596   Int_t isGamma = 0;
1597   if (fUseMCInfo) {
1598     Int_t pdgDg2prong[2] = {211, 211};
1599     Int_t labelK0S = v0part->MatchToMC(310,mcArray,2,pdgDg2prong);
1600     if (labelK0S>=0) isK0S = 1;
1601
1602     pdgDg2prong[0] = 211;
1603     pdgDg2prong[1] = 2212;
1604     Int_t lambdaLabel = v0part->MatchToMC(3122,mcArray,2,pdgDg2prong);
1605     if (lambdaLabel>=0) {
1606       AliAODMCParticle *lambdaTrack = (AliAODMCParticle*)mcArray->At(lambdaLabel);
1607       if (lambdaTrack->GetPdgCode()==3122) isLambda = 1;
1608       else if (lambdaTrack->GetPdgCode()==-3122) isLambdaBar = 1;
1609     }
1610
1611     pdgDg2prong[0] = 11;
1612     pdgDg2prong[1] = 11;
1613     Int_t gammaLabel = v0part->MatchToMC(22,mcArray,2,pdgDg2prong);
1614     if (gammaLabel>=0) {
1615       AliAODMCParticle *gammaTrack = (AliAODMCParticle*)mcArray->At(gammaLabel);
1616       if (gammaTrack->GetPdgCode()==22) isGamma = 1;
1617     }
1618   }
1619
1620   Int_t isV0ByMC = isK0S+isLambdaBar*2+isLambda*4+isGamma*8;
1621
1622   Int_t isBachelorSelected = (bachelor->TestFilterMask(BIT(4)))*1 + (!(bachelor->TestFilterMask(BIT(4))))*2;
1623   isBachelorSelected += (bachelor->GetLabel()<0)*4 + (bachelor->GetLabel()>=0)*8;
1624   if ( ( !(bachelor->HasPointOnITSLayer(0)) && !(bachelor->HasPointOnITSLayer(1)) ) )
1625     isBachelorSelected += 16;
1626   else {
1627     if ( bachelor->HasPointOnITSLayer(0) && !(bachelor->HasPointOnITSLayer(1)) )
1628       isBachelorSelected += 32;
1629     else if ( !(bachelor->HasPointOnITSLayer(0)) && bachelor->HasPointOnITSLayer(1) )
1630       isBachelorSelected += 64;
1631     else
1632       isBachelorSelected += 128;
1633   }
1634
1635   AliAODTrack *v0pos = (AliAODTrack*)part->Getv0PositiveTrack();
1636   AliAODTrack *v0neg = (AliAODTrack*)part->Getv0NegativeTrack(); 
1637
1638   Int_t areV0daughtersSelected = (v0pos->TestFilterMask(BIT(4)))*1 + (!(v0pos->TestFilterMask(BIT(4))))*2;
1639   areV0daughtersSelected += (v0pos->GetLabel()<0)*4 + (v0pos->GetLabel()>=0)*8;
1640   areV0daughtersSelected += (v0neg->TestFilterMask(BIT(4)))*16 + (!(v0neg->TestFilterMask(BIT(4))))*32;
1641   areV0daughtersSelected += (v0neg->GetLabel()<0)*64 + (v0neg->GetLabel()>=0)*128;
1642
1643   Double_t nSigmaITSpr=-999.;
1644   cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,4,nSigmaITSpr);
1645   Double_t nSigmaTPCpr=-999.;
1646   cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,4,nSigmaTPCpr);
1647   Double_t nSigmaTOFpr=-999.;
1648   cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,4,nSigmaTOFpr);
1649
1650   Double_t nSigmaITSpi=-999.;
1651   cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,2,nSigmaITSpi);
1652   Double_t nSigmaTPCpi=-999.;
1653   cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,2,nSigmaTPCpi);
1654   Double_t nSigmaTOFpi=-999.;
1655   cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,2,nSigmaTOFpi);
1656
1657   Double_t nSigmaITSka=-999.;
1658   cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,3,nSigmaITSka);
1659   Double_t nSigmaTPCka=-999.;
1660   cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,3,nSigmaTPCka);
1661   Double_t nSigmaTOFka=-999.;
1662   cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,3,nSigmaTOFka);
1663
1664
1665   Int_t flagToCheckCandidate = 1*(TMath::Abs(invmassK0S-mk0sPDG)<=0.050);
1666   flagToCheckCandidate+=2*((TMath::Abs(invmassLambdaBar-mLPDG)<=0.050) && (bachelor->Charge()==-1));
1667   flagToCheckCandidate+=4*((TMath::Abs(invmassLambda-mLPDG)<=0.050) && (bachelor->Charge()==+1));
1668   flagToCheckCandidate+=8*((TMath::Abs(invmassLambdaBar-mLPDG)<=0.050) && (bachelor->Charge()==+1));
1669   flagToCheckCandidate+=16*((TMath::Abs(invmassLambda-mLPDG)<=0.050) && (bachelor->Charge()==-1));
1670
1671   /*
1672   Bool_t areCutsUsingPID = cutsAnal->GetIsUsePID();
1673   cutsAnal->SetUsePID(kFALSE);
1674   Int_t aaa = cutsAnal->IsSelected(part,AliRDHFCuts::kCandidate);
1675   if ( (aaa&AliRDHFCutsLctoV0::kLcToK0Spr)==AliRDHFCutsLctoV0::kLcToK0Spr ) {
1676     if ( aaa==AliRDHFCutsLctoV0::kLcToK0Spr ) {
1677       flagToCheckCandidate = aaa; // Lc->K0S+p OK
1678     } else {
1679       if ( (aaa&AliRDHFCutsLctoV0::kLcToLpi)==AliRDHFCutsLctoV0::kLcToLpi ) {
1680         if (bachelor->Charge()==+1)
1681           flagToCheckCandidate = aaa; // Lc->Lambda+pi+
1682         else if (bachelor->Charge()==-1)
1683           flagToCheckCandidate =-aaa; // Lambda+pi- AS Lc->K0S+p candidate
1684       }
1685       if ( (aaa&AliRDHFCutsLctoV0::kLcToLBarpi)==AliRDHFCutsLctoV0::kLcToLBarpi ) {
1686         if (bachelor->Charge()==-1)
1687           flagToCheckCandidate = aaa; // Lc->LambdaBar+pi-
1688         else if (bachelor->Charge()==+1)
1689           flagToCheckCandidate =-aaa; // LambdaBar+pi+ AS Lc->K0S+p candidate
1690       }
1691     }
1692   } else {
1693     //if ( aaa==AliRDHFCutsLctoV0::kLcToK0Spr ) {
1694     //flagToCheckCandidate = -10-(AliRDHFCutsLctoV0::kLcToK0Spr); // NEVER
1695     //} else {
1696       if ( (aaa&AliRDHFCutsLctoV0::kLcToLpi)==AliRDHFCutsLctoV0::kLcToLpi ) {
1697         if (bachelor->Charge()==+1)
1698           flagToCheckCandidate = aaa; // Lc->Lambda+pi+ OK
1699         else if (bachelor->Charge()==-1)
1700           flagToCheckCandidate =-aaa; // Lambda+pi- AS Lc->Lambda+pi+ candidate
1701       }
1702       if ( (aaa&AliRDHFCutsLctoV0::kLcToLBarpi)==AliRDHFCutsLctoV0::kLcToLBarpi ) {
1703         if (bachelor->Charge()==-1)
1704           flagToCheckCandidate = aaa; // Lc->LambdaBar+pi- OK
1705         else if (bachelor->Charge()==+1)
1706           flagToCheckCandidate =-aaa; // LambdaBar+pi+ AS Lc->LambdaBar+pi- candidate
1707       }
1708       //}
1709   }
1710   cutsAnal->SetUsePID(areCutsUsingPID);
1711   */
1712
1713   fCandidateVariables[ 0] = fUseMCInfo+isLcByMC; // 0: real data; 1: bkg; 2: Lc->K0S+p; 3: Lc->LambdaBar+pbar; 5: Lc->Lambda+p
1714   fCandidateVariables[ 1] = fUseMCInfo+isV0ByMC; // 0: real data; 1: bkg; 2: K0S->pi+pi; 3: LambdaBar->pbar+pi+; 5: Lambda->p+pi-
1715   fCandidateVariables[ 2] = isBachelorSelected;
1716   fCandidateVariables[ 3] = areV0daughtersSelected;
1717   fCandidateVariables[ 4] = flagToCheckCandidate;
1718   fCandidateVariables[ 5] = invmassLc;
1719   fCandidateVariables[ 6] = invmassLc2Lpi;
1720   fCandidateVariables[ 7] = part->InvMass2Prongs(0,1,211,310); // D+ -> pi+ K0S
1721   fCandidateVariables[ 8] = part->InvMass2Prongs(0,1,321,310); // D+S -> K+ K0S
1722   fCandidateVariables[ 9] = invmassK0S;
1723   fCandidateVariables[10] = invmassLambda;
1724   fCandidateVariables[11] = invmassLambdaBar;
1725   fCandidateVariables[12] = v0part->InvMass2Prongs(0,1,11,11);
1726   fCandidateVariables[13] = part->GetDCA();
1727   fCandidateVariables[14] = dcaV0ptp;
1728   fCandidateVariables[15] = part->Getd0Prong(0);
1729   fCandidateVariables[16] = part->Getd0Prong(1);
1730   fCandidateVariables[17] = v0part->Getd0Prong(0);
1731   fCandidateVariables[18] = v0part->Getd0Prong(1);
1732   fCandidateVariables[19] = part->CosPointingAngle();
1733   fCandidateVariables[20] = part->CosV0PointingAngle();
1734   fCandidateVariables[21] = v0part->RadiusSecVtx();
1735   fCandidateVariables[22] = nSigmaITSpr;
1736   fCandidateVariables[23] = nSigmaITSpi;
1737   fCandidateVariables[24] = nSigmaITSka;
1738   fCandidateVariables[25] = nSigmaTPCpr;
1739   fCandidateVariables[26] = nSigmaTPCpi;
1740   fCandidateVariables[27] = nSigmaTPCka;
1741   fCandidateVariables[28] = nSigmaTOFpr;
1742   fCandidateVariables[29] = nSigmaTOFpi;
1743   fCandidateVariables[30] = nSigmaTOFka;
1744   fCandidateVariables[31] = part->Y(4122);
1745   fCandidateVariables[32] = bachelor->Eta();
1746   fCandidateVariables[33] = v0pos->Eta();
1747   fCandidateVariables[34] = v0neg->Eta();
1748   fCandidateVariables[35] = part->P();
1749   fCandidateVariables[36] = part->Pt();
1750   fCandidateVariables[37] = v0part->P();
1751   fCandidateVariables[38] = v0part->Pt();
1752   fCandidateVariables[39] = bachelor->P();
1753   fCandidateVariables[40] = bachelor->Pt();
1754   fCandidateVariables[41] = v0pos->P();
1755   fCandidateVariables[42] = v0pos->Pt();
1756   fCandidateVariables[43] = v0neg->P();
1757   fCandidateVariables[44] = v0neg->Pt();
1758   fCandidateVariables[45] = part->DecayLength();
1759   fCandidateVariables[46] = part->DecayLengthV0();
1760   fCandidateVariables[47] = part->CosPointingAngleXY();
1761   fCandidateVariables[48] = part->CosV0PointingAngleXY();
1762   fCandidateVariables[49] = part->DecayLengthXY();
1763   fCandidateVariables[50] = part->DecayLengthXYV0();
1764   fCandidateVariables[51] = part->NormalizedDecayLength();
1765   fCandidateVariables[52] = part->NormalizedV0DecayLength();
1766   fCandidateVariables[53] = part->NormalizedDecayLengthXY();
1767   fCandidateVariables[54] = part->NormalizedV0DecayLengthXY();
1768   Double_t xVtxLc=0, yVtxLc=0, zVtxLc=0;
1769   Double_t xLcMC=0,yLcMC=0,zLcMC=0;
1770   Double_t pxVtxBachelor=0, pyVtxBachelor=0, pzVtxBachelor=0;
1771   Double_t dcaForLc = PropagateToDCA(v0part,bachelor,fBzkG, xVtxLc, yVtxLc, zVtxLc, pxVtxBachelor, pyVtxBachelor, pzVtxBachelor);
1772   if (isLc) {
1773     Int_t pdgCand0 = 4122;
1774     Int_t pdgDgLctoV0bachelor0[2]={2212,310};
1775     Int_t pdgDgV0toDaughters0[2]={211,211};
1776     Int_t mcLabel0 = part->MatchToMC(pdgCand0,pdgDgLctoV0bachelor0[1],pdgDgLctoV0bachelor0,pdgDgV0toDaughters0,mcArray,kTRUE);
1777     AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel0));
1778     if(partLc){
1779       AliAODMCParticle *partLcDaug0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(partLc->GetDaughter(0)));
1780       if(partLcDaug0){
1781         xLcMC=partLcDaug0->Xv(), yLcMC=partLcDaug0->Yv(), zLcMC=partLcDaug0->Zv();
1782       }
1783     }
1784   } else if (isLc2LBarpi || isLc2Lpi) {
1785     AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1786     AliAODMCParticle *partLcDaug0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(partLc->GetDaughter(0)));
1787     xLcMC=partLcDaug0->Xv(), yLcMC=partLcDaug0->Yv(), zLcMC=partLcDaug0->Zv();
1788   } else if (isDp2K0Spi) {
1789     AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel2));
1790     AliAODMCParticle *partLcDaug0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(partLc->GetDaughter(0)));
1791     xLcMC=partLcDaug0->Xv(), yLcMC=partLcDaug0->Yv(), zLcMC=partLcDaug0->Zv();
1792   } else if (isDs2K0SK) {
1793     AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel3));
1794     AliAODMCParticle *partLcDaug0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(partLc->GetDaughter(0)));
1795     xLcMC=partLcDaug0->Xv(), yLcMC=partLcDaug0->Yv(), zLcMC=partLcDaug0->Zv();
1796   }
1797   fCandidateVariables[55]=dcaForLc;
1798   fCandidateVariables[56]=part->GetSecVtxX();
1799   fCandidateVariables[57]=part->GetSecVtxY();
1800   fCandidateVariables[58]=part->GetSecVtxZ();
1801   fCandidateVariables[59]=xVtxLc;
1802   fCandidateVariables[60]=yVtxLc;
1803   fCandidateVariables[61]=zVtxLc;
1804   fCandidateVariables[62]=xLcMC;
1805   fCandidateVariables[63]=yLcMC;
1806   fCandidateVariables[64]=zLcMC;
1807   fCandidateVariables[65]=bachelor->Px();
1808   fCandidateVariables[66]=bachelor->Py();
1809   fCandidateVariables[67]=pxVtxBachelor;
1810   fCandidateVariables[68]=pyVtxBachelor;
1811   fCandidateVariables[69]=v0part->Px();
1812   fCandidateVariables[70]=v0part->Py();
1813   fCandidateVariables[71]=fVtx1->GetX();
1814   fCandidateVariables[72]=fVtx1->GetY();
1815   fCandidateVariables[73]=fVtx1->GetZ();
1816   fCandidateVariables[74]=part->CosThetaStar(0,4122,2212,310);
1817   fCandidateVariables[75]=part->CosThetaStar(1,4122,2212,310);
1818   fCandidateVariables[76]=v0part->Eta();
1819   fCandidateVariables[77]=v0part->Y(310);
1820   fCandidateVariables[78]=pzVtxBachelor;
1821   fCandidateVariables[79]=v0part->Pz();
1822   fCandidateVariables[80]=bachelor->Charge();
1823   fCandidateVariables[81]=isMCparticleInFiducialAcceptance;
1824   if (fUseMCInfo) {
1825     if (bachelor->GetLabel()!=-1) {
1826       AliAODMCParticle *partBachelor = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(bachelor->GetLabel())));
1827       fCandidateVariables[82]=partBachelor->GetPdgCode();
1828     } else {
1829       fCandidateVariables[82]=0;
1830     }
1831   } else {
1832     fCandidateVariables[82]=-1;
1833   }
1834
1835   //fCandidateVariables[65] = bachelor->Px();
1836   //fCandidateVariables[66] = bachelor->Py();
1837   //fCandidateVariables[67] = bachelor->Pz();
1838   //fCandidateVariables[68] = v0pos->Px();
1839   //fCandidateVariables[69] = v0pos->Py();
1840   //fCandidateVariables[70] = v0pos->Pz();
1841   //fCandidateVariables[71] = v0neg->Px();
1842   //fCandidateVariables[72] = v0neg->Py();
1843   //fCandidateVariables[73] = v0neg->Pz();
1844   //fCandidateVariables[74] = part->PxProng(0);
1845   //fCandidateVariables[75] = part->PyProng(0);
1846   //fCandidateVariables[76] = part->PzProng(0);
1847   //fCandidateVariables[77] = part->PxProng(1);
1848   //fCandidateVariables[78] = part->PyProng(1);
1849   //fCandidateVariables[79] = part->PzProng(1);
1850   //fCandidateVariables[80] = v0part->PxProng(0);
1851   //fCandidateVariables[81] = v0part->PyProng(0);
1852   //fCandidateVariables[82] = v0part->PzProng(0);
1853   //fCandidateVariables[83] = v0part->PxProng(1);
1854   //fCandidateVariables[84] = v0part->PyProng(1);
1855   //fCandidateVariables[85] = v0part->PzProng(1);
1856   //fCandidateVariables[86] = part->QtProng(0);
1857   //fCandidateVariables[87] = part->Alpha();
1858
1859   fVariablesTree->Fill();
1860
1861   return;
1862 }
1863
1864 //-------------------------------------------------------------------------------
1865 void AliAnalysisTaskSELc2V0bachelor::DefineTreeVariables() {
1866   //
1867   // This is to define tree variables
1868   //
1869
1870   const char* nameoutput = GetOutputSlot(4)->GetContainer()->GetName();
1871   fVariablesTree = new TTree(nameoutput,"Candidates variables tree");
1872   Int_t nVar = 83;
1873   fCandidateVariables = new Float_t [nVar];
1874   TString * fCandidateVariableNames = new TString[nVar];
1875   fCandidateVariableNames[ 0]="isLcByMC";
1876   fCandidateVariableNames[ 1]="isV0ByMC";
1877   fCandidateVariableNames[ 2]="flagToCheckBachelor";
1878   fCandidateVariableNames[ 3]="flagToCheckV0daughters";
1879   fCandidateVariableNames[ 4]="flagToCheckCandidate";
1880   fCandidateVariableNames[ 5]="massLc2K0Sp";
1881   fCandidateVariableNames[ 6]="massLc2Lambdapi";
1882   fCandidateVariableNames[ 7]="massD2K0Spi"; // D+ -> pi+ K0S
1883   fCandidateVariableNames[ 8]="massDS2K0SK"; // D+S -> K+ K0S
1884   fCandidateVariableNames[ 9]="massK0S";
1885   fCandidateVariableNames[10]="massLambda";
1886   fCandidateVariableNames[11]="massLambdaBar";
1887   fCandidateVariableNames[12]="massGamma";
1888   fCandidateVariableNames[13]="dcaLcptp"; // DCA Lc prong-to-prong
1889   fCandidateVariableNames[14]="dcaV0ptp";
1890   fCandidateVariableNames[15]="tImpParBach";
1891   fCandidateVariableNames[16]="tImpParV0";
1892   fCandidateVariableNames[17]="dcaV0postoPV";
1893   fCandidateVariableNames[18]="dcaV0negtoPV";
1894   fCandidateVariableNames[19]="cosPALc";
1895   fCandidateVariableNames[20]="cosPAK0S";
1896   fCandidateVariableNames[21]="rhoV0";
1897   fCandidateVariableNames[22]="nSigmaITSpr";
1898   fCandidateVariableNames[23]="nSigmaITSpi";
1899   fCandidateVariableNames[24]="nSigmaITSka";
1900   fCandidateVariableNames[25]="nSigmaTPCpr";
1901   fCandidateVariableNames[26]="nSigmaTPCpi";
1902   fCandidateVariableNames[27]="nSigmaTPCka";
1903   fCandidateVariableNames[28]="nSigmaTOFpr";
1904   fCandidateVariableNames[29]="nSigmaTOFpi";
1905   fCandidateVariableNames[30]="nSigmaTOFka";
1906   fCandidateVariableNames[31]="yLc";
1907   fCandidateVariableNames[32]="etaBach"; // etaBachelor
1908   fCandidateVariableNames[33]="etaV0pos"; // etaV0pos
1909   fCandidateVariableNames[34]="etaV0neg"; // etaV0neg
1910   fCandidateVariableNames[35]="LcP"; // @ DCA
1911   fCandidateVariableNames[36]="LcPt"; // @ DCA
1912   fCandidateVariableNames[37]="v0P"; // @ V0 DCA
1913   fCandidateVariableNames[38]="v0Pt"; // @ V0 DCA
1914   fCandidateVariableNames[39]="bachelorP"; // @ prim vtx
1915   fCandidateVariableNames[40]="bachelorPt"; // @ prim vtx
1916   fCandidateVariableNames[41]="V0positiveP"; // @ prim vtx
1917   fCandidateVariableNames[42]="V0positivePt"; // @ prim vtx
1918   fCandidateVariableNames[43]="V0negativeP"; // @ prim vtx
1919   fCandidateVariableNames[44]="V0negativePt"; // @ prim vtx
1920   fCandidateVariableNames[45]="decayLengthLc";
1921   fCandidateVariableNames[46]="decayLengthV0";
1922   fCandidateVariableNames[47]="cosPALcXY"; // cosPA XY x Lc
1923   fCandidateVariableNames[48]="cosPAV0XY"; // cosPA XY x V0
1924   fCandidateVariableNames[49]="decayLengthLcXY"; // decay length XY x Lc
1925   fCandidateVariableNames[50]="decayLengthV0XY"; // decay length XY x V0
1926   fCandidateVariableNames[51]="normalizedDecayLengthLc"; // normalized decay length x Lc
1927   fCandidateVariableNames[52]="normalizedDecayLengthV0"; // normalized decay length x V0
1928   fCandidateVariableNames[53]="normalizedDecayLengthXYLc"; // normalized decay length XY x Lc
1929   fCandidateVariableNames[54]="normalizedDecayLengthXYV0"; // normalized decay length XY x V0
1930   fCandidateVariableNames[55]="newLcDCA";
1931   fCandidateVariableNames[56]="xVtxLcBad";
1932   fCandidateVariableNames[57]="yVtxLcBad";
1933   fCandidateVariableNames[58]="zVtxLcBad";
1934   fCandidateVariableNames[59]="xVtxLcGood";
1935   fCandidateVariableNames[60]="yVtxLcGood";
1936   fCandidateVariableNames[61]="zVtxLcGood";
1937   fCandidateVariableNames[62]="xVtxLcMC";
1938   fCandidateVariableNames[63]="yVtxLcMC";
1939   fCandidateVariableNames[64]="zVtxLcMC";
1940   fCandidateVariableNames[65]="pxVtxBachelorBad";
1941   fCandidateVariableNames[66]="pyVtxBachelorBad";
1942   fCandidateVariableNames[67]="pxVtxBachelorGood";
1943   fCandidateVariableNames[68]="pyVtxBachelorGood";
1944   fCandidateVariableNames[69]="pxVtxV0";
1945   fCandidateVariableNames[70]="pyVtxV0";
1946   fCandidateVariableNames[71]="xPvtx";
1947   fCandidateVariableNames[72]="yPvtx";
1948   fCandidateVariableNames[73]="zPvtx";
1949   fCandidateVariableNames[74]="cosThetaStarBachelor";
1950   fCandidateVariableNames[75]="cosThetaStarV0";
1951   fCandidateVariableNames[76]="etaV0";
1952   fCandidateVariableNames[77]="yV0";
1953   fCandidateVariableNames[78]="pzVtxBachelorGood";
1954   fCandidateVariableNames[79]="pzVtxV0";
1955   fCandidateVariableNames[80]="bachelorCharge";
1956   fCandidateVariableNames[81]="isMCparticleInFiducialAcceptance";
1957   fCandidateVariableNames[82]="pdgCodeBachelor";
1958
1959   //fCandidateVariableNames[65]="bachelorPx";
1960   //fCandidateVariableNames[66]="bachelorPy";
1961   //fCandidateVariableNames[67]="bachelorPz";
1962   //fCandidateVariableNames[68]="V0positivePx";
1963   //fCandidateVariableNames[69]="V0positivePy";
1964   //fCandidateVariableNames[70]="V0positivePz";
1965   //fCandidateVariableNames[71]="V0negativePx";
1966   //fCandidateVariableNames[72]="V0negativePy";
1967   //fCandidateVariableNames[73]="V0negativePz";
1968   //fCandidateVariableNames[74]="bachelorPxDCA";
1969   //fCandidateVariableNames[75]="bachelorPyDCA";
1970   //fCandidateVariableNames[76]="bachelorPzDCA";
1971   //fCandidateVariableNames[77]="v0PxDCA";
1972   //fCandidateVariableNames[78]="v0PyDCA";
1973   //fCandidateVariableNames[79]="v0PzDCA";
1974   //fCandidateVariableNames[80]="V0positivePxDCA";
1975   //fCandidateVariableNames[81]="V0positivePyDCA";
1976   //fCandidateVariableNames[82]="V0positivePzDCA";
1977   //fCandidateVariableNames[83]="V0negativePxDCA";
1978   //fCandidateVariableNames[84]="V0negativePyDCA";
1979   //fCandidateVariableNames[85]="V0negativePzDCA";
1980   //fCandidateVariableNames[86]="qtLc";
1981   //fCandidateVariableNames[87]="alphaLc";
1982
1983   for (Int_t ivar=0; ivar<nVar; ivar++) {
1984     fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
1985   }
1986
1987   return;
1988 }
1989
1990 //__________________________________________________________________________
1991 void  AliAnalysisTaskSELc2V0bachelor::DefineGeneralHistograms() {
1992   //
1993   // This is to define general histograms
1994   //
1995
1996   fCEvents = new TH1F("fCEvents","conter",18,0,18);
1997   fCEvents->SetStats(kTRUE);
1998   fCEvents->GetXaxis()->SetBinLabel(1,"X1");
1999   fCEvents->GetXaxis()->SetBinLabel(2,"Analyzed events");
2000   fCEvents->GetXaxis()->SetBinLabel(3,"AliAODVertex exists");
2001   fCEvents->GetXaxis()->SetBinLabel(4,"CascadesHF exists");
2002   fCEvents->GetXaxis()->SetBinLabel(5,"MCarray exists");
2003   fCEvents->GetXaxis()->SetBinLabel(6,"MCheader exists");
2004   fCEvents->GetXaxis()->SetBinLabel(7,"GetNContributors()>0");
2005   fCEvents->GetXaxis()->SetBinLabel(8,"IsEventSelected");
2006   fCEvents->GetXaxis()->SetBinLabel(9,"triggerClass!=CINT1");
2007   fCEvents->GetXaxis()->SetBinLabel(10,"triggerMask!=kAnyINT");
2008   fCEvents->GetXaxis()->SetBinLabel(11,"triggerMask!=kAny");
2009   fCEvents->GetXaxis()->SetBinLabel(12,"vtxTitle.Contains(Z)");
2010   fCEvents->GetXaxis()->SetBinLabel(13,"vtxTitle.Contains(3D)");
2011   fCEvents->GetXaxis()->SetBinLabel(14,"vtxTitle.Doesn'tContain(Z-3D)");
2012   fCEvents->GetXaxis()->SetBinLabel(15,Form("zVtx<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
2013   fCEvents->GetXaxis()->SetBinLabel(16,"!IsEventSelected");
2014   fCEvents->GetXaxis()->SetBinLabel(17,"triggerMask!=kAnyINT || triggerClass!=CINT1");
2015   fCEvents->GetXaxis()->SetBinLabel(18,Form("zVtxMC<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
2016   //fCEvents->GetXaxis()->SetTitle("");
2017   fCEvents->GetYaxis()->SetTitle("counts");
2018
2019   fOutput->Add(fCEvents);
2020   TString fillthis="";
2021
2022   if (fUseMCInfo && fAdditionalChecks) {
2023     fillthis="histMcStatLc";
2024     TH1F* mcStatisticLc = new TH1F(fillthis.Data(),"#Lambda_{c} generated and their decays",21,-10.5,10.5);
2025     fOutput->Add(mcStatisticLc);
2026   }
2027
2028   //fillthis="histopionV0SigmaVspTOF";
2029   //TH2F *hpionV0SigmaVspTOF=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2030   fillthis="histoprotonBachSigmaVspTOF";
2031   TH2F *hprotonBachSigmaVspTOF=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2032
2033   //fOutput->Add(hpionV0SigmaVspTOF);
2034   fOutput->Add(hprotonBachSigmaVspTOF);
2035
2036   //fillthis="histopionV0SigmaVspTPC";
2037   //TH2F *hpionV0SigmaVspTPC=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2038   fillthis="histoprotonBachSigmaVspTPC";
2039   TH2F *hprotonBachSigmaVspTPC=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2040
2041   //fOutput->Add(hpionV0SigmaVspTPC);
2042   fOutput->Add(hprotonBachSigmaVspTPC);
2043
2044   if (fUseMCInfo) {
2045
2046     //fillthis="histopionV0SigmaVspTOFsgn";
2047     //TH2F *hpionV0SigmaVspTOFsgn=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2048     fillthis="histoprotonBachSigmaVspTOFsgn";
2049     TH2F *hprotonBachSigmaVspTOFsgn=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2050
2051     //fOutput->Add(hpionV0SigmaVspTOFsgn);
2052     fOutput->Add(hprotonBachSigmaVspTOFsgn);
2053
2054     //fillthis="histopionV0SigmaVspTPCsgn";
2055     //TH2F *hpionV0SigmaVspTPCsgn=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2056     fillthis="histoprotonBachSigmaVspTPCsgn";
2057     TH2F *hprotonBachSigmaVspTPCsgn=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2058
2059     //fOutput->Add(hpionV0SigmaVspTPCsgn);
2060     fOutput->Add(hprotonBachSigmaVspTPCsgn);
2061
2062
2063     //fillthis="histopionV0SigmaVspTOFbkg";
2064     //TH2F *hpionV0SigmaVspTOFbkg=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2065     fillthis="histoprotonBachSigmaVspTOFbkg";
2066     TH2F *hprotonBachSigmaVspTOFbkg=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2067
2068     //fOutput->Add(hpionV0SigmaVspTOFbkg);
2069     fOutput->Add(hprotonBachSigmaVspTOFbkg);
2070
2071     //fillthis="histopionV0SigmaVspTPCbkg";
2072     //TH2F *hpionV0SigmaVspTPCbkg=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2073     fillthis="histoprotonBachSigmaVspTPCbkg";
2074     TH2F *hprotonBachSigmaVspTPCbkg=new TH2F(fillthis.Data(),fillthis.Data(),300,0.,30.,100,-5.,5.);
2075
2076     //fOutput->Add(hpionV0SigmaVspTPCbkg);
2077     fOutput->Add(hprotonBachSigmaVspTPCbkg);
2078
2079   }
2080
2081   if (fAdditionalChecks) {
2082
2083     TH1F *hZ2 = new TH1F("hZ2","",100,-50.,50.);
2084     fOutput->Add(hZ2);
2085     TH1F *hZ3 = new TH1F("hZ3","",100,-50.,50.);
2086     fOutput->Add(hZ3);
2087     TH1F *hZ4 = new TH1F("hZ4","",100,-50.,50.);
2088     fOutput->Add(hZ4);
2089     TH1F *hZ5 = new TH1F("hZ5","",100,-50.,50.);
2090     fOutput->Add(hZ5);
2091     TH1F *hZ6 = new TH1F("hZ6","",100,-50.,50.);
2092     fOutput->Add(hZ6);
2093     TH1F *hZ7 = new TH1F("hZ7","",100,-50.,50.);
2094     fOutput->Add(hZ7);
2095     TH1F *hZ8 = new TH1F("hZ8","",100,-50.,50.);
2096     fOutput->Add(hZ8);
2097     TH1F *hZ9 = new TH1F("hZ9","",100,-50.,50.);
2098     fOutput->Add(hZ9);
2099     TH1F *hZ10 = new TH1F("hZ10","",100,-50.,50.);
2100     fOutput->Add(hZ10);
2101     TH1F *hZ11 = new TH1F("hZ11","",100,-50.,50.);
2102     fOutput->Add(hZ11);
2103     TH1F *hZ12 = new TH1F("hZ12","",100,-50.,50.);
2104     fOutput->Add(hZ12);
2105     TH1F *hZ13 = new TH1F("hZ13","",100,-50.,50.);
2106     fOutput->Add(hZ13);
2107     TH1F *hZ14 = new TH1F("hZ14","",100,-50.,50.);
2108     fOutput->Add(hZ14);
2109     TH1F *hZ15 = new TH1F("hZ15","",100,-50.,50.);
2110     fOutput->Add(hZ15);
2111     TH1F *hZ16 = new TH1F("hZ16","",100,-50.,50.);
2112     fOutput->Add(hZ16);
2113   }
2114
2115   TH1F *hCandidateSelection = new TH1F("hCandidateSelection","",10,-0.5,9.5);
2116   hCandidateSelection->GetXaxis()->SetBinLabel(1,"IsEventSelected");
2117   hCandidateSelection->GetXaxis()->SetBinLabel(2,"IsSecondaryVtx");
2118   hCandidateSelection->GetXaxis()->SetBinLabel(3,"V0toPosNeg");
2119   hCandidateSelection->GetXaxis()->SetBinLabel(4,"offlineV0");
2120   hCandidateSelection->GetXaxis()->SetBinLabel(5,"isInFiducialAcceptance");
2121   hCandidateSelection->GetXaxis()->SetBinLabel(6,"analCuts::kTracks");
2122   hCandidateSelection->GetXaxis()->SetBinLabel(7,"analCuts::kCandidateNoPID");
2123   hCandidateSelection->GetXaxis()->SetBinLabel(8,"analCuts::kPID");
2124   hCandidateSelection->GetXaxis()->SetBinLabel(9,"analCuts::kCandidateWithPID");
2125   hCandidateSelection->GetXaxis()->SetBinLabel(10,"analCuts::kAll");
2126   fOutput->Add(hCandidateSelection);
2127
2128   TH1F *hEventsWithCandidates = new TH1F("hEventsWithCandidates","conter",11,5.5,16.5);
2129   hEventsWithCandidates->GetXaxis()->SetBinLabel(1,"GetNContributors()>0");
2130   hEventsWithCandidates->GetXaxis()->SetBinLabel(2,"IsEventSelected");
2131   hEventsWithCandidates->GetXaxis()->SetBinLabel(3,"triggerClass!=CINT1");
2132   hEventsWithCandidates->GetXaxis()->SetBinLabel(4,"triggerMask!=kAnyINT");
2133   hEventsWithCandidates->GetXaxis()->SetBinLabel(5,"triggerMask!=kAny");
2134   hEventsWithCandidates->GetXaxis()->SetBinLabel(6,"vtxTitle.Contains(Z)");
2135   hEventsWithCandidates->GetXaxis()->SetBinLabel(7,"vtxTitle.Contains(3D)");
2136   hEventsWithCandidates->GetXaxis()->SetBinLabel(8,"vtxTitle.Doesn'tContain(Z-3D)");
2137   hEventsWithCandidates->GetXaxis()->SetBinLabel(9,Form("zVtx<=%2.0fcm",fAnalCuts->GetMaxVtxZ()));
2138   hEventsWithCandidates->GetXaxis()->SetBinLabel(10,"!IsEventSelected");
2139   hEventsWithCandidates->GetXaxis()->SetBinLabel(11,"triggerMask!=kAnyINT || triggerClass!=CINT1");
2140   fOutput->Add(hEventsWithCandidates);
2141
2142   if (fAdditionalChecks) {
2143
2144     TH1F *hZ6a = new TH1F("hZ6a","",100,-50.,50.);
2145     fOutput->Add(hZ6a);
2146     TH1F *hZ7a = new TH1F("hZ7a","",100,-50.,50.);
2147     fOutput->Add(hZ7a);
2148     TH1F *hZ8a = new TH1F("hZ8a","",100,-50.,50.);
2149     fOutput->Add(hZ8a);
2150     TH1F *hZ9a = new TH1F("hZ9a","",100,-50.,50.);
2151     fOutput->Add(hZ9a);
2152     TH1F *hZ10a = new TH1F("hZ10a","",100,-50.,50.);
2153     fOutput->Add(hZ10a);
2154     TH1F *hZ11a = new TH1F("hZ11a","",100,-50.,50.);
2155     fOutput->Add(hZ11a);
2156     TH1F *hZ12a = new TH1F("hZ12a","",100,-50.,50.);
2157     fOutput->Add(hZ12a);
2158     TH1F *hZ13a = new TH1F("hZ13a","",100,-50.,50.);
2159     fOutput->Add(hZ13a);
2160     TH1F *hZ14a = new TH1F("hZ14a","",100,-50.,50.);
2161     fOutput->Add(hZ14a);
2162     TH1F *hZ15a = new TH1F("hZ15a","",100,-50.,50.);
2163     fOutput->Add(hZ15a);
2164     TH1F *hZ16a = new TH1F("hZ16a","",100,-50.,50.);
2165     fOutput->Add(hZ16a);
2166   }
2167
2168   TH1F *hSwitchOnCandidates1 = new TH1F("hSwitchOnCandidates1","",15,-7.5,7.5);
2169   fOutput->Add(hSwitchOnCandidates1);
2170   TH1F *hSwitchOnCandidates2 = new TH1F("hSwitchOnCandidates2","",15,-7.5,7.5);
2171   fOutput->Add(hSwitchOnCandidates2);
2172   TH1F *hSwitchOnCandidates3 = new TH1F("hSwitchOnCandidates3","",15,-7.5,7.5);
2173   fOutput->Add(hSwitchOnCandidates3);
2174   TH1F *hSwitchOnCandidates4 = new TH1F("hSwitchOnCandidates4","",15,-7.5,7.5);
2175   fOutput->Add(hSwitchOnCandidates4);
2176
2177   return;
2178 }
2179
2180 //________________________________________________________________________
2181 void  AliAnalysisTaskSELc2V0bachelor::DefineAnalysisHistograms() {
2182   //
2183   // This is to define analysis histograms
2184   //
2185
2186   if (fIsK0SAnalysis) DefineK0SHistos();// hK0S histos declarations
2187
2188   return;
2189 }
2190
2191 //________________________________________________________________________
2192 void  AliAnalysisTaskSELc2V0bachelor::FillAnalysisHistograms(AliAODRecoCascadeHF *part, Bool_t isBachelorID, TString appendthis) {
2193   //
2194   // This is to fill analysis histograms
2195   //
2196
2197   TString fillthis="";
2198
2199   Double_t invmassLc = part->InvMassLctoK0sP();
2200   Double_t lambdacpt = part->Pt();
2201
2202   AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
2203   Double_t momBach  = bachelor->P();
2204
2205   AliAODv0 * v0part = (AliAODv0*)part->Getv0();
2206   Double_t momK0S = v0part->P();
2207   Double_t ptK0S = v0part->Pt();
2208   Double_t dcaV0ptp = v0part->GetDCA();
2209   Double_t invmassK0S = v0part->MassK0Short();
2210
2211     fillthis="histK0SMass"+appendthis;
2212     cout << fillthis << endl;
2213     ((TH2F*)(fOutputAll->FindObject(fillthis)))->Fill(invmassK0S,ptK0S);
2214     if (isBachelorID)  ((TH2F*)(fOutputPIDBach->FindObject(fillthis)))->Fill(invmassK0S,ptK0S);
2215
2216     fillthis="histpK0Svsp"+appendthis;
2217     cout << fillthis << endl;
2218     ((TH2F*)(fOutputAll->FindObject(fillthis)))->Fill(momBach,momK0S);
2219     if (isBachelorID)  ((TH2F*)(fOutputPIDBach->FindObject(fillthis)))->Fill(momBach,momK0S);
2220
2221     fillthis="histDCAtoPVvspK0S"+appendthis;
2222     cout << fillthis << endl;
2223     ((TH2F*)(fOutputAll->FindObject(fillthis)))->Fill(momK0S,dcaV0ptp);
2224     if (isBachelorID)  ((TH2F*)(fOutputPIDBach->FindObject(fillthis)))->Fill(momK0S,dcaV0ptp);
2225
2226     fillthis="histArmPodK0S"+appendthis;
2227     cout << fillthis << endl;
2228     FillArmPodDistribution(v0part,fillthis,fOutputAll);
2229     if (isBachelorID) FillArmPodDistribution(v0part,fillthis,fOutputPIDBach);
2230
2231     fillthis="histLcMassByK0S"+appendthis;
2232     cout << fillthis << endl;
2233     ((TH2F*)(fOutputAll->FindObject(fillthis)))->Fill(invmassLc,lambdacpt);
2234     if (isBachelorID)((TH2F*)(fOutputPIDBach->FindObject(fillthis)))->Fill(invmassLc,lambdacpt);
2235
2236     return;
2237 }
2238 //---------------------------
2239 Double_t AliAnalysisTaskSELc2V0bachelor::PropagateToDCA(AliAODv0 *v, AliAODTrack *bachelor, Double_t b,
2240                                                         Double_t &xVtxLc, Double_t &yVtxLc, Double_t &zVtxLc,
2241                                                         Double_t &pxVtxBachelor, Double_t &pyVtxBachelor, Double_t &pzVtxBachelor) {
2242   //--------------------------------------------------------------------
2243   // This function returns the DCA between the V0 and the track
2244   // This is a copy of AliCascadeVertexer::PropagateToDCA(...) method
2245   //--------------------------------------------------------------------
2246
2247   // Get AliExternalTrackParam out of the AliAODTracks                                                      
2248   Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign;
2249   bachelor->PxPyPz(pxpypz);
2250   bachelor->XvYvZv(xyz);
2251   bachelor->GetCovarianceXYZPxPyPz(cv);
2252   sign=bachelor->Charge();
2253   AliExternalTrackParam *t = new AliExternalTrackParam(xyz,pxpypz,cv,sign);
2254
2255   Double_t alpha=t->GetAlpha(), cs1=TMath::Cos(alpha), sn1=TMath::Sin(alpha);
2256   //Double_t alpha = GetAlpha(xyz,pxpypz), cs1=TMath::Cos(alpha), sn1=TMath::Sin(alpha);
2257
2258   // position and momentum of bachelor
2259   Double_t x1=xyz[0], y1=xyz[1], z1=xyz[2];
2260   Double_t px1=pxpypz[0], py1=pxpypz[1], pz1=pxpypz[2];
2261
2262   // position and momentum of V0
2263   Double_t x2=v->DecayVertexV0X(),
2264     y2=v->DecayVertexV0Y(),
2265     z2=v->DecayVertexV0Z();
2266   Double_t px2=v->Px(),
2267     py2=v->Py(),
2268     pz2=v->Pz();
2269
2270   /*
2271   AliAODTrack *trackP = (AliAODTrack*) v->GetDaughter(0);
2272   //Double_t pxpypzP[3]; trackP->PxPyPz(pxpypzP);
2273   //Double_t xyzP[3]; trackP->XvYvZv(xyzP);
2274   Double_t cvP[21]; trackP->GetCovarianceXYZPxPyPz(cvP);
2275   //Short_t signP=trackP->Charge();
2276   //AliExternalTrackParam *tP = new AliExternalTrackParam(xyzP,pxpypzP,cvP,signP);
2277
2278   // Get AliExternalTrackParam out of the AliAODTrack
2279   AliAODTrack *trackN = (AliAODTrack*) v->GetDaughter(1);
2280   //Double_t pxpypzN[3]; trackN->PxPyPz(pxpypzN);
2281   //Double_t xyzN[3]; trackN->XvYvZv(xyzN);
2282   Double_t cvN[21]; trackN->GetCovarianceXYZPxPyPz(cvN);
2283   //Short_t signN=trackN->Charge();
2284   //AliExternalTrackParam *tN = new AliExternalTrackParam(xyzN,pxpypzN,cvN,signN);
2285
2286   Double_t xyzV0[3]={x2,y2,z2};
2287   Double_t pxpypzV0[3]={px2,py2,pz2};
2288   Double_t cvV0[21]; for (Int_t ii=0; ii<21; ii++) cvV0[ii]=cvP[ii]+cvN[ii];
2289   AliNeutralTrackParam *trackV0 = new AliNeutralTrackParam(xyzV0,pxpypzV0,cvV0,0);
2290   */
2291  
2292   // calculation dca
2293   Double_t dd= Det(x2-x1,y2-y1,z2-z1,px1,py1,pz1,px2,py2,pz2);
2294   Double_t ax= Det(py1,pz1,py2,pz2);
2295   Double_t ay=-Det(px1,pz1,px2,pz2);
2296   Double_t az= Det(px1,py1,px2,py2);
2297
2298   Double_t dca=TMath::Abs(dd)/TMath::Sqrt(ax*ax + ay*ay + az*az);
2299
2300   // bachelor point @ the DCA
2301   Double_t t1 = Det(x2-x1,y2-y1,z2-z1,px2,py2,pz2,ax,ay,az)/
2302     Det(px1,py1,pz1,px2,py2,pz2,ax,ay,az);
2303   x1 += px1*t1; y1 += py1*t1; z1 += pz1*t1;
2304
2305   //propagate track to the point of DCA
2306   Double_t rho1=x1*cs1 + y1*sn1;
2307   if (!t->PropagateTo(rho1,b)) {
2308     Error("PropagateToDCA","Propagation failed !");
2309     delete t; t=NULL;
2310     return 1.e+33;
2311   }
2312
2313   Double_t pBachelorDCA[3]; t->GetPxPyPz(pBachelorDCA);
2314   pxVtxBachelor=pBachelorDCA[0], pyVtxBachelor=pBachelorDCA[1], pzVtxBachelor=pBachelorDCA[2];
2315
2316   delete t; t=NULL;
2317
2318   // V0 point @ the DCA
2319   Double_t t2 = Det(x1-x2,y1-y2,z1-z2,px1,py1,pz1,ax,ay,az)/
2320     Det(px2,py2,pz2,px1,py1,pz1,ax,ay,az);
2321   x2 += px2*t2; y2 += py2*t2; z2 += pz2*t2;
2322
2323
2324   // Lc decay vtx
2325   xVtxLc = 0.5*(x1+x2);
2326   yVtxLc = 0.5*(y1+y2);
2327   zVtxLc = 0.5*(z1+z2);
2328   
2329   return dca;
2330
2331 }
2332
2333 //---------------------------
2334 Double_t AliAnalysisTaskSELc2V0bachelor::GetAlpha(Double_t xyz[3],Double_t pxpypz[3])
2335 {
2336   //
2337   // To estimate alpha according to what done in the AliExternalTrackParam::Set(...) method
2338   //
2339
2340   Double_t alpha = 0.;
2341
2342   const double kSafe = 1e-5;
2343   Double_t radPos2 = xyz[0]*xyz[0]+xyz[1]*xyz[1];
2344   Double_t radMax  = 45.; // approximately ITS outer radius
2345   if (radPos2 < radMax*radMax) { // inside the ITS
2346     alpha = TMath::ATan2(pxpypz[1],pxpypz[0]);
2347   } else { // outside the ITS
2348     Float_t phiPos = TMath::Pi()+TMath::ATan2(-xyz[1], -xyz[0]);
2349      alpha =
2350        TMath::DegToRad()*(20*((((Int_t)(phiPos*TMath::RadToDeg()))/20))+10);
2351   }
2352
2353   Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
2354   // protection: avoid alpha being too close to 0 or +-pi/2
2355   if (TMath::Abs(sn)<2*kSafe) {
2356     if (alpha>0) alpha += alpha< TMath::Pi()/2. ?  2*kSafe : -2*kSafe;
2357     else         alpha += alpha>-TMath::Pi()/2. ? -2*kSafe :  2*kSafe;
2358     cs=TMath::Cos(alpha);
2359     sn=TMath::Sin(alpha);
2360   }
2361   else if (TMath::Abs(cs)<2*kSafe) {
2362     if (alpha>0) alpha += alpha> TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
2363     else         alpha += alpha>-TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
2364     cs=TMath::Cos(alpha);
2365     sn=TMath::Sin(alpha);
2366   }
2367
2368
2369   return alpha;
2370 }
2371
2372 //---------------------------
2373 Double_t AliAnalysisTaskSELc2V0bachelor::Det(Double_t a00, Double_t a01,
2374                                              Double_t a10, Double_t a11) const {
2375   //--------------------------------------------------------------------
2376   // This function calculates locally a 2x2 determinant.
2377   // This is a copy of the AliCascadeVertexer::Det(...) method
2378   //--------------------------------------------------------------------
2379   return a00*a11 - a01*a10;
2380 }
2381
2382 //---------------------------
2383 Double_t AliAnalysisTaskSELc2V0bachelor::Det(Double_t a00,Double_t a01,Double_t a02,
2384                                              Double_t a10,Double_t a11,Double_t a12,
2385                                              Double_t a20,Double_t a21,Double_t a22) const {
2386   //--------------------------------------------------------------------
2387   // This function calculates locally a 3x3 determinant
2388   // This is a copy of the AliCascadeVertexer::Det(...) method
2389   //--------------------------------------------------------------------
2390   return  a00*Det(a11,a12,a21,a22)-a01*Det(a10,a12,a20,a22)+a02*Det(a10,a11,a20,a21);
2391 }
2392
2393 //----------------------------------------------------------------------------
2394 Int_t AliAnalysisTaskSELc2V0bachelor::MatchToMClabelC(AliAODRecoCascadeHF *candidate,
2395                                                       TClonesArray *mcArray)
2396 {
2397   //
2398   // Check if this candidate is matched to a MC signal  Lc -> p K0S + X
2399   // If no, return -1
2400   // If yes, return label (>=0) of the AliAODMCParticle
2401   // 
2402
2403   AliAODv0 *theV0 = dynamic_cast<AliAODv0*>(candidate->Getv0()); // the V0
2404   AliVTrack *trk = dynamic_cast<AliVTrack*>(candidate->GetBachelor()); // the bachelor
2405   if (!trk || !theV0) return -1;
2406
2407   if (trk->GetLabel()==-1) return -1;
2408   Int_t bachLabels = TMath::Abs(trk->GetLabel());
2409   AliAODMCParticle*bachelorMC = dynamic_cast<AliAODMCParticle*>(mcArray->At(bachLabels));
2410   if (!bachelorMC) return -1;
2411   if (TMath::Abs(bachelorMC->GetPdgCode())!=2212) return -1;
2412   Int_t indexMotherBach = bachelorMC->GetMother();
2413   if (indexMotherBach==-1) return -1;
2414
2415   Int_t pdgDg2prong[2] = {211,211};
2416   Int_t lab2Prong = theV0->MatchToMC(310,mcArray,2,pdgDg2prong); // the V0
2417   if(lab2Prong<0) return -1;
2418   AliAODMCParticle*partK0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(lab2Prong));
2419   if (!partK0S) return -1;
2420   Int_t indexMotherK0S = partK0S->GetMother();
2421   if (indexMotherK0S==-1) return -1;
2422   AliAODMCParticle*partK0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(indexMotherK0S));
2423   if (!partK0) return -1;
2424   Int_t indexMotherK0 = partK0->GetMother();
2425   if (indexMotherK0==-1) return -1;
2426
2427   if (indexMotherBach!=indexMotherK0) return -1; // p e K0S sono fratelli
2428
2429   AliAODMCParticle*partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(indexMotherK0));
2430   if (!partLc) return -1;
2431   Int_t ndg2 = partLc->GetDaughter(1)-partLc->GetDaughter(0)+1;
2432   if (ndg2==2) return -1;
2433
2434   Char_t stringaCheck[100];
2435   sprintf(stringaCheck,">>>>>>>> %d -> ",partLc->GetPdgCode());
2436   for(Int_t ii=0; ii<ndg2; ii++) {
2437     AliAODMCParticle* partDau=(AliAODMCParticle*)(mcArray->At(partLc->GetDaughter(0)+ii));
2438     sprintf(stringaCheck," %s %d",stringaCheck,partDau->GetPdgCode());
2439   }
2440   printf("%s \n",stringaCheck);
2441
2442   return indexMotherBach;
2443
2444 }
2445 //-----------------------------------------------------------------------------