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