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