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