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