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 //! ! ! ! !  KINK FROM HERE --------------->
 139                   Double_t gPt7Comb[48] = {
 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.0
 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 }