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