PWGLF/SPECTRA/Kinks/AliAnalysisKinkESDMC.cxx

   1 /**************************************************************************
   2  * Authors: Martha Spyropoulou-Stassinaki and the  members
   3  * of the Greek group at Physics Department of Athens University
   4  * Paraskevi Ganoti, Anastasia Belogianni and Filimon Roukoutakis
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-----------------------------------------------------------------
  17 //                 AliAnalysisKinkESDMC class
  18 //       Example of an analysis task for kink topology study
  19 //      Kaons from kink topology are 'identified' in this code
  20 //-----------------------------------------------------------------
  21
  22 #include "TChain.h"
  23 #include "TTree.h"
  24 #include "TH1F.h"
  25 #include "TH2F.h"
  26 #include "TH3F.h"
  27 #include "TH1D.h"
  28 #include "TH2D.h"
  29 #include "TParticle.h"
  30 #include <TVector3.h>
  31 #include "TF1.h"
  32
  33 #include "AliAnalysisTask.h"
  34 #include "AliAnalysisManager.h"
  35 #include "AliTrackReference.h"
  36
  37 #include "AliVEvent.h"
  38 #include "AliESDEvent.h"
  39 #include "AliMCEvent.h"
  40 #include "AliAnalysisKinkESDMC.h"
  41 #include "AliStack.h"
  42 #include "AliESDpid.h"
  43 #include "AliESDkink.h"
  44 #include "AliESDtrack.h"
  45 #include "AliESDtrackCuts.h"
  46 #include "AliPhysicsSelectionTask.h"
  47 #include "AliInputEventHandler.h"
  48 #include "AliESDInputHandler.h"
  49 #include "AliAnalysisManager.h"
  50 #include "AliVEvent.h"
  51  #include "AliPIDResponse.h"
  52 ///#include "AliTPCpidESD.h"
  53
  54 ClassImp(AliAnalysisKinkESDMC)
  55
  56
  57 //________________________________________________________________________
  58 AliAnalysisKinkESDMC::AliAnalysisKinkESDMC(const char *name)
  59   : AliAnalysisTaskSE(name), fHistPtESD(0),fHistPt(0),fHistQtAll(0),fHistQt1(0),fHistQt2(0)
  60   , fHistPtKaon(0),fHistPtKPDG(0),fHistEta(0),fHistEtaK(0),fptKMC(0),fMultiplMC(0),fESDMult(0),frad(0),
  61   fradMC(0), fKinkKaon(0), fKinkKaonBg(0), fM1kaon(0),  fgenPtEtR(0),fPtKink(0),
  62    fcodeH(0), fdcodeH(0), fAngMomK(0),fAngMomPi(0), fAngMomKC(0),  fMultESDK(0), fMultMCK(0),
  63  fSignPtNcl(0), fSignPtEta(0), fSignPtEtaMC(0), fSignPtMC(0),  fEtaNcl(0), fSignPt(0), fChi2NclTPC(0), fRatChi2Ncl(0),
  64   fRadiusNcl(0),  fTPCSgnlP(0), fTPCSgnlPa(0),  fSignPtGen(0),
  65   fRpr(0),fZpr(0), fdcatoVxXY(0),    fMCEtaKaon(0),
  66  fZvXv(0),fZvYv(0),fXvYv(0),fPtPrKink(0),fgenPtEtRP(0),fgenPtEtRN(0),fkinkKaonP(0),fkinkKaonN(0),
  67    frapidESDK(0), frapidKMC(0),  fPtKPlMC(0), fPtKMnMC(0),
  68      fHistPtKaoP(0), fHistPtKaoN(0), fHiPtKPDGP(0), fHiPtKPDGN(0),fKinKBGP(0),fKinKBGN(0),
  69     fQtKMu(0),fQtKPi(0),fQtKEl(0),fFakepipi(0), fFakeKPi(0),
  70      fDCAkink(0), fDCAkinkBG(0), fPosiKink(0),  fPosiKinkK(0), fPosiKinKXZ(0), fPosiKinKYZ(0),  fPosiKinKBgZY(0),
  71     fcode2(0), fcode4(0), fZkinkZDau(0),
  72          fQtKMuMC(0),  fQtKElMC(0), fQtKPiMC(0),   fQtK3PiP(0),fQtK3PiM(0),  fmaxAngMomKmu(0),
  73            fPosiKinKBgZX(0), fPosiKinKBgXY(0),  fMinvPi(0),fMinvKa(0),fMinvPr(0),
  74                         fTPCSgnlPtpc(0),
  75        fTPCMomNSgnl(0),  fMothKinkMomSgnl(0), fNSigmTPC(0),  fTPCSgnlKinkDau(0),fcodeDau1(0),fcodeDau2(0), fMothKinkMomSgnlD(0),
  76      fInvMassMuNuAll(0),   fInvMassMuNuPt(0), fInvMassMuNuPtAll(0), fRatioCrossedRows(0), fRatioCrossedRowsKink(0), fRadiusPt(0), fRadiusPtcln(0),
  77       fRadiusPtKaon(0), fRadiusPtPion(0), fRadiusPtFake(0), fPtCut1(0), fPtCut2(0), fPtCut3(0),  fAngMomKKinks(0),
  78   flengthMCK(0), flifetiMCK(0), flifetim2(0), fLHelESDK(0),flifeInt(0), flifeYuri(0), flenYuri(0), flenTrRef(0),flifeSmall(0), flifetime(0),flifTiESDK(0),
  79     flifeKink(), flenHelx(0), fradPtRapMC(0), fradPtRapDC(0), fradPtRapESD(0), fRadNclcln(0),
  80     f1(0), f2(0),
  81   fListOfHistos(0),fLowMulcut(-1),fUpMulcut(-1), fKinkRadUp(200),fKinkRadLow(130), fLowCluster(20), fLowQt(.12),  fCutsMul(0),fMaxDCAtoVtxCut(0), fPIDResponse(0)
  82
  83 {
  84   // Constructor
  85
  86   // Define input and output slots here
  87
  88         DefineOutput(1, TList::Class());
  89 }
  90
  91 //________________________________________________________________________
  92 void AliAnalysisKinkESDMC::UserCreateOutputObjects()
  93 {
  94
  95   // Input slot #0 works with a TChain
  96   // DefineInput(0, TChain::Class());liESDtrackCuts("Mul","Mul");
  97         fCutsMul=new AliESDtrackCuts("Mul","Mul");
  98         fCutsMul->SetMinNClustersTPC(70);
  99         fCutsMul->SetMaxChi2PerClusterTPC(4);
 100         fCutsMul->SetAcceptKinkDaughters(kFALSE);
 101         fCutsMul->SetRequireTPCRefit(kTRUE);
 102         // ITS
 103         fCutsMul->SetRequireITSRefit(kTRUE);
 104         fCutsMul->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 105                                                 AliESDtrackCuts::kAny);
 106         fCutsMul->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
 107
 108         fCutsMul->SetMaxDCAToVertexZ(2);
 109         fCutsMul->SetDCAToVertex2D(kFALSE);
 110         fCutsMul->SetRequireSigmaToVertex(kFALSE);
 111
 112         fCutsMul->SetEtaRange(-0.8,+0.8);
 113         fCutsMul->SetPtRange(0.15, 1e10);
 114 //
 115                    fMaxDCAtoVtxCut=new AliESDtrackCuts("fMaxDCA", "fMaxDCA");
 116        fMaxDCAtoVtxCut->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
 117     fMaxDCAtoVtxCut->SetMaxChi2TPCConstrainedGlobal(36);
 118
 119         // Create histograms
 120         // Called once
 121
 122         f1=new TF1("f1","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",1.1,10.0);
 123         f1->SetParameter(0,0.493677);
 124         f1->SetParameter(1,0.9127037);
 125         f1->SetParameter(2,TMath::Pi());
 126
 127
 128         f2=new TF1("f2","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",0.1,10.0);
 129         f2->SetParameter(0,0.13957018);
 130         f2->SetParameter(1,0.2731374);
 131         f2->SetParameter(2,TMath::Pi());
 132 //
 133       Double_t gPt7K0[45] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8,0.9,1.0,
 134                         1.1, 1.2, 1.3, 1.4, 1.5, 1.6,1.7,1.8,1.9,  2.0,
 135                          2.2, 2.4, 2.6, 2.8,  3.0,   3.3, 3.6, 3.9,
 136                          4.2, 4.6,5.0, 5.4, 5.9,  6.5,   7.0,7.5, 8.0,8.5,  9.2, 10., 11., 12., 13.5,15.0 };  // David K0
 137 //
 138                   Double_t gPt7Comb[48] = {
 139 //! ! ! ! !  KINK FROM HERE --------------->
 140 0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0, 1.1, 1.2,
 141 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9,
 142 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8,5.0, 5.5, 6
 143  };  // 25/11/2013 from Francesco
 144 /*
 145    Double_t gPt7TOF[47] = { 0.2,0.25, 0.3,0.35,  0.4,0.45,  0.5,0.55,  0.6,0.65,  0.7,0.75,  0.8, 0.85, 0.9, 0.95, 1.0,
 146                         1.1, 1.2, 1.3, 1.4, 1.5, 1.6,1.7,1.8,1.9,  2.0,
 147                          2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7,2.8, 2.9, 3.0,
 148                          3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8,5.0 };  //  Barbara TOF Kch
 149 */
 150         fHistPtESD = new TH1F("fHistPtESD", "P_{T} distribution",200, 0.0,10.0);
 151         fHistPt = new TH1F("fHistPt", "P_{T} distribution",200, 0.0,10.0);
 152         fHistQtAll = new TH1F("fHistQtAll", "Q_{T} distr All Kinks ",100, 0.0,.300);
 153         fHistQt1= new TH1F("fHistQt1", "Q_{T} distribution",100, 0.0,.300);
 154         fHistQt2= new TH1F("fHistQt2", "Q_{T} distribution",100, 0.0,.300);
 155 //      fHistPtKaon = new TH1F("fHistPtKaon", "P_{T}Kaon distribution",200, 0.0,10.0);
 156         fHistPtKaon = new TH1F("fHistPtKaon", "P_{T}Kaon distribution",47,gPt7Comb );
 157         fHistPtKPDG = new TH1F("fHistPtKPDG", "P_{T}Kaon distribution",47, gPt7Comb   );
 158         fHistEta= new TH1F("fHistEta", "Eta distribution", 26,-1.3, 1.3);
 159         fHistEtaK= new TH1F("fHistEtaK", "EtaK distribution", 26,-1.3, 1.3);
 160         fptKMC= new TH1F("fptKMC", "P_{T}Kaon generated",47,  gPt7Comb  );
 161         fMultiplMC= new TH1F("fMultiplMC", " charge particle multipl",100, 0.0,2500.);
 162         fESDMult= new TH1F("fESDMult", "charge multipliESD",100, 0.0,100.);
 163         frad= new TH1F("frad", "radius  K ESD recon",100,0.,1000.);
 164         fradMC= new TH1F("fradMC", "radius  K generated",100,0.,1000.);
 165         fKinkKaon= new TH1F("fKinkKaon", "P_{T}Kaon kinks identi", 47, gPt7Comb  );
 166         fKinkKaonBg= new TH1F("fKinkKaonBg", "P_{T}Kaon kinks backgr",47 , gPt7Comb );
 167         //fM1kaon= new TH1F("fM1kaon","Invar m(kaon) from kink->mu+netrino decay",180,0.1, 1.0);
 168         fM1kaon= new TH1F("fM1kaon","Invar m(kaon) from kink->mu+netrino decay",600,0.1, 0.7);
 169         fgenPtEtR= new TH1F("fgenPtEtR", "P_{T}Kaon distribution", 47, gPt7Comb  );
 170         fPtKink= new TH1F("fPtKink", "P_{T}Kaon Kink  bution",47, gPt7Comb  );
 171         fcodeH   = new TH2F("fcodeH", "code vrs dcode dist. kinks,K",100,0.,2500.,100,0.,2500.);
 172         fdcodeH = new TH2F("fdcodeH", "code vrs dcode dist. kinks,K",100,0.,2500.,100,0.,2500.);
 173         fAngMomK= new TH2F("fAngMomK","Decay angle vrs Mother Mom,K",100,0.0,10.0,80,0.,80.);
 174         fAngMomPi= new TH2F("fAngMomPi","Decay angle vrs Mother Mom,Pi",100,0.0,10.0,80,0.,80.);
 175         fAngMomKC= new TH2F("fAngMomKC","Decay angle vrs Mother Mom,K",100,0.0,10.0,80,0.,80.);
 176         fMultESDK=new TH1F("fMultESDK", "charge multipliESD kaons",100, 0.,2500.);
 177         fMultMCK=new TH1F("fMultMCK", "charge multipli MC kaons",100, 0.,2500.);
 178         fSignPtNcl= new TH2F("fSignPtNcl","SignPt vrs Ncl,K",80,-4.,4.0,70,20.,160.);
 179         fSignPtEta= new TH2F("fSignPtEta","SignPt vrs Eta,K",80,-4.0,4.0,30,-1.5,1.5);
 180         fSignPtEtaMC= new TH2F("fSignPtEtaMC","SignPt vrs Eta,K",80,-4.0,4.0,30,-1.5,1.5);
 181         fSignPtMC= new TH1F("fSignPtMC","SignPt ,K",100,-5.,5.0);
 182         fEtaNcl= new TH2F("fEtaNcl","Eta vrs Ncl,K",26,-1.3,1.3,70,20.,160.);
 183         fSignPt= new TH1F("fSignPt","SignPt ,K",100,-5.,5.0);
 184         fChi2NclTPC= new TH2F("fChi2NclTPC","Chi2vrs nclust,K",100,0.,500., 70,20, 160);
 185         fRatChi2Ncl= new TH1F("fRatChi2Ncl","Ratio chi2/nclusters in TPC,K",50,0.0,5.0);
 186         fRadiusNcl = new TH2F("fRadiusNcl","KinkRadius Ncl,K",75,100.,250., 80,0, 160);
 187         fTPCSgnlP = new TH2F("fTPCSgnlP","Kink TCP de/dx,K",300,0.,15.,150,0.,300);
 188         fTPCSgnlPa= new TH2F("fTPCSgnlPa","Kink TCP de/dx,a",300,0.,15.,150,0.,300);
 189         fSignPtGen= new TH1F("fSignPtGen","SignPtGen ,K",100,-5.0,5.0);
 190         fRpr = new TH1D("fRpr", "rad distribution  PID pr",100,-1.0,1.0);
 191         fZpr = new TH1D("fZpr", "z distribution PID pr  ",60,-15.,15.);
 192         fdcatoVxXY = new TH1D("fdcatoVxXY", "dca  distribution PID  ",20,-1.,1.);
 193         fMCEtaKaon = new TH1F("fMCEtaKaon"," Hist of Eta K -Kink Selecied",26,-1.3,1.3);
 194         fZvXv= new TH2F("fZvXv","Xv-Zv main vtx",60,-0.5,0.5, 60, -15., 15.);
 195         fZvYv= new TH2F("fZvYv","Yv-Zv main vtx",60,-0.5,0.5, 60, -15., 15.);
 196         fXvYv= new TH2F("fXvYv","Xv-Yv main vtx", 60,-1.5,1.5, 60, -1.5, 1.5);
 197         fPtPrKink=new TH1F("fPtPrKink","pt of ESD  kaonKink tracks",300, 0.0,15.0);
 198         fgenPtEtRP= new TH1F("fgenPtEtRP", "P_{T}Kaon distribution positive", 47, gPt7Comb  );
 199         fgenPtEtRN= new TH1F("fgenPtEtRN", "P_{T}Kaon distribution negative", 47, gPt7Comb  );
 200         fkinkKaonP= new TH1F("fKinkKaonP", "P_{T}Kaon distribution positive", 47, gPt7Comb  );
 201         fkinkKaonN= new TH1F("fKinkKaonN", "P_{T}Kaon distribution negative", 47, gPt7Comb  );
 202         frapidESDK= new TH1F("frapidESDK", "rapidity distribution", 26,-1.3, 1.3);
 203         frapidKMC = new TH1F("frapidKMC ", "rapidity distri  MC  ",26,-1.3, 1.3);
 204         fPtKPlMC= new TH1F("fPtKPlMC", "P_{T}Kaon Pos  generated", 47, gPt7Comb  );
 205         fPtKMnMC= new TH1F("fPtKMnMC", "P_{T}Kaon Minusgenerated",47 , gPt7Comb  );
 206         fHistPtKaoP= new TH1F("fHistPtKaoP", "P_{T}Kaon Pos ESD", 47, gPt7Comb   );
 207         fHistPtKaoN= new TH1F("fHistPtKaoN", "P_{T}Kaon Neg ESD", 47, gPt7Comb   );
 208         fHiPtKPDGP= new TH1F("fHiPtKPDGP", "P_{T}Kaon Pos ESD", 47,  gPt7Comb  );
 209         fHiPtKPDGN= new TH1F("fHiPtKPDGN", "P_{T}Kaon neg ESD", 47, gPt7Comb   );
 210         fKinKBGP  = new TH1F("fKinKBGP  ", "P_{T}Kaon Pos ESD", 47, gPt7Comb   );
 211         fKinKBGN= new TH1F("fKinKBGN", "P_{T}Kaon neg ESD", 47, gPt7Comb   );
 212         fQtKMu= new TH1F("fQtKMu", "Q_{T} distribution  K to mu ",100, 0.0,.300);
 213         fQtKPi= new TH1F("fQtKPi", "Q_{T} distribution K to pi",100, 0.0,.300);
 214         fQtKEl= new TH1F("fQtKEl", "Q_{T} distribution   K to elec",100, 0.0,.300);
 215         fFakepipi = new TH1F("fFakepipi", "P_{T}fake pipi   ",47 , gPt7Comb   );
 216         fFakeKPi = new TH1F("fFakeKPi", "P_{T}fake Kpi   ", 47, gPt7Comb   );
 217         fDCAkink = new TH1F("fDCAkink", "DCA kink vetrex ",50, 0.0,1.0);
 218         fDCAkinkBG = new TH1F("fDCAkinkBG", "DCA kink vetrex ",50, 0.0,1.0);
 219         fPosiKink= new TH2F("fPosiKink", "Y vrx kink Vrex ",100, -300.0,300.0,100, -300, 300.);
 220         fPosiKinkK= new TH2F("fPosiKinkK", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
 221         fPosiKinKXZ= new TH2F("fPosiKinKXZ", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
 222         fPosiKinKYZ= new TH2F("fPosiKinKYZ", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
 223         fPosiKinKBgZY= new TH2F("fPosiKinKBgZY", "Y vrx Z kink Vrexbg ",100, -300.0,300.0,100, -300, 300.);
 224         fcode2   = new TH2F("fcode2", "code vrs dcode dist. kinks,K",100,0.,2500.,100,0.,2500.);
 225         fcode4   = new TH2F("fcode4", "code vrs dcode dist. kinks,K",100,0.,2500.,100,0.,2500.);
 226         fZkinkZDau = new TH2F("fZkinkZDau", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
 227   fQtKMuMC= new TH1F("fQtKMuMC", "Q_{T} distribution  K to mu MC",100, 0.0,.300);
 228   fQtKPiMC= new TH1F("fQtKPiMC", "Q_{T} distribution K to pi MC",100, 0.0,.300);
 229   fQtKElMC= new TH1F("fQtKElMC", "Q_{T} distribution   K to elec MC",100, 0.0,.300);
 230   fQtK3PiP= new TH1F("fQtK3PiP", "Q_{T} distribution K to 3pi ",100, 0.0,.300);
 231   fQtK3PiM= new TH1F("fQtK3PiM", "Q_{T} distribution K to 3pi ",100, 0.0,.300);
 232   fmaxAngMomKmu= new TH2F("fmaxAngMomKmu","Decay angle vrs Mother Mom,Kmu",100,0.0,10.0,80,0.,80.);
 233         fPosiKinKBgZX= new TH2F("fPosiKinKBgZX", "X vrx Z kink Vrexbg ",100, -20.0,20.0,100, 0., 300.);
 234         fPosiKinKBgXY= new TH2F("fPosiKinKBgXY", "Y vrx X kink Vrexbg ",100, -300.0,300.0,100, -300, 300.);
 235         fMinvPi= new TH1F("fMinvPi","Invar m(kaon) from kink-> decay",100,0.0, 1.2);
 236         fMinvKa= new TH1F("fMinvKa","Invar m(kaon) from kink-> decay",100,0.0, 2.0);
 237         fMinvPr= new TH1F("fMinvPr","Invar m(kaon) from kink-> decay",100,0.0, 1.2);
 238                 fTPCSgnlPtpc = new TH2F("fTPCSgnlPtpc","TPC signal de/dx Mom TPC,K  ",100,0.0,8.0,100, 0., 250.    );
 239     fTPCMomNSgnl = new TH2F("fTPCMomNsgnl","TPC signal de/dx Mom TPC,K  ",100,0.0,8.0,20 , -10., 10.);
 240     fMothKinkMomSgnl  = new TH2F("fMothKinkMomSgnl","TPC signal de/dx Mom TPC,Kink  ",100,0.0,250.,100, 0., 250.    );
 241     fNSigmTPC    = new TH1F("fNSigmTPC","TPC Nsigma  de/dx  TPC,K  ", 30 , -7.5, 7.5);
 242     fTPCSgnlKinkDau = new TH2F("fTPCSgnlKinkDau","TPC signal de/dx Mom,K",100,0.0,8.0,100,0.,250.);
 243         fcodeDau1   = new TH2F("fcodeDau1", "code vrs dcode dist. kinks,K",100,0.,500.,100,0.,500.);
 244         fcodeDau2   = new TH2F("fcodeDau2", "code vrs dcode dist. kinks,K",100,0.,500.,100,0.,500.);
 245     fMothKinkMomSgnlD = new TH2F("fMothKinkMomSgnlD","TPC signal de/dx Mom TPC,Kink  ",100,0.0,250.,100, 0., 250.    );
 246         //fInvMassMuNuAll = new TH1F("fInvMassMuNuAll","Invar from kink->mu+netrino decay",180,0.1, 1.0);
 247         fInvMassMuNuAll = new TH1F("fInvMassMuNuAll","Invar from kink->mu+netrino decay",600,0.1, 0.7); //  22/8/2013
 248         //fInvMassMuNuPt  = new TH2F("fInvMassMuNuPt","Invar from kink->mu+netrino decay vs Pt",180,0.1, 1.0, 100, 0. , 10.);
 249         fInvMassMuNuPt  = new TH2F("fInvMassMuNuPt","Invar from kink->mu+netrino decay vs Pt",600,0.1, 0.7, 100, 0. , 10.); //  22/8/2013
 250         fInvMassMuNuPtAll  = new TH2F("fInvMassMuNuPtAll","Invar from kink->mu+netrino decay vs Pt",600,0.1, 0.7, 100, 0. , 10.); //  22/8/2013
 251       fRatioCrossedRows = new TH1F("fRatioCrossedRows","Ratio crossed rows  in TPC",20,0.0,1.0);
 252   fRatioCrossedRowsKink = new TH1F("fRatioCrossedRowsKink","Ratio crossed rows  in TPC for kinks",20,0.0,1.0);
 253   fRadiusPt =new TH2F("fRadiusPt","radius vs pt  ",80, 90.,250.,100, 0.,10.              );
 254   fRadiusPtcln =new TH2F("fRadiusPtcln","radius vs pt clean ",80, 90.,250.,100, 0.,10.              );
 255   fRadiusPtKaon  =new TH2F("fRadiusPtKaon","radius vs pt Kaon PDG ",80, 90.,250.,100, 0.,10.              );
 256   fRadiusPtPion  =new TH2F("fRadiusPtPion","radius vs pt Pion PDG ",80, 90.,250.,100, 0.,10.              );
 257   fRadiusPtFake  =new TH2F("fRadiusPtFake","radius vs pt Pion Fake ",80, 90.,250.,100, 0.,10.              );
 258   fPtCut1 = new TH1F("fPtCut1", "P_{T}Kaon distribution",300, 0.0,15.0);
 259   fPtCut2 = new TH1F("fPtCut2", "P_{T}Kaon distribution",300, 0.0,15.0);
 260   fPtCut3 = new TH1F("fPtCut3", "P_{T}Kaon distribution",300, 0.0,15.0);
 261   fAngMomKKinks = new TH2F("fAngMomKKinks","Decay angle vrs Mother Mom,Kinks",300,0.0,15.0,100,0.,100.);
 262    flengthMCK=new TH1F("flengthMCK", "length of K  MCref decay ",100,0.,1000.0);
 263   flifetiMCK=new TH1F("flifetiMCK", "lifetime ref K   Decay   ",100,0.,1000.0);
 264   flifetim2 =new TH1F("flifetim2", "lifetime ref K   Decay   ",100,0.,1000.0);
 265   fLHelESDK =new TH1F("fLHelESDK", "lifetime ref K   Decay   ",100,0.,1000.0);
 266   flifeInt =new TH1F("flifeInt", "lifetime ref K   Decay   ",100,0.,1000.0);
 267   flifeYuri=new TH1F("flifeYuri","lifetime ref K   Decay   ",100,0.,1000.0);
 268   flenYuri=new TH1F("flenYuri","lifetime ref K   Decay   ",100,0.,1000.0);
 269   flenTrRef =new TH1F("flenTrRef","lifetime ref K   Decay   ",100,0.,1000.0);
 270   flifeSmall=new TH1F("flifeSmall","lifetime ref K   Decay   ",100,0.,1000.0);
 271   flifetime =new TH1F("flifetime","lifetime ref K   Decay   ",100,0.,1000.0);
 272   flifTiESDK=new TH1F("flifTiESDK","lifetime ref K   Decay   ",100,0.,1000.0);
 273   flifeKink =new TH1F("flifeKink", "lifetime ref K   Decay   ",100,0.,1000.0);
 274   flenHelx  =new TH1F("flenHelx", "lifetime ref K   Decay   ",100,0.,1000.0);
 275      fradPtRapMC=new TH3F("fradPtRapMC","rad pt rap dat",28,100.,240., 20, 0., 5., 20, -1., 1. );
 276      fradPtRapDC=new TH3F("fradPtRapDC","rad pt rap dat",28,100.,240., 20, 0., 5., 20, -1., 1. );
 277      fradPtRapESD=new TH3F("fradPtRapESD","rad pt rap dat",28,100.,240., 20, 0., 5., 20, -1., 1. );
 278         fRadNclcln = new TH2F("fRadNclcln","KinkRadius Ncl,K",75,100.,250., 80,0, 160);
 279
 280
 281
 282
 283    fListOfHistos=new TList();
 284    fListOfHistos->SetOwner();
 285    fListOfHistos->Add(fHistPtESD);
 286    fListOfHistos->Add(fHistPt);
 287    fListOfHistos->Add(fHistQtAll);
 288    fListOfHistos->Add(fHistQt1);
 289    fListOfHistos->Add(fHistQt2);
 290    fListOfHistos->Add(fHistPtKaon);
 291    fListOfHistos->Add(fHistPtKPDG);
 292    fListOfHistos->Add(fHistEta);
 293    fListOfHistos->Add(fHistEtaK);
 294    fListOfHistos->Add(fptKMC);
 295    fListOfHistos->Add(fMultiplMC);
 296    fListOfHistos->Add(fESDMult);
 297    fListOfHistos->Add(frad);
 298    fListOfHistos->Add(fradMC);
 299    fListOfHistos->Add(fKinkKaon);
 300    fListOfHistos->Add(fKinkKaonBg);
 301    fListOfHistos->Add(fM1kaon);
 302    fListOfHistos->Add(fgenPtEtR);
 303    fListOfHistos->Add(fPtKink);
 304    fListOfHistos->Add(fcodeH);
 305    fListOfHistos->Add(fdcodeH);
 306    fListOfHistos->Add(fAngMomK);
 307    fListOfHistos->Add(fAngMomPi);
 308    fListOfHistos->Add(fAngMomKC);
 309    fListOfHistos->Add(fMultESDK);
 310    fListOfHistos->Add(fMultMCK);
 311    fListOfHistos->Add(fSignPtNcl);
 312    fListOfHistos->Add(fSignPtEta);
 313    fListOfHistos->Add(fSignPtEtaMC);
 314    fListOfHistos->Add(fSignPtMC);
 315    fListOfHistos->Add(fEtaNcl);
 316    fListOfHistos->Add(fSignPt);
 317    fListOfHistos->Add(fChi2NclTPC);
 318    fListOfHistos->Add(fRatChi2Ncl);
 319    fListOfHistos->Add(fRadiusNcl);
 320    fListOfHistos->Add(fTPCSgnlP);
 321    fListOfHistos->Add(fTPCSgnlPa);
 322    fListOfHistos->Add(fSignPtGen);
 323    fListOfHistos->Add(fRpr);
 324    fListOfHistos->Add(fZpr);
 325    fListOfHistos->Add(fdcatoVxXY);
 326    fListOfHistos->Add(fMCEtaKaon);
 327    fListOfHistos->Add(fZvXv);
 328    fListOfHistos->Add(fZvYv);
 329    fListOfHistos->Add(fXvYv);
 330    fListOfHistos->Add(fPtPrKink);
 331    fListOfHistos->Add(fgenPtEtRP);
 332    fListOfHistos->Add(fgenPtEtRN);
 333    fListOfHistos->Add(fkinkKaonP);
 334    fListOfHistos->Add(fkinkKaonN);
 335    fListOfHistos->Add(frapidESDK);
 336    fListOfHistos->Add(frapidKMC);
 337    fListOfHistos->Add(fPtKPlMC);
 338    fListOfHistos->Add(fPtKMnMC);
 339    fListOfHistos->Add(fHistPtKaoP);
 340    fListOfHistos->Add(fHistPtKaoN);
 341    fListOfHistos->Add(fHiPtKPDGP);
 342    fListOfHistos->Add(fHiPtKPDGN);
 343    fListOfHistos->Add(fKinKBGP);
 344    fListOfHistos->Add(fKinKBGN);
 345    fListOfHistos->Add(fQtKMu);
 346    fListOfHistos->Add(fQtKPi);
 347    fListOfHistos->Add(fQtKEl);
 348    fListOfHistos->Add(fFakepipi);
 349    fListOfHistos->Add(fFakeKPi);
 350    fListOfHistos->Add(fDCAkink);
 351    fListOfHistos->Add(fDCAkinkBG);
 352    fListOfHistos->Add(fPosiKink);
 353    fListOfHistos->Add(fPosiKinkK);
 354    fListOfHistos->Add(fPosiKinKXZ);
 355    fListOfHistos->Add(fPosiKinKYZ);
 356    fListOfHistos->Add(fPosiKinKBgZY);
 357    fListOfHistos->Add(fcode2);
 358    fListOfHistos->Add(fcode4);
 359    fListOfHistos->Add(fZkinkZDau);
 360    fListOfHistos->Add(fQtKMuMC);
 361    fListOfHistos->Add(fQtKPiMC);
 362    fListOfHistos->Add(fQtKElMC);
 363    fListOfHistos->Add(fQtK3PiP);
 364    fListOfHistos->Add(fQtK3PiM);
 365    fListOfHistos->Add(fmaxAngMomKmu);
 366    fListOfHistos->Add(fPosiKinKBgZX);
 367    fListOfHistos->Add(fPosiKinKBgXY);
 368    fListOfHistos->Add(fMinvPi);
 369    fListOfHistos->Add(fMinvKa);
 370    fListOfHistos->Add(fMinvPr);
 371    fListOfHistos->Add(fTPCSgnlPtpc);
 372    fListOfHistos->Add(fTPCMomNSgnl);
 373    fListOfHistos->Add(fMothKinkMomSgnl);
 374    fListOfHistos->Add(fNSigmTPC);
 375    fListOfHistos->Add(fTPCSgnlKinkDau);
 376    fListOfHistos->Add(fcodeDau1);
 377    fListOfHistos->Add(fcodeDau2);
 378    fListOfHistos->Add(fMothKinkMomSgnlD);
 379    fListOfHistos->Add(fInvMassMuNuAll);
 380    fListOfHistos->Add(fInvMassMuNuPt);
 381    fListOfHistos->Add(fInvMassMuNuPtAll);
 382    fListOfHistos->Add(fRatioCrossedRows);
 383    fListOfHistos->Add(fRatioCrossedRowsKink);
 384    fListOfHistos->Add(fRadiusPt);
 385    fListOfHistos->Add(fRadiusPtcln);
 386    fListOfHistos->Add(fRadiusPtKaon);
 387    fListOfHistos->Add(fRadiusPtPion);
 388    fListOfHistos->Add(fRadiusPtFake);
 389    fListOfHistos->Add(fPtCut1);
 390    fListOfHistos->Add(fPtCut2);
 391    fListOfHistos->Add(fPtCut3);
 392    fListOfHistos->Add(fAngMomKKinks);
 393    fListOfHistos->Add(flengthMCK);
 394    fListOfHistos->Add(flifetiMCK);
 395    fListOfHistos->Add(flifetim2);
 396    fListOfHistos->Add(fLHelESDK);
 397    fListOfHistos->Add(flifeInt);
 398    fListOfHistos->Add(flifeYuri);
 399    fListOfHistos->Add(flenYuri);
 400    fListOfHistos->Add(flenTrRef);
 401    fListOfHistos->Add(flifeSmall);
 402    fListOfHistos->Add(flifetime);
 403    fListOfHistos->Add(flifTiESDK);
 404    fListOfHistos->Add(flifeKink);
 405    fListOfHistos->Add(flenHelx);
 406    fListOfHistos->Add(fradPtRapMC);
 407    fListOfHistos->Add(fradPtRapDC);
 408    fListOfHistos->Add(fradPtRapESD);
 409    fListOfHistos->Add(fRadNclcln);
 410
 411   PostData(1, fListOfHistos);
 412 }
 413
 414 //________________________________________________________________________
 415 void AliAnalysisKinkESDMC::UserExec(Option_t *)
 416 {
 417   // Main loop
 418   // Called for each event
 419
 420   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 421   // This handler can return the current MC event
 422
 423    AliVEvent *event = InputEvent();
 424   if (!event) {
 425      Printf("ERROR: Could not retrieve event");
 426      return;
 427   }
 428
 429   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 430   if (!esd) {
 431      Printf("ERROR: Could not retrieve esd");
 432      return;
 433 }
 434
 435   AliMCEvent* mcEvent = MCEvent();
 436   if (!mcEvent) {
 437      Printf("ERROR: Could not retrieve MC event");
 438      return;
 439   }
 440      fMultMCK->Fill(mcEvent->GetNumberOfTracks() );
 441 //
 442
 443
 444 //   multiplicity selection
 445    Float_t refmultiplicity=fCutsMul->CountAcceptedTracks(esd);
 446         if(fLowMulcut>-1)
 447         {
 448                 if(refmultiplicity<fLowMulcut)
 449                         return;
 450         }
 451         if(fUpMulcut>-1)
 452         {
 453                 if(refmultiplicity>fUpMulcut)
 454                         return;
 455         }
 456
 457
 458        fMultESDK->Fill(refmultiplicity);
 459
 460
 461
 462
 463   Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
 464
 465   AliStack* stack=mcEvent->Stack();
 466
 467 //primary tracks  in MC
 468          Int_t  nPrim = stack->GetNprimary();
 469 //
 470
 471   // loop over mc particles
 472
 473   // variable to count tracks
 474   Int_t nMCTracks = 0;
 475   Int_t nMCKinkKs = 0;
 476
 477 //   15/2/12 OK   for (Int_t iMc = 0; iMc < nPrim; iMc++)
 478  //for (Int_t iMc = 0; iMc < stack->GetNtrack(); iMc++)
 479  for (Int_t iMc = 0; iMc < mcEvent->GetNumberOfTracks(); iMc++)
 480   {
 481
 482     TParticle* particle = stack->Particle(iMc);
 483
 484     if (!particle)
 485     {
 486     //  AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
 487       continue;
 488     }
 489 // keep only primaries
 490   if (!stack->IsPhysicalPrimary(iMc)) continue;
 491
 492  //
 493
 494             Double_t ptK = particle->Pt();
 495
 496      if( ptK <0.200) continue;       //    12/2/2012
 497 //
 498
 499                 Float_t charg=0;
 500       Float_t code = particle->GetPdgCode();
 501             Int_t  mcProcess=-1011;
 502 //---------------------------------------kaon selection
 503       if ((code==321)||(code==-321)){
 504
 505
 506           Double_t   etracKMC= TMath::Sqrt(particle->P() *particle->P()  + 0.493677 *0.493677  );
 507          Double_t rapidiKMC = 0.5 * (TMath::Log(  (etracKMC +particle->Pz())/( etracKMC-particle->Pz() )) )  ;
 508
 509      if ( TMath::Abs( rapidiKMC) > 0.7) continue;   //
 510             frapidKMC ->Fill(rapidiKMC) ;  //18/feb rapiddistr of PDG kink ESD  kaons
 511
 512
 513 //  maximum angle    vrs momentum
 514   //   Double_t maxDecAnKmu=f1->Eval(particle->P(),      0.,0.,0.);
 515  //       fmaxAngMomKmu->Fill(particle->P() , maxDecAnKmu);
 516
 517                if(code > 0 ) charg =1;
 518                if(code < 0 ) charg =-1;
 519                          Float_t chargPt= charg*ptK;
 520
 521               fptKMC->Fill(ptK);
 522          fSignPtGen->Fill(chargPt);// kaon gensign pt
 523               if (charg==1 )  fPtKPlMC->Fill( ptK );
 524                if ( charg==-1) fPtKMnMC->Fill( ptK  );
 525 // primary   vertex
 526             //      Double_t mVx=particle->Vx();
 527             //     Double_t mVy=particle->Vy();
 528                  Double_t mVz=particle->Vz();
 529 // 25/11/2012  ???????back 10/1/2013
 530                                                 TClonesArray* trArray=0;
 531                                                 TParticle* tempParticle=0;
 532
 533                                             TVector3 DecayMomentumK(0,0,0);
 534                                               Float_t lengthKMC=0;
 535                       if (mcEvent->GetParticleAndTR(iMc, tempParticle, trArray) != -1) {
 536                                                 AliTrackReference* MCtrackReference = static_cast<AliTrackReference*>(trArray->Last());
 537                                                    lengthKMC = MCtrackReference->GetLength();
 538
 539 //                      DecayMomentumK.SetXYZ(MCtrackReference->Px(), MCtrackReference->Py(), MCtrackReference->Pz());
 540                                               }
 541                    flenTrRef ->Fill(lengthKMC);
 542          flifetime ->Fill(  (lengthKMC*0.493667  /particle->P()));  // 19/7
 543                  if ((lengthKMC>100.)&& (lengthKMC<300.) )  flifeSmall->Fill( (lengthKMC*0.493667/particle->P() ) );
 544
 545        Int_t nMCKpi =0;
 546        Int_t mcProc4 =0;
 547        Int_t mcProc13=0;
 548         Double_t radiusD=0;
 549    //    Double_t lengthK =0.;
 550        Double_t LengthK =0.;
 551        Double_t lenYuri =0.;
 552         Double_t MCQt =0.;
 553         Double_t MCQt3[2];
 554             Int_t firstD=particle->GetFirstDaughter();
 555             Int_t lastD=particle->GetLastDaughter();
 556
 557              if( (lastD<=0) || (firstD<=0)  ) continue;
 558
 559                      if ( lastD > mcEvent->GetNumberOfTracks() ) continue;
 560                      if (firstD > mcEvent->GetNumberOfTracks() ) continue;
 561 // 25/112012
 562 //                                              TClonesArray* trArray=0;
 563 //                                              TParticle* tempParticle=0;
 564  //                                           TVector3 DecayMomentumK(0,0,0);
 565 //loop on secondaries
 566              for (Int_t k=firstD;k<=lastD;k++) {
 567               if ( k > 0 )    {
 568              TParticle*    daughter1=stack->Particle(k);   // 27/8
 569              Float_t dcode = daughter1->GetPdgCode();
 570
 571 //     mother momentum trackrefs    and QtMC     // 17/9/2010,  test Feb 2011
 572                                                 if (mcEvent->GetParticleAndTR(iMc, tempParticle, trArray) != -1) {
 573                                                 AliTrackReference* MCtrackReference = static_cast<AliTrackReference*>(trArray->Last());
 574                         DecayMomentumK.SetXYZ(MCtrackReference->Px(), MCtrackReference->Py(), MCtrackReference->Pz());
 575                               }
 576                                                 const TVector3 MCP3d(daughter1->Px(), daughter1->Py(), daughter1->Pz()); //MC daughter's momentum
 577                                                  MCQt = MCP3d.Perp(DecayMomentumK); //MC daughter's transverse momentum in mother's frame
 578 //
 579                   Double_t MCKinkAngle = TMath::ASin(MCQt/daughter1->P() );
 580                 Double_t  MCKinkAngle2= TMath::RadToDeg() * MCKinkAngle; // in degrees
 581              //    fmaxAngMomKmu->Fill(particle->P() , MCKinkAngle2);// MC
 582 //
 583              mcProcess=daughter1->GetUniqueID();
 584                   radiusD=daughter1->R();
 585 // secondary vertex
 586                //   Double_t hVx=daughter1->Vx();
 587               //  Double_t hVy=daughter1->Vy();
 588                   Double_t hVz=daughter1->Vz();
 589
 590            LengthK = TMath::Sqrt( radiusD*radiusD  + ( mVz-hVz) * (mVz-hVz) );  //   19/7/2010 mss
 591
 592 //          lengthK = TMath::Sqrt( (mVx -hVx)*( mVx    -hVx)  + ( mVy    -hVy)* (mVy    -hVy ) + ( mVz   -hVz ) * (mVz     -hVz) );
 593             lenYuri  = (TMath::Abs( mVz-hVz))* (TMath::Sqrt( 1.+ ( ptK*ptK)/ (particle->Pz() * particle->Pz()) )) ;
 594
 595        if(mcProcess==13) {
 596     mcProc13++;
 597
 598       if(mcProc13==1)     flifeInt  ->Fill(  (lengthKMC*0.493667  /particle->P()));  // 19/7
 599
 600    //            if( (charg==1)&&(mcProc13==1 )) fradIntKP->Fill(daughter1->R());
 601
 602      //            if( ( charg ==-1)&&(mcProc13==1))  fradIntKM->Fill(daughter1->R());
 603   }
 604
 605
 606 //
 607      if (mcProcess==4) {
 608
 609     mcProc4++;
 610    if ( mcProc4==1)  {
 611    // 10/1/13                if( (radiusD >120.)&&( radiusD< 210.) )  {
 612           flifeYuri ->Fill(  (lenYuri  *0.493667  /particle->P()));  // 19/7
 613                   flifetiMCK->Fill(LengthK);
 614                   flenYuri  ->Fill(lenYuri);
 615 // 10/1/13                    }
 616
 617   //                    flengthMCK->Fill(lengthK);  //
 618   //      flifetime ->Fill(  (lengthK*0.493667  /particle->P()));  // 19/7
 619                       flengthMCK->Fill(lengthKMC);  //
 620         flifetime ->Fill(  (lengthKMC*0.493667  /particle->P()));  // 19/7
 621             fradPtRapMC->Fill( radiusD, 1./ptK, rapidiKMC);  // systematics 26/8
 622   }
 623
 624
 625 //
 626     if (  ( ( code==321 )&&( dcode ==-13  ))||( ( code ==-321)&&(dcode== 13) ) || ( ( code==321 )&&( dcode ==-11  )) || ( (code ==-321)&&(dcode== 11))) {
 627                       flifetim2 ->Fill(  (lengthKMC*0.493667  /particle->P()));
 628            //   8/2/2013allgh radius    if( (radiusD >120.)&&( radiusD< 210.) )
 629            //   14/2/13 if( (radiusD >130.)&&( radiusD< 200.) )
 630            if( (radiusD >fKinkRadLow )&&( radiusD< fKinkRadUp) )
 631                  fmaxAngMomKmu->Fill(particle->P() , MCKinkAngle2);// MC
 632               }
 633
 634       if (( (TMath::Abs(code)==321 )&&(TMath::Abs(dcode)  ==211  ))&& ( mcProc4<2)) flifetim2->Fill( lengthKMC *0.493667 /particle->P()) ;//19/7
 635
 636      // test feb 2013                    if ((TMath::Abs(hVz)<0.5) || (TMath::Abs(hVz )>225)) continue;
 637 ///   inside radius region ----------------------------------------------
 638                        if(MCKinkAngle2 < 2.) continue;  // as in ESD
 639                         //          ======  8/2/13 if (((daughter1->R())>120)&&((daughter1->R())<210)&& (MCQt>0.120)  ){
 640                         if (((daughter1->R())> fKinkRadLow )&&((daughter1->R())< fKinkRadUp )&& (MCQt>fLowQt)  ){
 641
 642         if ( ( code==321 )&&( dcode ==-13  ))   {
 643             fradPtRapDC->Fill( radiusD, 1./ptK, rapidiKMC);  // systematics 26/8
 644                   fradMC->Fill(daughter1->R());
 645             fQtKMuMC ->Fill(MCQt );//to muon
 646             fgenPtEtR->Fill( ptK );//to muon
 647             fgenPtEtRP->Fill( ptK );//to muon
 648          fMCEtaKaon  ->Fill(rapidiKMC );//to muon
 649          fSignPtEtaMC->Fill(ptK,rapidiKMC );//to muon
 650          fSignPtMC->Fill(ptK);//to muon
 651                  flifetiMCK->Fill(   lenYuri*0.4933667/particle->P()  );
 652                     //  flifetiMCK->Fill(   LengthK*0.4933667/   ptK
 653     } //  positive kaon
 654         if (  ( code ==-321)&&(dcode== 13)){
 655         fgenPtEtR->Fill(  ptK   );//to muon
 656             fQtKMuMC ->Fill(MCQt );//to muon
 657                fgenPtEtRN->Fill(ptK);  //
 658         fSignPtEtaMC->Fill(chargPt,rapidiKMC );//to muon
 659         fMCEtaKaon  ->Fill(rapidiKMC );//to muon
 660             fradPtRapDC->Fill( radiusD, 1./ptK, rapidiKMC);  // systematics 26/8
 661                   fradMC->Fill(daughter1->R());
 662          fSignPtMC->Fill(chargPt);//to muon
 663                           flifetiMCK->Fill(   lenYuri*0.4933667/particle->P() ) ;
 664           //           flifetiMCK->Fill(   LengthK*0.4933667/  ptK   ) ;
 665
 666     } //  negative code
 667       //  if (( ( code==321 )&& ( dcode ==211  ))|| (( code == -321 )&& ( dcode ==-211))) fgenPtEtR->Fill( ptK );//to pion
 668         if ( ( code==321 )&&( dcode ==-11  ))   {
 669             fQtKElMC ->Fill(MCQt );//to muon
 670             fgenPtEtR->Fill( ptK );//to electron
 671             fgenPtEtRP->Fill( ptK );//to muon
 672          fMCEtaKaon  ->Fill(rapidiKMC );//to electron
 673             fradPtRapDC->Fill( radiusD, 1./ptK, rapidiKMC);  // systematics 26/8
 674                   fradMC->Fill(daughter1->R());
 675          fSignPtEtaMC->Fill(ptK,rapidiKMC );//to electron
 676          fSignPtMC->Fill(ptK);
 677                       flifetiMCK->Fill(   lenYuri*0.4933667/particle->P()  );
 678         //             flifetiMCK->Fill(   LengthK*0.4933667/ ptK      );
 679
 680     } //  positive kaon
 681         if (  ( code ==-321)&&(dcode== 11)){
 682         fgenPtEtR->Fill(   ptK  );//to electron
 683             fQtKElMC ->Fill(MCQt );//to muon
 684                fgenPtEtRN->Fill(ptK);  //
 685         fSignPtEtaMC->Fill(chargPt,rapidiKMC  );//to electron
 686         fMCEtaKaon  ->Fill(rapidiKMC  );//to electron
 687             fradPtRapDC->Fill( radiusD, 1./ptK, rapidiKMC);  // systematics 26/8
 688                   fradMC->Fill(daughter1->R());
 689          fSignPtMC->Fill(chargPt);//to electron
 690                               flifetiMCK->Fill(   lenYuri*0.4933667/particle->P()  );
 691       //               flifetiMCK->Fill(   LengthK*0.4933667/   ptK         );
 692
 693     } //  negative code
 694
 695         if (( ( code==321 )&& ( dcode ==211  ))|| (( code == -321 )&& ( dcode ==-211)))    nMCKpi++ ;
 696         if (( ( code==321 )&& ( dcode ==211  ))|| (( code == -321 )&& ( dcode ==-211)))   {
 697                  if ( nMCKpi > 0) {
 698               MCQt3[nMCKpi-1] = MCQt ;// k to pipipi
 699 }
 700                        }
 701     nMCKinkKs++;
 702        }
 703
 704                     }//    decay
 705          } // positive k
 706        }//  daughters
 707
 708
 709
 710                  if ( code > 0) {
 711               if( nMCKpi == 1) fgenPtEtR->Fill(ptK);  //  k to pipi
 712               if( nMCKpi == 1) fgenPtEtRP->Fill(ptK);  //  k to pipi
 713               if( nMCKpi == 1) fSignPtEtaMC->Fill(ptK,rapidiKMC );  //  k to pipi
 714               if( nMCKpi == 1) fSignPtMC->Fill(ptK);  //  k to pipi
 715               if( nMCKpi == 1) fMCEtaKaon->Fill(rapidiKMC );  //  k to pipi
 716                   if(nMCKpi==1) fradMC->Fill(radiusD       );
 717                   if(nMCKpi==1) fQtKPiMC->Fill( MCQt   );
 718                 if (nMCKpi==1)  fradPtRapDC->Fill( radiusD, 1./ptK, rapidiKMC);  // systematics 26/8
 719                   if(nMCKpi==1)     flifetiMCK->Fill(   lenYuri*0.4933667/particle->P()  );
 720              //   if(nMCKpi==1)     flifetiMCK->Fill(   LengthK*0.4933667/      ptK      );
 721    }   //positive kaon
 722      //
 723            if ( code < 0) {
 724
 725               if( nMCKpi == 1) fgenPtEtR->Fill(   ptK );  //  k to pipi
 726               if( nMCKpi == 1) fgenPtEtRN->Fill(ptK);  //  k to pipi
 727               if( nMCKpi == 1) fSignPtEtaMC->Fill(chargPt,rapidiKMC );  //  k to pipi
 728               if( nMCKpi == 1) fSignPtMC->Fill(chargPt);  //  k to pipi
 729               if( nMCKpi == 1) fMCEtaKaon->Fill(rapidiKMC );  //  k to pipi
 730                   if(nMCKpi==1) fradMC->Fill(radiusD       );
 731                   if(nMCKpi==1) fQtKPiMC->Fill( MCQt   );
 732                              if( nMCKpi== 1) fradPtRapDC->Fill( radiusD, 1./ptK, rapidiKMC);  // systematics 26/8
 733                  if(nMCKpi==1)     flifetiMCK->Fill(   lenYuri*0.4933667/particle->P()  );
 734          //       if(nMCKpi==1)     flifetiMCK->Fill(   LengthK*0.4933667/      ptK      );
 735
 736    }   //negative K
 737
 738       }   /// kaons  loop
 739
 740
 741
 742     nMCTracks++;
 743   }// end of mc particle
 744
 745 //  Phys sel    2012 EFF calculation
 746 Bool_t isSelected =
 747 ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kMB;
 748
 749          if ( isSelected ==kFALSE) return;   //  24/6/11 apo MF
 750 //
 751        fMultiplMC->Fill(nPrim);
 752 //=======================================================================================
 753 //            main vertex selection
 754   const AliESDVertex *vertex=GetEventVertex(esd);
 755     if(!vertex) return;
 756 ///
 757 //       fMultiplMC->Fill(nPrim);
 758 //
 759 /*  / apo Alexander  Feb 2012
 760     AliESDpid *fESDpid = new AliESDpid();
 761     if (!fESDpid) fESDpid =
 762   ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
 763 */       ///========================================================================
 764 //               apo Eftihi
 765                   if(!fPIDResponse) {
 766     AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 767     AliInputEventHandler* inputHandler =
 768 (AliInputEventHandler*)(man->GetInputEventHandler());
 769     fPIDResponse = inputHandler->GetPIDResponse();
 770   }
 771
 772
 773     Double_t vpos[3];
 774      vertex->GetXYZ(vpos);
 775     fZpr->Fill(vpos[2]);
 776     if ( TMath::Abs( vpos[2])  > 10. ) return;  ///  it is applied on ESD and generation
 777
 778   Double_t vtrack[3], ptrack[3];
 779
 780
 781  // Int_t nESDTracK = 0;
 782
 783    Int_t nGoodTracks =  esd->GetNumberOfTracks();
 784 //
 785     fESDMult->Fill(nGoodTracks);
 786
 787 //
 788       Double_t nsigmall = 100.0;
 789        Double_t nsigma = 100.0;
 790        Double_t nsigmaPion =-100.0;
 791        Double_t nsigmaPi=-100.0;
 792   //     Double_t dEdxDauMC  =   0.0;
 793
 794 //
 795 for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 796
 797     AliESDtrack* trackD = esd->GetTrack(iTrack);
 798     if (!trackD) {
 799       Printf("ERROR: Could not receive track %d", iTrack);
 800       continue;
 801     }
 802
 803          //       Int_t indexKinkDau=trackD->GetKinkIndex(0);
 804 // daughter kink
 805           nsigmaPion     = (fPIDResponse->NumberOfSigmasTPC(trackD  , AliPID::kPion));// 26/10 eftihis
 806  //   nsigmaPion= (fESDpid->NumberOfSigmasTPC(trackD,AliPID::kPion));
 807  //   22/11/12 if((indexKinkDau >0)&& (nsigmaPion>1.2)) fTPCSgnlKinkDau->Fill(trackD->P(), (trackD->GetTPCsignal()  ) ) ;  //  daughter kink
 808  //if((indexKinkDau >0)) fTPCSgnlKinkDau->Fill(trackD->P(), (trackD->GetTPCsignal()  ) ) ;  //  daughter kink
 809 // if((indexKinkDau >0))    dEdxDauMC   = trackD->GetTPCsignal()     ;  //  daughter kink
 810    }
 811
 812 // loop on ESD tracks
 813
 814 //
 815    for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
 816
 817     AliESDtrack* track = esd->GetTrack(iTracks);
 818     if (!track) {
 819       Printf("ERROR: Could not receive track %d", iTracks);
 820       continue;
 821     }
 822
 823
 824
 825 //    sigmas
 826              nsigmall  = (fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon));
 827 //    nsigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon));
 828  //   nsigmall = (fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon));
 829        nsigma  = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon));
 830
 831 //==================================
 832   Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 833   Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 834   if (track->GetTPCNclsF()>0) {
 835     ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC / track->GetTPCNclsF();
 836     fRatioCrossedRows->Fill(ratioCrossedRowsOverFindableClustersTPC);
 837    }
 838
 839     fHistPt->Fill(track->Pt());
 840
 841
 842   Double_t tpcNCl = track->GetTPCclusters(0);  //
 843   //Int_t tpcNCl = track->GetTPCclusters(0);  //
 844      Double_t tpcSign = track->GetSign();  //
 845
 846   Int_t label = track->GetLabel();
 847
 848     label = TMath::Abs(label);
 849
 850      if(label > mcEvent->GetNumberOfTracks()) continue; //
 851
 852     TParticle * part = stack->Particle(label);
 853     if (!part) continue;
 854 // loop only on Primary tracks
 855      if (label > nPrim) continue; /// primary tracks only   ,21/3/10  EFF study
 856
 857      //    pt cut
 858       if ( (track->Pt())<.200)continue;   //12/2/2012     -------------------------------------------
 859
 860     UInt_t status=track->GetStatus();
 861
 862      if((status&AliESDtrack::kITSrefit)==0) continue;
 863     if((status&AliESDtrack::kTPCrefit)==0) continue;   //6 feb
 864      //  if((track->GetTPCchi2()/track->GetTPCclusters(0))>3.8) continue;
 865      //if((track->GetTPCchi2()/track->GetTPCclusters(0))>4.5) continue;   // test 10/3/2012
 866      if((track->GetTPCchi2()/track->GetTPCclusters(0))>4.0) continue;   // test 10/3/2012
 867
 868       Double_t extCovPos[15];
 869       track->GetExternalCovariance(extCovPos);
 870 //
 871
 872
 873     track->GetXYZ(vtrack);
 874  fXvYv->Fill(vtrack[0],vtrack[1]);
 875  fZvYv->Fill(vtrack[0],vtrack[2]);
 876  fZvXv->Fill(vtrack[1],vtrack[2]);
 877
 878 // track momentum
 879      track->GetPxPyPz(ptrack);
 880
 881     TVector3 trackMom(ptrack[0],ptrack[1],ptrack[2]);
 882
 883           Double_t   etracK= TMath::Sqrt(trackMom.Mag()*trackMom.Mag() + 0.493677 *0.493677  );
 884          Double_t rapiditK = 0.5 * (TMath::Log(  (etracK + ptrack[2]  ) / ( etracK - ptrack[2])  ))  ;
 885     Double_t trackEta=trackMom.Eta();
 886    // Double_t trMoment= trackMom.Mag();
 887     Double_t trackPt = track->Pt();
 888
 889
 890     Float_t bpos[2];
 891     Float_t bCovpos[3];
 892     track->GetImpactParameters(bpos,bCovpos);
 893
 894     if (bCovpos[0]<=0 || bCovpos[2]<=0) {
 895      Printf("Estimated b resolution lower or equal zero!");
 896      bCovpos[0]=0; bCovpos[2]=0;
 897     }
 898
 899     Float_t dcaToVertexXYpos = bpos[0];
 900 //    Float_t dcaToVertexZpos = bpos[1];
 901
 902     //fRpr->Fill(dcaToVertexZpos);
 903     fRpr->Fill(dcaToVertexXYpos);
 904
 905           //if((TMath::Abs(dcaToVertexXYpos) >0.3)||(TMath::Abs(dcaToVertexZpos)>2.5)) continue;  //   22/7/11
 906    //       if((TMath::Abs(dcaToVertexXYpos) >0.24)||(TMath::Abs(dcaToVertexZpos)>2.5)) continue;  //   12/2/13
 907                      if (!fMaxDCAtoVtxCut->AcceptTrack(track)) continue;
 908
 909
 910  //  cut on eta
 911  //       if(  (TMath::Abs(trackEta )) > 0.9 ) continue;
 912         if(  (TMath::Abs(rapiditK  )) > 0.7 ) continue; ////   rapid K, Feb 20
 913     fHistPtESD->Fill(track->Pt());
 914
 915    // Add Kink analysis
 916
 917                 Int_t indexKinkPos=track->GetKinkIndex(0);
 918 //  loop on mother kinks
 919                 if(indexKinkPos<0){
 920                fPtKink->Fill(track->Pt()); /// pt from track
 921
 922                     fRatioCrossedRowsKink->Fill(ratioCrossedRowsOverFindableClustersTPC);
 923
 924         // select kink class
 925
 926           AliESDkink *kink=esd->GetKink(TMath::Abs(indexKinkPos)-1);
 927 //
 928
 929 // DCA kink
 930           Double_t  Dist2 = kink->GetDistance();
 931   //        fDCAkink->Fill( Dist2   );
 932 //
 933           Int_t eSDfLabel1=kink->GetLabel(0);
 934           TParticle *particle1= stack->Particle(TMath::Abs(eSDfLabel1));
 935           Int_t code1= particle1->GetPdgCode();
 936  //         Int_t mcProcssMo= particle1->GetUniqueID();
 937
 938           Int_t eSDfLabeld=kink->GetLabel(1);
 939           TParticle *particled= stack->Particle(TMath::Abs(eSDfLabeld));
 940           Int_t dcode1= particled->GetPdgCode();
 941           Int_t mcProcssDa= particled->GetUniqueID();
 942 //
 943 //    loop on MC daugthres for 3Pi    24/9/2010
 944        Int_t nESDKpi =0;
 945        if(mcProcssDa==4) {
 946             Int_t firstDa=particle1->GetFirstDaughter();
 947             Int_t lastDa=particle1->GetLastDaughter();
 948
 949              if( (lastDa<=0) || (firstDa<=0)  ) continue;
 950
 951                      if ( lastDa > mcEvent->GetNumberOfTracks() ) continue;
 952                      if (firstDa > mcEvent->GetNumberOfTracks() ) continue;
 953 //loop on secondaries
 954              for (Int_t kk=firstDa;kk<=lastDa;kk++) {
 955               if ( kk > 0 )    {
 956              TParticle*    daughter3=stack->Particle(kk);   // 24/9
 957              Float_t dcode3= daughter3->GetPdgCode();
 958         if (( ( code1==321 )&& ( dcode3==211  ))|| (( code1 == -321 )&& ( dcode3==-211)))    nESDKpi++ ;
 959              }
 960              }
 961              }
 962 //---------------------------edw telos  9/2010
 963                   Double_t hVzdau=particled->Vz();
 964
 965            const TVector3 motherMfromKink(kink->GetMotherP());
 966            const TVector3 daughterMKink(kink->GetDaughterP());
 967            Float_t qT=kink->GetQt();
 968     //   Float_t motherPt=motherMfromKink.Pt();
 969 // Kink  mother momentum
 970 //     Double_t trMomTPCKink=motherMfromKink.Mag();
 971 // TPC mother momentun
 972   //   Double_t trMomTPC=track->GetTPCmomentum();
 973        //     Float_t etaMother=motherMfromKink.Eta();
 974
 975
 976            fHistQtAll->Fill(qT) ;  //  Qt   distr
 977
 978            const TVector3 vposKink(kink->GetPosition());
 979           fPosiKink ->Fill( vposKink[0], vposKink[1]  );
 980
 981                    Double_t dxKink = vpos[0]-vposKink[0], dyKink=vpos[1]-vposKink[1], dzKink=vpos[2]-vposKink[2] ;
 982   //   Double_t  dzKink=vpos[2]-vposKink[2] ;    ///  ??
 983    Double_t lifeKink= TMath::Sqrt( dxKink*dxKink + dyKink*dyKink + dzKink*dzKink ) ;
 984
 985              Double_t  tanLamda   = track-> GetTgl()    ;//  ??
 986                                   if (tanLamda ==0 )  continue;//   ??
 987         Double_t lenHelx = (TMath::Abs(dzKink    )   ) *(TMath::Sqrt( 1. + tanLamda *tanLamda  ) ) / ( TMath::Abs( tanLamda))  ;// ??
 988
 989
 990
 991
 992            Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
 993
 994 //       fake kinks, small Qt and small kink angle
 995     if(( kinkAngle<1.))  fHistQt1  ->Fill(qT) ;  //  Qt   distr
 996 //
 997     if(       ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==321))) fFakeKPi->Fill(track->Pt());
 998     if(       ( (TMath::Abs(code1)==211)&&(TMath::Abs(dcode1)==211))) fFakepipi->Fill(track->Pt());
 999 //
1000
1001 //remove the double taracks
1002            if( (kinkAngle<2.)  ) continue;    // test 15/7 2010 , it removes 3  from   10000 good Kaons
1003          //  BG  ?????==============
1004               if ( TMath::Abs(vposKink[2]) >  225. ) continue ;
1005               if ( TMath::Abs(vposKink[2]) <  0.5 ) continue ;
1006
1007                     fPtCut1   ->Fill(trackPt );
1008                   fChi2NclTPC->Fill( (track->GetTPCchi2() ) ,  tpcNCl );
1009
1010              Float_t signPt= tpcSign*trackPt;
1011 //
1012            if  ( (TMath::Abs(code1)==211)&&(TMath::Abs(dcode1)==13))
1013            fRadiusPtPion->Fill( kink->GetR(), track->Pt()); //
1014
1015          if(       ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==13))||
1016            ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==11))  ||
1017            ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==211))  ) {
1018            fRadiusPtKaon->Fill( kink->GetR(), track->Pt()); //
1019               }
1020 //
1021          // ======8/1/13 if((kink->GetR()>120.)&&(kink->GetR()<210.)&&(TMath::Abs(rapiditK)<0.7)&&(label<nPrim)) {
1022          if((kink->GetR()> fKinkRadLow )&&(kink->GetR()< fKinkRadUp )&&(TMath::Abs(rapiditK)<0.7)&&(label<nPrim)) {
1023     if(       ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==13))) fQtKMu->Fill(qT);
1024     if     ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==11))   fQtKEl->Fill(qT);
1025     if     ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==211))   fQtKPi->Fill(qT);
1026     if  (( nESDKpi>1) &&    ( ((code1)==321)&&((dcode1)==211)) )   fQtK3PiP->Fill(qT);
1027     if  (( nESDKpi>1) &&    ( ((code1)==-321)&&((dcode1)==-211)) )   fQtK3PiM->Fill(qT);
1028          if(       ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==13))||
1029            ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==11))  ||
1030            ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==211))  ) {
1031          if(qT>fLowQt)        fHistPtKPDG->Fill(track->Pt());  // ALL KAONS (pdg) inside ESD  kink sample
1032         if(qT>fLowQt )    {
1033       if(code1>0.)  fHiPtKPDGP->Fill(trackPt             ); // 26/feb  // ALL KAONS (pdg) inside ESD  kink sample
1034       if(code1<0.)  fHiPtKPDGN->Fill(      trackPt       ); // 26/feb  // ALL KAONS (pdg) inside ESD  kink sample
1035                    }
1036             fHistEta->Fill(trackEta) ;  //   Eta distr of PDG kink ESD  kaons
1037             frapidESDK->Fill(rapiditK) ;  //18/feb rapiddistr of PDG kink ESD  kaons
1038       if( qT > fLowQt  )  fHistQt2->Fill(qT);  // PDG ESD kaons
1039            fRadiusPt->Fill( kink->GetR(), track->Pt()); //
1040      }
1041      }
1042
1043
1044            //Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
1045            Double_t maxDecAngKmu=f1->Eval(track->P(),      0.,0.,0.);
1046            Double_t maxDecAngpimu=f2->Eval(track->P(),   0.,0.,0.);
1047 //  two dimensional plot
1048                 if(TMath::Abs(code1)==321) fAngMomK->Fill(track->P(),            kinkAngle);
1049                 if(TMath::Abs(code1)==211) fAngMomPi->Fill( track->P(), kinkAngle);
1050 //_______
1051
1052 //
1053 // invariant mass of mother track decaying to mu
1054          Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
1055          Float_t energyDaughterPi=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.139570*0.139570);
1056          Float_t energyDaughterKa=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.493677*0.493677);
1057 //       Float_t energyDaughterPr=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.938658*0.938658);
1058          Float_t p1XM= motherMfromKink.Px();
1059          Float_t p1YM= motherMfromKink.Py();
1060          Float_t p1ZM= motherMfromKink.Pz();
1061          Float_t p2XM= daughterMKink.Px();
1062          Float_t p2YM= daughterMKink.Py();
1063          Float_t p2ZM= daughterMKink.Pz();
1064          Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
1065          Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
1066          Double_t invariantMassKpi= TMath::Sqrt((energyDaughterPi+p3Daughter)*(energyDaughterPi+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
1067          Double_t invariantMassKK = TMath::Sqrt((energyDaughterKa+p3Daughter)*(energyDaughterKa+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
1068          fInvMassMuNuAll ->Fill(invariantMassKmu);
1069          fInvMassMuNuPtAll ->Fill(invariantMassKmu,  trackPt);
1070            //    20/4 testRadiusPt->Fill( kink->GetR(), track->Pt()); //
1071
1072
1073                if (qT> fLowQt )     fSignPtNcl->Fill( signPt  ,   tpcNCl   );
1074
1075 //
1076     //  if((qT>0.12)&&((kink->GetR()>120.)&&(kink->GetR()<210.))&&(TMath::Abs(rapiditK )<0.7)) {
1077     if((qT> fLowQt )&&((kink->GetR()> fKinkRadLow )&&(kink->GetR()< fKinkRadUp ))&&(TMath::Abs(rapiditK )<0.7)) {
1078          fM1kaon->Fill(invariantMassKmu);
1079          fMinvPi->Fill(invariantMassKpi);
1080          fMinvKa->Fill(invariantMassKK);
1081          fRadiusNcl->Fill( (kink->GetR()) ,(track->GetTPCclusters(0)  ) ) ;
1082                         }
1083
1084 //
1085          //  if ( tpcNCl<30  ) continue;
1086       //   test for systematics  , march 13   if ( tpcNCl<50. ) continue;
1087             if ( tpcNCl<fLowCluster  ) continue;
1088       //   if (  ( tpcNCl<20. )|| ( tpcNCl > 100.)  )  continue;//   test   , den edwse kati shmantiko
1089
1090                              //if( ( ( track->GetTPCclusters(0) ) / (kink->GetR() ) ) > 0.63 ) continue;
1091       Double_t tpcNClHigh = -51.67+ (11./12.)  *( kink->GetR() ) ;
1092          if (tpcNCl >  tpcNClHigh  )   fcodeH->Fill(TMath::Abs(code1), TMath::Abs(dcode1));
1093         if (  tpcNCl  >tpcNClHigh)   fZkinkZDau->Fill( vposKink[2],hVzdau              );
1094            if (tpcNCl > tpcNClHigh) fRadiusPtFake->Fill( kink->GetR(), track->Pt()); //
1095
1096      Double_t tpcNClMin  = -85.5 + (65./95.)  *( kink->GetR() ) ;
1097                // if ( tpcNClMin < tpcNCl ) continue;
1098                if ( tpcNCl > tpcNClHigh) continue;
1099                if ( tpcNCl < tpcNClMin ) continue;
1100  //
1101
1102 //  kaon selection from kinks
1103
1104    //=====  8/2/13 if((kinkAngle>maxDecAngpimu)&&(qT>0.12)&&(qT<0.30)&&((kink->GetR()>120.)&&(kink->GetR()<210.))&&(TMath::Abs(rapiditK )<0.7)&&(invariantMassKmu<0.8)) {
1105    if((kinkAngle>maxDecAngpimu)&&(qT>fLowQt )&&(qT<0.30)&&((kink->GetR()> fKinkRadLow )&&(kink->GetR()< fKinkRadUp ))&&(TMath::Abs(rapiditK )<0.7)&&(invariantMassKmu<0.8)) {
1106 // 29092010     if((kinkAngle>maxDecAngpimu)&&(qT>0.120)&&(qT<0.25)&&((kink->GetR()>120.)&&(kink->GetR()<210.))&&(TMath::Abs(rapiditK )<0.7)&&(invariantMassKmu<0.6)) {
1107 //  if((kinkAngle>maxDecAngpimu)&&(qT>0.04)&&(qT<0.30)&&((kink->GetR()>133.)&&(kink->GetR()<179.))&&(TMath::Abs(rapiditK )<0.5)&&(invariantMassKmu<0.6)) {
1108
1109                 fAngMomKKinks->Fill(track->P(), kinkAngle);
1110             fPtCut2   ->Fill(trackPt );
1111
1112                  if ( (kinkAngle>maxDecAngKmu*0.98)&& ( track->P() > 1.2 ) ) continue; // maximum angle s
1113                  if ( (kinkAngle<maxDecAngpimu*1.20)  )  continue; // maximum angle s
1114
1115                 fPtCut3   ->Fill(trackPt );
1116      //fTPCSgnlPa->Fill(track->P(),track->GetTPCsignal());
1117      fTPCSgnlPa->Fill(track->GetInnerParam()->GetP(),track->GetTPCsignal());
1118                 //    if( nsigma > 3.5 )      fcode2->Fill(TMath::Abs(code1), TMath::Abs(dcode1));
1119                   if(nsigma  > 3.5) continue;  // 1/11/12
1120       // test 16/2/13           if(nsigma  > 3.5) continue;  // 15/2/13
1121                //  if(nsigma  > 4.0) continue; // test 17/2/2011  4% or more ? bg?
1122 //
1123    fTPCSgnlP->Fill(track->GetInnerParam()->GetP(), (track->GetTPCsignal()  ) ) ;
1124
1125          fInvMassMuNuPt ->Fill(invariantMassKmu,  trackPt);
1126 //  loop on kink daughters inside  mother's loop
1127           Int_t ESDLabelM   =  0. ;
1128           Int_t ESDLabelD   =  0. ;
1129        Double_t dEdxDauMC  =   0.0;
1130         Double_t raDAU=0.;
1131         Int_t Ikink =0;
1132         Int_t IRkink =0;
1133 for (Int_t jTrack = 0; jTrack < esd->GetNumberOfTracks(); jTrack++) {
1134
1135     AliESDtrack* trackDau = esd->GetTrack(jTrack);
1136     if (!trackDau) {
1137       Printf("ERROR: Could not receive track %d", jTrack);
1138       continue;
1139     }
1140                 Int_t indexKinkDAU =trackDau->GetKinkIndex(0);
1141                      if (indexKinkDAU <0  ){
1142           AliESDkink *kinkDau=esd->GetKink(TMath::Abs(indexKinkDAU)-1);
1143                raDAU= kinkDau->GetR();
1144            ESDLabelM=kinkDau->GetLabel(0);   //  mothers's label
1145                 ESDLabelM = TMath::Abs(ESDLabelM);
1146            ESDLabelD=kinkDau->GetLabel(1);   //  Daughter'slabel
1147                 ESDLabelD = TMath::Abs(ESDLabelD);
1148                      if ( kink->GetR() == kinkDau->GetR() ) IRkink++;
1149                      if ( ESDLabelM == label ) Ikink++  ;
1150    }
1151   //           if (( ESDLabelM== eSDfLabel1))   {
1152              if (   (Ikink >0)  && (IRkink>0 )       )   {
1153 // daughter kink
1154      //if(indexKinkDAU >0)     nsigmaPi     = (fPIDResponse->NumberOfSigmasTPC(trackDau,AliPID::kPion));// 26/10 eftihis
1155      if(indexKinkDAU >0)     nsigmaPi     = (fPIDResponse->NumberOfSigmasTPC(trackDau,AliPID::kKaon));// 26/10 eftihis
1156  //   nsigmaPion= (fESDpid->NumberOfSigmasTPC(trackD,AliPID::kPion));
1157  //   22/11/12 if((indexKinkDau >0)&& (nsigmaPion>1.2)) fTPCSgnlKinkDau->Fill(trackD->P(), (trackD->GetTPCsignal()  ) ) ;  //  daughter kink
1158  //if((indexKinkDau >0)) fTPCSgnlKinkDau->Fill(trackD->P(), (trackD->GetTPCsignal()  ) ) ;  //  daughter kink
1159  if((indexKinkDAU      >0))    dEdxDauMC   = trackDau->GetTPCsignal()     ;  //  daughter kink
1160    }
1161    }
1162 // end internal loop for kink daughters
1163
1164       //   fTPCSgnlP->Fill(trMomTPC  , (track->GetTPCsignal()  ) ) ;
1165       //   fTPCSgnlPtpc->Fill(trMomTPC  , (track->GetTPCsignal()  ) ) ;
1166          //fMothKinkMomSgnl ->Fill(trMomTPCKink  , (track->GetTPCsignal()  ) ) ;
1167      //    fMothKinkMomSgnl ->Fill( dEdxDauMC     , (track->GetTPCsignal()  ) ) ;
1168          // fTPCMomNSgnl->Fill(trMomTPC ,pidResponse->NumberOfSigmasTPC(track, AliPID::kKaon)  );
1169          fNSigmTPC   ->Fill(nsigmall );
1170 //  daughter selection
1171   //fTPCSgnlKinkDau->Fill( daughterMKink.Mag() ,     dEdxDauMC  ) ;  //  daughter kink
1172            // if(  nsigmaPion > 1.0 )        fcode4->Fill(TMath::Abs(code1), TMath::Abs(dcode1));
1173        //    if( dEdxDauMC > 1.5 *(track->GetTPCsignal()   ) )       fcode4->Fill(TMath::Abs(code1), TMath::Abs(dcode1));
1174 //
1175          //fTPCMomNSgnl->Fill(trMomTPC ,nsigmall );
1176          fTPCMomNSgnl->Fill(track->GetInnerParam()->GetP() ,nsigmall );
1177 //
1178
1179           fRadNclcln->Fill( (kink->GetR()) ,(track->GetTPCclusters(0)  ) ) ;
1180            fRadiusPtcln->Fill( kink->GetR(), trackPt); //
1181
1182              fAngMomKC->Fill(track->P()           , kinkAngle);
1183                                     fHistPtKaon->Fill(track->Pt()         );   //all PID kink-kaon
1184         if( code1>0.)     fHistPtKaoP->Fill(track->Pt()         );   //all PID kink-kaon
1185         if( code1<0.)    fHistPtKaoN->Fill(track->Pt()         );   //all PID kink-kaon
1186 // systematics
1187                  fradPtRapESD->Fill(kink->GetR(), 1./ track->Pt(), rapiditK);
1188 //
1189                 fHistEtaK->Fill(rapiditK );
1190                  frad->Fill( kink->GetR() );
1191                flenHelx->Fill( lenHelx   );  //??
1192                 flifeKink ->Fill(lifeKink      );//??
1193                         fLHelESDK ->Fill(  ( lenHelx /track->P() )*0.493677);// for all 'PID' kaons  31/7/11// ??
1194                       flifTiESDK->Fill(  ( lifeKink    /track->P() )*0.493677);  //  ??
1195
1196
1197
1198                   fSignPtNcl->Fill( signPt  ,   tpcNCl   );
1199                   fSignPtEta->Fill(signPt  , rapiditK  );
1200                   fEtaNcl->Fill(rapiditK, tpcNCl    );
1201                   fSignPt->Fill( signPt   );
1202        fRatChi2Ncl-> Fill((track->GetTPCchi2()/track->GetTPCclusters(0) )) ;
1203     fdcatoVxXY->Fill(dcaToVertexXYpos);
1204
1205 //  kaons from k to mun and k to pipi and to e decay
1206          if(       ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==13))||
1207            ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==11))  ||
1208            ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==211))  ) {
1209
1210          if ( label<=nPrim ) fPtPrKink->Fill(track->Pt());
1211
1212                                   //         flifetim2  ->Fill(  ( lenHelx /track->P() )*0.493667); // to compare with fLHelESDK
1213                       flifTiESDK->Fill(  ( lifeKink    /track->P() )*0.493667);
1214
1215
1216
1217               fKinkKaon->Fill(track->Pt());
1218           fDCAkink->Fill( Dist2   );
1219           fPosiKinkK->Fill( vposKink[0], vposKink[1]  );
1220           fPosiKinKXZ->Fill( vposKink[0], vposKink[2]  );
1221           fPosiKinKYZ->Fill( vposKink[1], vposKink[2]  );
1222         if( code1>0.)     fkinkKaonP->Fill(trackPt);  //                  kPtPID kink-kaon
1223         if( code1<0.)    fkinkKaonN->Fill(trackPt);    //    PID kink-kaon
1224 //     daughters
1225            if((((nsigmaPi) > 0.)&& ( dEdxDauMC > 70.  ) ))       fcode4->Fill(TMath::Abs(code1), TMath::Abs(dcode1));
1226           if ( TMath::Abs(dEdxDauMC - track->GetTPCsignal() ) <  2)  fcode2->Fill( TMath::Abs(code1), TMath::Abs(dcode1));
1227        //  fTPCSgnlPtpc->Fill(trMomTPC  , (track->GetTPCsignal()  ) ) ;
1228   fTPCSgnlPtpc   ->Fill( daughterMKink.Mag() ,     dEdxDauMC  ) ;  //  daughter kink
1229          fMothKinkMomSgnlD->Fill( dEdxDauMC     , (track->GetTPCsignal()  ) ) ;
1230                              }
1231          else {
1232               fKinkKaonBg->Fill(track->Pt());
1233          fMothKinkMomSgnl ->Fill( dEdxDauMC     , (track->GetTPCsignal()  ) ) ;
1234 //  fTPCSgnlKinkDau->Fill( daughterMKink.Mag() ,     dEdxDauMC  ) ;  //  daughter kink
1235    //        if( dEdxDauMC > 1.5 *(track->GetTPCsignal()   ) )       fcode4->Fill(TMath::Abs(code1), TMath::Abs(dcode1));
1236       //     if(( (track->P())<1. )&& ( dEdxDauMC > 1.5 *(track->GetTPCsignal()   ) ))   fcodeDau1  ->Fill(TMath::Abs(code1), TMath::Abs(dcode1));
1237            if( (( nsigmaPi) > 0. ) && ((  dEdxDauMC > 70  ) ))   fcodeDau1  ->Fill(TMath::Abs(code1), TMath::Abs(dcode1));
1238           if ( TMath::Abs(dEdxDauMC - track->GetTPCsignal() ) <  2)  fcodeDau2->Fill( TMath::Abs(code1), TMath::Abs(dcode1));
1239   fTPCSgnlKinkDau->Fill( daughterMKink.Mag() ,     dEdxDauMC  ) ;  //  daughter kink
1240 //
1241          fMinvPr->Fill(invariantMassKmu);
1242 //
1243           fDCAkinkBG->Fill( Dist2   );
1244           fPosiKinKBgXY->Fill( vposKink[0], vposKink[1]  );
1245           fPosiKinKBgZY->Fill( vposKink[2], vposKink[1]  );
1246           fPosiKinKBgZX->Fill( vposKink[2], kink->GetR() );  //  31/7/11
1247         if( code1>0.)     fKinKBGP  ->Fill(   trackPt          );   //all PID kink-kaon
1248         if( code1<0.)     fKinKBGN  ->Fill( trackPt          );   //all PID kink-kaonl
1249  fdcodeH->Fill( TMath::Abs(code1), TMath::Abs(dcode1));   // put it here,  22/10/2009
1250  //        if (eSDfLabel1==eSDfLabeld)   fcodeH->Fill(TMath::Abs(code1), TMath::Abs(dcode1));
1251    //      if (eSDfLabeld>nPrim     )   fZkinkZDau->Fill( vposKink[2],hVzdau              );
1252
1253           }   // primary and all +BG
1254
1255         }  //  kink selection
1256
1257
1258         }  //End Kink Information
1259
1260
1261   } //track loop
1262
1263 //                                       } // vx 10 cm only on  esd
1264   PostData(1, fListOfHistos);
1265
1266 }
1267
1268 //________________________________________________________________________
1269 void AliAnalysisKinkESDMC::Terminate(Option_t *)
1270 {
1271   // Draw result to the screen
1272   // Called once at the end of the query
1273
1274 }
1275
1276 const AliESDVertex* AliAnalysisKinkESDMC::GetEventVertex(const AliESDEvent* esd) const
1277   // Get the vertex from the ESD and returns it if the vertex is valid
1278
1279 {
1280   // Get the vertex
1281
1282 // 10/4  const AliESDVertex* vertex = esd->GetPrimaryVertex();
1283   const AliESDVertex* vertex = esd->GetPrimaryVertexTracks();
1284
1285   //if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
1286   if((vertex->GetStatus()==kTRUE)) return vertex;
1287   else
1288   {
1289      vertex = esd->GetPrimaryVertexSPD();
1290       if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>0)) return vertex;
1291   //   if((vertex->GetStatus()==kTRUE)) return vertex;
1292      else
1293      return 0;
1294   }
1295 }