PWGLF/SPECTRA/Kinks/AliAnalysisKinkESDat.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  * The method is applied in pp and Pb-Pb real data.                       *
   6  *                                                                        *
   7  **************************************************************************/
   8
   9 //-----------------------------------------------------------------
  10 //                 AliAnalysisKinkESDat class
  11 //       Example of an analysis task for kink topology study
  12 //      Kaons from kink topology are 'identified' in this code
  13 //-----------------------------------------------------------------
  14
  15 #include "TChain.h"
  16 #include "TTree.h"
  17 #include "TH1F.h"
  18 #include "TH2F.h"
  19 #include "TH3F.h"
  20 #include "TH1D.h"
  21 #include "TH2D.h"
  22 #include "TParticle.h"
  23 #include <TVector3.h>
  24 #include "TF1.h"
  25
  26 #include "AliAnalysisTask.h"
  27 #include "AliAnalysisManager.h"
  28
  29 #include "AliVEvent.h"
  30 #include "AliESDEvent.h"
  31 #include "AliMCEvent.h"
  32 #include "AliAnalysisKinkESDat.h"
  33 #include "AliStack.h"
  34 #include "AliESDpid.h"
  35 #include "AliPID.h"
  36 #include "AliESDkink.h"
  37 #include "AliESDtrack.h"
  38 #include "AliPhysicsSelectionTask.h"
  39 #include "AliInputEventHandler.h"
  40 #include "AliESDInputHandler.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliVEvent.h"
  43 #include "AliESDtrackCuts.h"
  44 #include "AliPIDResponse.h"
  45 /////////#include "AliTENDERSupplies.h"
  46 ClassImp(AliAnalysisKinkESDat)
  47
  48
  49 //________________________________________________________________________
  50 AliAnalysisKinkESDat::AliAnalysisKinkESDat(const char *name)
  51   : AliAnalysisTaskSE(name), fHistPtESD(0),fHistPt(0),fHistQtAll(0),fHistQt1(0),fHistQt2(0)
  52   , fHistPtKaon(0),fHistPtKPDG(0),fHistEta(0),fHistEtaK(0),fptKMC(0),fMultiplMC(0),fESDMult(0),fgenpt(0),frad(0),
  53   fKinkKaon(0),fKinKRbn(0), fKinkKaonBg(0), fM1kaon(0),  fPtKink(0),  fptKink(0),
  54     fAngMomK(0),fAngMomPi(0), fAngMomKC(0),  fMultESDK(0), fMultMCK(0),
  55  fSignPtNcl(0), fSignPtEta(0), fEtaNcl(0), fSignPt(0), fChi2NclTPC(0), fRatChi2Ncl(0), fRadiusNcl(0), fTPCSgnlP(0),
  56    fTPCSgnlPa(0), fRpr(0),fZpr(0), fdcatoVxXY(0), fnSigmToVx(0), fKinkMothDau(0),
  57  fZvXv(0),fZvYv(0), fXvYv(0),  fHistPtKaoP(0), fHistPtKaoN(0),frapiKESD(0), flifetime(), fradLK(0),
  58     fradPtRpDt(0), fInvMuNuAll(0), fQtInvM(0),
  59          fDCAkink(0), fPosiKink(0),  fPosiKinkK(0),fPosiKinKXZ(0), fPosiKinKYZ(0), fPosiKinKBg(0), fQtMothP(0), fTPCSgnlPtpc(0),
  60        fTPCMomNSgnl(0),  fMothKinkMomSgnl(0), fNSigmTPC(0),  fTPCSgnlKinkDau(0), fPtKinkPos(0), fPtKinkNeg(0),  fRadNclCln(0),
  61        fRatioCrossedRows(0), fRatioCrossedRowsKink(0),fRadiusPt(0), fRadiusPtcln(0),  fInvMassMuNuPt(0), fInvMassMuNuPtAll(0),fPtCut1(0), fPtCut2(0),
  62       fPtCut3(0), fAngMomKKinks(0),
  63  f1(0), f2(0),
  64       fListOfHistos(0),fLowMulcut(-1),fUpMulcut(-1), fKinkRadUp(200), fKinkRadLow(130), fLowCluster(20), fLowQt(.12), fCutsMul(0),  fMaxDCAtoVtxCut(0),  fPIDResponse(0)
  65 {
  66   // Constructor
  67
  68   // Define input and output slots here
  69   // Input slot #0 works with a TChain
  70  // DefineInput(0, TChain::Class());
  71 //----------------------Marek multiplicity bins
  72  fCutsMul=new AliESDtrackCuts("Mul","Mul");
  73         fCutsMul->SetMinNClustersTPC(70);
  74         fCutsMul->SetMaxChi2PerClusterTPC(4);
  75         fCutsMul->SetAcceptKinkDaughters(kFALSE);
  76         fCutsMul->SetRequireTPCRefit(kTRUE);
  77         // ITS
  78         fCutsMul->SetRequireITSRefit(kTRUE);
  79         fCutsMul->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
  80                                                 AliESDtrackCuts::kAny);
  81         fCutsMul->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
  82
  83         fCutsMul->SetMaxDCAToVertexZ(2);
  84         fCutsMul->SetDCAToVertex2D(kFALSE);
  85         fCutsMul->SetRequireSigmaToVertex(kFALSE);
  86
  87         fCutsMul->SetEtaRange(-0.8,+0.8);
  88         fCutsMul->SetPtRange(0.15, 1e10);
  89
  90         fMaxDCAtoVtxCut=new AliESDtrackCuts("fMaxDCA", "fMaxDCA");
  91        fMaxDCAtoVtxCut->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
  92     fMaxDCAtoVtxCut->SetMaxChi2TPCConstrainedGlobal(36);
  93 //        esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
  94   //  esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
  95
  96
  97
  98   DefineOutput(1, TList::Class());
  99 }
 100
 101 //________________________________________________________________________
 102 void AliAnalysisKinkESDat::UserCreateOutputObjects()
 103 {
 104   // Create histograms
 105   // Called once
 106
 107    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);
 108    f1->SetParameter(0,0.493677);
 109    f1->SetParameter(1,0.9127037);
 110    f1->SetParameter(2,TMath::Pi());
 111
 112
 113    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);
 114    f2->SetParameter(0,0.13957018);
 115    f2->SetParameter(1,0.2731374);
 116    f2->SetParameter(2,TMath::Pi());
 117    //Open file  1= CAF
 118     //OpenFile(1);
 119    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,
 120                         1.1, 1.2, 1.3, 1.4, 1.5, 1.6,1.7,1.8,1.9,  2.0,
 121                          2.2, 2.4, 2.6, 2.8,  3.0,   3.3, 3.6, 3.9,
 122                          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
 123
 124    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,
 125                         1.1, 1.2, 1.3, 1.4, 1.5, 1.6,1.7,1.8,1.9,  2.0,
 126                          2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7,2.8, 2.9, 3.0,
 127                          3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8,5.0 };  //  Barbara TOF Kch
 128
 129   fHistPtESD = new TH1F("fHistPtESD", "P_{T} distribution",300, 0.0,15.0);
 130   fHistPtESD->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 131   fHistPtESD->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 132   fHistPtESD->SetMarkerStyle(kFullCircle);
 133   fHistPt = new TH1F("fHistPt", "P_{T} distribution",300, 0.0,15.0);
 134   fHistQtAll = new TH1F("fHistQtAll", "Q_{T} distr All Kinks ",100, 0.0,.300);
 135   fHistQt1= new TH1F("fHistQt1", "Q_{T} distribution",100, 0.0,.300);
 136   fHistQt2= new TH1F("fHistQt2", "Q_{T} distribution",100, 0.0,.300);
 137   fHistPtKaon = new TH1F("fHistPtKaon", "P_{T}Kaon distribution",300, 0.0,15.0);
 138   fHistPtKPDG = new TH1F("fHistPtKPDG", "P_{T}Kaon distribution",300, 0.0,15.0);
 139   fHistEta= new TH1F("fHistEta", "Eta distribution", 26,-1.3, 1.3);
 140   fHistEtaK= new TH1F("fHistEtaK", "EtaK distribution", 26,-1.3, 1.3);
 141   fptKMC= new TH1F("fptKMC", "P_{T}Kaon generated",300, 0.0,15.0);
 142   fMultiplMC= new TH1F("fMultiplMC", "charge multiplicity MC",100, 0.0,300.0);
 143   fESDMult= new TH1F("fESDMult", "charge multipliESD",100, 0.0,300.0);
 144   fgenpt= new TH1F("fgenpt", "genpt   K distribution",300, 0.0,15.0);
 145    //frad= new TH1F("frad", "radius  K generated",100, 50., 250.0);
 146    frad= new TH1F("frad", "radius  K generated",100, 0.,1000.0);
 147   fKinkKaon= new TH1F("fKinkKaon", "P_{T}Kaon kinks identi",300, 0.0,15.0);
 148   fKinKRbn= new TH1F("fKinKRbn", "p_{t}Kaon kinks identi[GeV/c],Entries",46,gPt7TOF);
 149   fKinkKaonBg= new TH1F("fKinkKaonBg", "P_{T}Kaon kinks backgr",300, 0.0,15.0);
 150   //fM1kaon= new TH1F("fM1kaon","Invar m(kaon) from kink->mu+netrino decay",180,0.10, 1.0);
 151   fM1kaon= new TH1F("fM1kaon","Invar m(kaon) from kink->mu+netrino decay",600,0.10, 0.7); //  23/8/2013
 152   fPtKink= new TH1F("fPtKink", "P_{T}Kaon Kink  distribution, counts",44, gPt7K0);
 153   fptKink= new TH1F("fptKink", "P_{T}Kaon Kink  bution",300, 0.0,15.0);
 154   fAngMomK= new TH2F("fAngMomK","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.);
 155   fAngMomPi= new TH2F("fAngMomPi","Decay angle vrs Mother Mom,Pi",100,0.0,5.0,80,0.,80.);
 156   fAngMomKC= new TH2F("fAngMomKC","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.);
 157   fMultESDK=new TH1F("fMultESDK", "charge multipliESD kaons",100, 0.0,100.0);
 158   fMultMCK=new TH1F("fMultMCK", "charge multipli MC kaons",100, 0.0,100.0);
 159   fSignPtNcl= new TH2F("fSignPtNcl","SignPt vrs Ncl,K",80,-4.,4.0,70,20.,160.);
 160   fSignPtEta= new TH2F("fSignPtEta","SignPt vrs Eta,K",80,-4.0,4.0,30,-1.5,1.5);
 161   fEtaNcl= new TH2F("fEtaNcl","Eta vrs nclust,K",30,-1.5,1.5, 70,20, 160);
 162   fSignPt= new TH1F("fSignPt","SignPt ,K",80,-4.0,4.0);
 163   fChi2NclTPC= new TH2F("fChi2NclTPC","Chi2vrs nclust,K",100,0.,500., 70,20, 160);
 164   fRatChi2Ncl= new TH1F("fRatChi2Ncl","Ratio chi2/nclusters in TPC,K",50,0.0,5.0);
 165   fRadiusNcl = new TH2F("fRadiusNcl","kink radius vrs Nclust,K",75,100.,250., 80,0, 160);
 166     fTPCSgnlP = new TH2F("fTPCSgnlP","TPC signal de/dx Mom,K",1000,0.0,20.0,150,0.,300.);
 167   fTPCSgnlPa= new TH2F("fTPCSgnlPa","TPC signal de/dx Mom,K",1000,0.0,20.,150, 0.,300.);
 168   fRpr = new TH1D("fRpr", "rad distribution  PID pr",100,-10.0, 10.0);
 169   fZpr = new TH1D("fZpr", "z distribution PID pr  ",80,-20.,20.);
 170   fdcatoVxXY = new TH1D("fdcatoVxXY", "dca  distribution PID  ",20,-1.,1.);
 171   fnSigmToVx = new TH1D("fnSigmToVx", "dca  distribution PID  ",80,0.,8.);
 172   fKinkMothDau= new TH2F("fKinkMothDau","TPC kink Moth Daugh ,K",50,0.0,2.5,50, 0., 2.5);
 173   fZvXv= new TH2F("fZvXv","Xv-Zv main vtx",60,-0.5,0.5,60, -15., 15.0);
 174   fZvYv= new TH2F("fZvYv","Yv-Zv main vtx",60,-0.5,0.5, 60, -15., 15.);
 175   fXvYv= new TH2F("fXvYv","Xv-Yv main vtx", 60,-1.5,1.5, 60, -1.5, 1.5);
 176   fHistPtKaoP = new TH1F("fHistPtKaoP", "P_{T}KaonP  distribution",300, 0.0,15.0);
 177   fHistPtKaoN = new TH1F("fHistPtKaoN", "P_{T}KaonN  distribution",300, 0.0,15.0);
 178   frapiKESD=new TH1F("frapiKESD","rapid Kdistribution", 26,-1.3, 1.3);
 179   flifetime= new TH1F("flifetime", "ct study of K-kinks",100,0.,1000.);
 180   fradLK= new TH1F("fradLK", "Length of   K generated",100,0.,1000.);
 181   fradPtRpDt=new TH3F("fradPtRpDt","rad pt rap dat",28,100.,240., 20, 0., 5., 20, -1., 1. );
 182   //fInvMuNuAll= new TH1F("fInvMuNuAll", " Inv Mass MuNu all kink",180,0.1,1.0);
 183   fInvMuNuAll= new TH1F("fInvMuNuAll", " Inv Mass MuNu all kink",600,0.1,0.7); //  23/8/2013
 184   fQtInvM= new TH2F("fQtInvM", "Q_{T} Versus Inv MuNu ",80, 0., 0.80 , 100 , 0., 0.300);
 185     fDCAkink = new TH1F("fDCAkink ", "DCA kink vetrex ",50, 0.0,1.0);
 186   fPosiKink= new TH2F("fPosiKink", "Y vrx kink Vrex ",100, -300.0,300.0,100, -300, 300.);
 187   fPosiKinkK= new TH2F("fPosiKinkK", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
 188   fPosiKinKXZ= new TH2F("fPosiKinKXZ", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
 189   fPosiKinKYZ= new TH2F("fPosiKinKYZ", "Y vrx kink VrexK ",100, -300.0,300.0,100, -300, 300.);
 190   fPosiKinKBg= new TH2F("fPosiKinKBg", "z vrx kink rad    ",100, -300.0,300.0,100,100., 300.);
 191   fQtMothP = new TH2F("fQtMothP", " Qt vrs Mother P", 100, 0., 5.0,100, 0.,0.300);
 192     fTPCSgnlPtpc = new TH2F("fTPCSgnlPtpc","TPC signal de/dx Mom TPC,K  ",300,0.0,15.0,100, 0., 250.    );
 193     fTPCMomNSgnl = new TH2F("fTPCMomNsgnl","TPC signal de/dx Mom TPC,K  ",300,0.0,15.0,20 , -10., 10.);
 194     fMothKinkMomSgnl  = new TH2F("fMothKinkMomSgnl","TPC signal de/dx Mom TPC,Kink  ",100,0.0,250.0,100, 0., 250.    );
 195     fNSigmTPC    = new TH1F("fNSigmTPC","TPC Nsigma  de/dx  TPC,K  ", 30 , -7.5, 7.5);
 196     fTPCSgnlKinkDau = new TH2F("fTPCSgnlKinkDau","TPC signal de/dx Mom,K",500,0.0,10.0,100,0.,250.);
 197   fPtKinkPos= new TH1F("fPtKinkPos", "Pos P_{T}Kaon Kink  distribution, counts",44, gPt7K0);
 198   fPtKinkNeg= new TH1F("fPtKinkNeg", "Neg P_{T}Kaon Kink  distribution, counts",44, gPt7K0);
 199   fRadNclCln = new TH2F("fRadNclCln","kink radius vrs Nclust,K Clean ",75,100.,250., 80,0, 160);
 200   fRatioCrossedRows = new TH1F("fRatioCrossedRows","Ratio crossed rows  in TPC",20,0.0,1.0);
 201   fRatioCrossedRowsKink = new TH1F("fRatioCrossedRowsKink","Ratio crossed rows  in TPC for kinks",20,0.0,1.0);
 202   fRadiusPt =new TH2F("fRadiusPt","radius vs pt  ",80, 90.,250.,100, 0.,10.              );
 203   fRadiusPtcln =new TH2F("fRadiusPtcln","radius vs pt clean ",80, 90.,250.,100, 0.,10.              );
 204   //fInvMassMuNuPt =new TH2F("fInvMassMuNuPt","Invariant mass-munu  vs pt  ",180, 0.10, 1.00, 100, 0.0, 10.0  );
 205   fInvMassMuNuPt =new TH2F("fInvMassMuNuPt","Invariant mass-munu  vs pt  ",600, 0.10, 0.7, 100, 0.0, 10.0  );// 23/8/2013
 206   fInvMassMuNuPtAll =new TH2F("fInvMassMuNuPtAll","Invariant mass-munu  vs pt  ",600, 0.10, 0.7, 100, 0.0, 10.0  );// 23/8/2013
 207   fPtCut1 = new TH1F("fPtCut1", "P_{T}Kaon distribution",300, 0.0,15.0);
 208   fPtCut2 = new TH1F("fPtCut2", "P_{T}Kaon distribution",300, 0.0,15.0);
 209   fPtCut3 = new TH1F("fPtCut3", "P_{T}Kaon distribution",300, 0.0,15.0);
 210   fAngMomKKinks = new TH2F("fAngMomKKinks","Decay angle vrs Mother Mom,Kinks",300,0.0,15.0,100,0.,100.);
 211
 212    fListOfHistos=new TList();
 213
 214    fListOfHistos->Add(fHistPtESD);
 215    fListOfHistos->Add(fHistPt);
 216    fListOfHistos->Add(fHistQtAll);
 217    fListOfHistos->Add(fHistQt1);
 218    fListOfHistos->Add(fHistQt2);
 219    fListOfHistos->Add(fHistPtKaon);
 220    fListOfHistos->Add(fHistPtKPDG);
 221    fListOfHistos->Add(fHistEta);
 222    fListOfHistos->Add(fHistEtaK);
 223    fListOfHistos->Add(fptKMC);
 224    fListOfHistos->Add(fMultiplMC);
 225    fListOfHistos->Add(fESDMult);
 226    fListOfHistos->Add(fgenpt);
 227    fListOfHistos->Add(frad);
 228    fListOfHistos->Add(fKinkKaon);
 229    fListOfHistos->Add(fKinKRbn);
 230    fListOfHistos->Add(fKinkKaonBg);
 231    fListOfHistos->Add(fM1kaon);
 232    fListOfHistos->Add(fPtKink);
 233    fListOfHistos->Add(fptKink);
 234    fListOfHistos->Add(fAngMomK);
 235    fListOfHistos->Add(fAngMomPi);
 236    fListOfHistos->Add(fAngMomKC);
 237    fListOfHistos->Add(fMultESDK);
 238    fListOfHistos->Add(fMultMCK);
 239    fListOfHistos->Add(fSignPtNcl);
 240    fListOfHistos->Add(fSignPtEta);
 241    fListOfHistos->Add(fEtaNcl);
 242    fListOfHistos->Add(fSignPt);
 243    fListOfHistos->Add(fChi2NclTPC);
 244    fListOfHistos->Add(fRatChi2Ncl);
 245    fListOfHistos->Add(fRadiusNcl);
 246    fListOfHistos->Add(fTPCSgnlP);
 247    fListOfHistos->Add(fTPCSgnlPa);
 248    fListOfHistos->Add(fRpr);
 249    fListOfHistos->Add(fZpr);
 250    fListOfHistos->Add(fdcatoVxXY);
 251    fListOfHistos->Add(fnSigmToVx);
 252    fListOfHistos->Add(fKinkMothDau);
 253    fListOfHistos->Add(fZvXv);
 254    fListOfHistos->Add(fZvYv);
 255    fListOfHistos->Add(fXvYv);
 256    fListOfHistos->Add(fHistPtKaoP);
 257    fListOfHistos->Add(fHistPtKaoN);
 258    fListOfHistos->Add(frapiKESD);
 259    fListOfHistos->Add(flifetime);
 260    fListOfHistos->Add(fradLK);
 261    fListOfHistos->Add(fradPtRpDt);
 262    fListOfHistos->Add(fInvMuNuAll);
 263    fListOfHistos->Add(fQtInvM);
 264    fListOfHistos->Add(fDCAkink);
 265    fListOfHistos->Add(fPosiKink);
 266    fListOfHistos->Add(fPosiKinkK);
 267    fListOfHistos->Add(fPosiKinKXZ);
 268    fListOfHistos->Add(fPosiKinKYZ);
 269    fListOfHistos->Add(fPosiKinKBg);
 270    fListOfHistos->Add(fQtMothP);
 271    fListOfHistos->Add(fTPCSgnlPtpc);
 272    fListOfHistos->Add(fTPCMomNSgnl);
 273    fListOfHistos->Add(fMothKinkMomSgnl);
 274    fListOfHistos->Add(fNSigmTPC);
 275    fListOfHistos->Add(fTPCSgnlKinkDau);
 276    fListOfHistos->Add(fPtKinkPos);
 277    fListOfHistos->Add(fPtKinkNeg);
 278    fListOfHistos->Add(fRadNclCln);
 279    fListOfHistos->Add(fRatioCrossedRows);
 280    fListOfHistos->Add(fRatioCrossedRowsKink);
 281    fListOfHistos->Add(fRadiusPt);
 282    fListOfHistos->Add(fRadiusPtcln);
 283    fListOfHistos->Add(fInvMassMuNuPt);
 284    fListOfHistos->Add(fInvMassMuNuPtAll);
 285    fListOfHistos->Add(fPtCut1);
 286    fListOfHistos->Add(fPtCut2);
 287    fListOfHistos->Add(fPtCut3);
 288   fListOfHistos->Add(fAngMomKKinks);
 289
 290 }
 291 //________________________________________________________________________
 292 void AliAnalysisKinkESDat::UserExec(Option_t *)
 293 {
 294   // Main loop
 295   // Called for each event
 296
 297   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 298   // This handler can return the current MC event
 299
 300    AliVEvent *event = InputEvent();
 301   if (!event) {
 302      Printf("ERROR: Could not retrieve event");
 303      return;
 304   }
 305
 306   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 307   if (!esd) {
 308      Printf("ERROR: Could not retrieve esd");
 309      return;
 310   }
 311 // Number ESD tracks
 312    Int_t nESDTracks =  esd->GetNumberOfTracks();
 313       fMultiplMC->Fill(nESDTracks);
 314 //
 315 //==================check of Physics selection?
 316        Bool_t isSelected =
 317 ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()&AliVEvent::kMB;
 318
 319          if ( isSelected ==kFALSE) return;   //  24/6/11 apo MF
 320 //
 321       fMultMCK->Fill(nESDTracks);
 322 //
 323 //===============Marek's  multiplicity
 324    Float_t refmultiplicity=fCutsMul->CountAcceptedTracks(esd);
 325         if(fLowMulcut>-1)
 326         {
 327                 if(refmultiplicity<fLowMulcut)
 328                         return;
 329         }
 330         if(fUpMulcut>-1)
 331         {
 332                 if(refmultiplicity>fUpMulcut)
 333                         return;
 334         }
 335
 336
 337
 338        fMultESDK->Fill(refmultiplicity);
 339
 340 //
 341
 342
 343
 344   const AliESDVertex *vertex=GetEventVertex(esd);    // 22/8
 345   if(!vertex) return;
 346 //
 347   Double_t vpos[3];
 348   vertex->GetXYZ(vpos);
 349     fZpr->Fill(vpos[2]);
 350       if (TMath::Abs( vpos[2] ) > 10. ) return;
 351
 352
 353
 354   Double_t vtrack[3], ptrack[3];
 355
 356
 357  Int_t nESDTracK = 0;
 358 // Int_t nESDTrKink = 0;
 359
 360    Int_t nGoodTracks =  esd->GetNumberOfTracks();
 361     fESDMult->Fill(nGoodTracks);
 362
 363        Double_t nsigmall = 100.0;
 364        Double_t nsigma = 100.0;
 365        Double_t nsigmaPion =-100.0;
 366   //     Double_t nsigmaDau  =-100.0;
 367        Double_t dEdxKinkDau =0.0;
 368        Double_t KinkDauCl   =0.0;
 369 // apo Eftihi
 370                   if(!fPIDResponse) {
 371     AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 372     AliInputEventHandler* inputHandler =
 373 (AliInputEventHandler*)(man->GetInputEventHandler());
 374     fPIDResponse = inputHandler->GetPIDResponse();
 375   }
 376
 377 // loop on kink daughters
 378    for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 379
 380     AliESDtrack* trackD = esd->GetTrack(iTrack);
 381     if (!trackD) {
 382       Printf("ERROR: Could not receive track %d", iTrack);
 383       continue;
 384     }
 385 //
 386                 Int_t indexKinkDau=trackD->GetKinkIndex(0);
 387 // daughter kink
 388 //        AliESDkink *kink=esd->GetKink(TMath::Abs(indexKinkDau)-1);
 389       if ( indexKinkDau > 0 )    {
 390     Int_t labelD = trackD->GetLabel();
 391     labelD = TMath::Abs(labelD);
 392    nsigmaPion     = (fPIDResponse->NumberOfSigmasTPC(trackD  , AliPID::kPion));// 26/10 eftihis
 393      dEdxKinkDau =  (trackD->GetTPCsignal()  )  ;  //  daughter kink  dEdx
 394    KinkDauCl=(trackD->GetTPCclusters(0)  )  ;
 395  }
 396 //if((indexKinkDau >0)&& (nsigmaPion>1.2)) fTPCSgnlKinkDau->Fill(trackD->P(), (trackD->GetTPCsignal()  ) ) ;  //  daughter kink
 397 // Ayto mexri 26/11/2012     if(indexKinkDau >0) fTPCSgnlKinkDau->Fill(trackD->P(), (trackD->GetTPCsignal()  ) ) ;  //  daughter kink
 398    }
 399
 400 // track loop
 401 //
 402    for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
 403
 404     AliESDtrack* track = esd->GetTrack(iTracks);
 405     if (!track) {
 406       Printf("ERROR: Could not receive track %d", iTracks);
 407       continue;
 408     }
 409
 410     fHistPt->Fill(track->Pt());
 411
 412
 413      //    sigmas
 414      nsigmall  = (fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon));
 415   //  nsigmaPion= (fESDpid->NumberOfSigmasTPC(track,AliPID::kPion));
 416       nsigma  = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon));
 417
 418 //=======================new
 419      Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 420   Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 421   if (track->GetTPCNclsF()>0) {
 422     ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC / track->GetTPCNclsF();
 423     fRatioCrossedRows->Fill(ratioCrossedRowsOverFindableClustersTPC);
 424    }
 425 //_______
 426
 427                 Int_t indexKinkPos=track->GetKinkIndex(0);   // kink index
 428
 429       Int_t tpcNCl = track->GetTPCclusters(0);
 430       Double_t tpcSign = track->GetSign();
 431
 432     Int_t label = track->GetLabel();
 433     label = TMath::Abs(label);
 434
 435
 436     UInt_t status=track->GetStatus();
 437
 438     if((status&AliESDtrack::kITSrefit)==0) continue;
 439     if((status&AliESDtrack::kTPCrefit)==0) continue;
 440       if((track->GetTPCchi2()/track->GetTPCclusters(0))>4.0) continue;
 441
 442       Double_t extCovPos[15];
 443       track->GetExternalCovariance(extCovPos);
 444
 445
 446     track->GetXYZ(vtrack);
 447  fXvYv->Fill(vtrack[0],vtrack[1]);
 448  fZvYv->Fill(vtrack[0],vtrack[2]);
 449  fZvXv->Fill(vtrack[1],vtrack[2]);
 450
 451 // track momentum, rapidity calculation
 452      track->GetPxPyPz(ptrack);
 453
 454     TVector3 trackMom(ptrack[0],ptrack[1],ptrack[2]);
 455
 456 // K-rapidity calcualtion
 457           Double_t   etracK= TMath::Sqrt(trackMom.Mag()*trackMom.Mag() + 0.493677 *0.493677  );
 458          Double_t rapiditK = 0.5 * (TMath::Log(  (etracK + ptrack[2]  ) / ( etracK - ptrack[2])  ))  ;
 459
 460     Double_t trackEta=trackMom.Eta();
 461     Double_t trackPt = track->Pt();
 462
 463
 464
 465     Float_t bpos[2];
 466     Float_t bCovpos[3];
 467     track->GetImpactParameters(bpos,bCovpos);
 468
 469     if (bCovpos[0]<=0 || bCovpos[2]<=0) {
 470      Printf("Estimated b resolution lower or equal zero!");
 471      bCovpos[0]=0; bCovpos[2]=0;
 472     }
 473
 474     Float_t dcaToVertexXYpos = bpos[0];
 475     Float_t dcaToVertexZpos = bpos[1];
 476
 477     fRpr->Fill(dcaToVertexZpos);
 478
 479 //  14/2/13 /================/   if((TMath::Abs(dcaToVertexXYpos)>0.3)||(TMath::Abs(dcaToVertexZpos)>2.5))  nESDTrKink++;  //  count of second  23Jul11
 480  //if((TMath::Abs(dcaToVertexXYpos)>0.3)||(TMath::Abs(dcaToVertexZpos)>2.5))
 481 //          continue;   //    allagi  23Jul11
 482
 483                     if (!fMaxDCAtoVtxCut->AcceptTrack(track)) continue;
 484
 485     fdcatoVxXY->Fill(dcaToVertexXYpos);
 486 //
 487
 488 //  track Mult. after selection
 489     nESDTracK++;
 490   //
 491 //=========================================
 492     fHistPtESD->Fill(track->Pt());
 493
 494    // Add Kink analysis           =============================
 495
 496 // daughter kink
 497 //if(indexKinkPos >0)fTPCSgnlKinkDau->Fill(track->P(), (track->GetTPCsignal()  ) ) ;  //  daughter kink
 498
 499 //  loop on kinks
 500                 if(indexKinkPos<0){     ////mother kink
 501
 502              fptKMC   ->Fill(  track->Pt()    );  // Pt from tracks , all kinks
 503
 504     fRatioCrossedRowsKink->Fill(ratioCrossedRowsOverFindableClustersTPC);
 505         // select kink class
 506
 507           AliESDkink *kink=esd->GetKink(TMath::Abs(indexKinkPos)-1);
 508 //
 509          // DCA kink
 510           Double_t  Dist2 = kink->GetDistance();
 511          //   fDCAkink->Fill( Dist2   );
 512             //   if (Dist2 > 0.2) continue; //  systematics 11/8/11
 513 //
 514
 515 // TPC mother momentum
 516
 517            const TVector3 vposKink(kink->GetPosition());
 518  fPosiKink ->Fill( vposKink[0], vposKink[1]  );
 519    Double_t  dzKink=vpos[2]-vposKink[2];
 520 //
 521 //   lifitime
 522             Double_t tanLamda = track->GetTgl();  // 25/6/2010
 523
 524    Double_t lifeKink= (TMath::Abs( dzKink ))*( TMath::Sqrt(1.+ tanLamda*tanLamda) ) / (TMath::Abs( tanLamda)) ;
 525 //
 526            const TVector3 motherMfromKink(kink->GetMotherP());
 527            const TVector3 daughterMKink(kink->GetDaughterP());
 528
 529            Float_t qT=kink->GetQt();
 530             Float_t motherPt=motherMfromKink.Pt();
 531 // Kink  mother momentum
 532 //     Double_t trMomTPCKink=motherMfromKink.Mag();
 533 // TPC mother momentun
 534      Double_t trMomTPC=track->GetTPCmomentum();
 535   //     fTPCSgnlKinkDau->Fill( daughterMKink.Mag(), dEdxKinkDau  ) ;  //  daughter kink
 536   //
 537            fHistQtAll->Fill(qT) ;  //  Qt   distr
 538
 539            fptKink->Fill(motherMfromKink.Pt()); /// pt from kink
 540
 541            Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
 542
 543 //   rapiditya nd pt selection
 544           if(  (TMath::Abs(rapiditK )) > 0.7 ) continue;
 545         if ( (track->Pt())<.200)continue;  // new Feb 2012
 546
 547                 fQtMothP->Fill( track->P(), qT);
 548
 549         if ( qT> fLowQt )  fHistQt1  ->Fill(qT) ;  //  Qt   distr
 550
 551
 552
 553             fHistEta->Fill(trackEta) ;  //   Eta distr
 554           fHistQt2->Fill(qT);  //
 555             fKinkKaonBg->Fill(motherPt);
 556
 557 //          maximum decay angle at a given mother momentum
 558            //Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
 559            Double_t maxDecAngKmu=f1->Eval(track->P()          ,0.,0.,0.);
 560            Double_t maxDecAngpimu=f2->Eval(    track->P(),       0.,0.,0.);
 561
 562 //  fake kinks are removed
 563          if( (kinkAngle<2.)  ) continue;
 564
 565
 566       //  BG  ?????==============
 567               if ( TMath::Abs(vposKink[2]) >  225. ) continue ;
 568               if ( TMath::Abs(vposKink[2]) <  0.5 ) continue ;
 569 //
 570             fPtCut1   ->Fill(trackPt );
 571             fAngMomPi->Fill( track->P(),           kinkAngle);
 572 //
 573 // invariant mass of mother track decaying to mu
 574          Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
 575          Float_t p1XM= motherMfromKink.Px();
 576          Float_t p1YM= motherMfromKink.Py();
 577          Float_t p1ZM= motherMfromKink.Pz();
 578          Float_t p2XM= daughterMKink.Px();
 579          Float_t p2YM= daughterMKink.Py();
 580          Float_t p2ZM= daughterMKink.Pz();
 581          Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
 582          Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
 583            fQtInvM -> Fill ( invariantMassKmu,  qT);
 584            fInvMuNuAll->Fill(invariantMassKmu);
 585            fInvMassMuNuPtAll ->Fill(invariantMassKmu,  trackPt);
 586 //
 587          fRadiusPt->Fill( kink->GetR(), trackPt); //
 588  //  radius and Minv selection
 589        //if( ( kink->GetR()> 120 ) && ( kink->GetR() < 210 )  )  {
 590        if( ( kink->GetR()> fKinkRadLow ) && ( kink->GetR() <fKinkRadUp   )  )  {
 591     //  for systematics   if( ( kink->GetR()> 130 ) && ( kink->GetR() < 200 )  )  {
 592       if (qT>fLowQt )  fAngMomKC->Fill(track->P(), kinkAngle);
 593           if ( qT> fLowQt ) fM1kaon->Fill(invariantMassKmu);
 594              if ( qT > fLowQt)
 595          fRadiusNcl->Fill( (kink->GetR()) ,(track->GetTPCclusters(0)  ) ) ;
 596   }
 597 //  tails cleaning
 598                if(  ( tpcNCl<fLowCluster) ) continue;  // test 27 feb 2012 ,, OK
 599           //  edw iatn !!!    if(  ( tpcNCl<50 ) ) continue;  // test 15 March  13,, OK
 600 // cleaning BG in tails
 601       Int_t tpcNClHigh = -51.67+ (11./12.)  *( kink->GetR() ) ;
 602                if ( tpcNCl > tpcNClHigh) continue;
 603
 604       Int_t tpcNClMin  = -85.5 + (65./95.)  *( kink->GetR() ) ;
 605                if ( tpcNCl < tpcNClMin ) continue;
 606
 607    //  back, 20/1/2013   if (ratioCrossedRowsOverFindableClustersTPC< 0.5) continue;// test 14/1/2013
 608 //
 609                fHistPtKPDG->Fill(track->Pt());  // ALL  K-candidates until now
 610     //  if((kinkAngle>maxDecAngpimu)&&(qT>0.12)&&(qT<0.30)&&((kink->GetR()>=120.)&&(kink->GetR()<=210.))&&(TMath::Abs(rapiditK)<0.7)&&(invariantMassKmu<0.6)){
 611      //if((kinkAngle>maxDecAngpimu)&&(qT>0.12)&&(qT<0.30)&&((kink->GetR()>=120.)&&(kink->GetR()<=210.))&&(TMath::Abs(rapiditK)<0.7)&&(invariantMassKmu<0.8)){
 612      //if((kinkAngle>maxDecAngpimu)&&(qT>0.04)&&(qT<0.30)&&((kink->GetR()>= fKinkRadLow )&&(kink->GetR()<= fKinkRadUp ))&&(TMath::Abs(rapiditK)<0.7)&&(invariantMassKmu<0.8)){
 613      if((kinkAngle>maxDecAngpimu)&&(qT>fLowQt)&&(qT<0.30)&&((kink->GetR()>= fKinkRadLow )&&(kink->GetR()<= fKinkRadUp ))&&(TMath::Abs(rapiditK)<0.7)&&(invariantMassKmu<0.8)){
 614   // systematics   if((kinkAngle>maxDecAngpimu)&&(qT>0.12)&&(qT<0.30)&&((kink->GetR()>=130.)&&(kink->GetR()<=200.))&&(TMath::Abs(rapiditK)<0.7)&&(invariantMassKmu<0.8)){
 615 //
 616         fAngMomKKinks->Fill(track->P(), kinkAngle);
 617             fPtCut2   ->Fill(trackPt );
 618 //  maximum angles selection with some error cut
 619         if( (kinkAngle<maxDecAngpimu*1.2) ) continue;
 620                  if ( (kinkAngle>maxDecAngKmu*.98) && ( track->P() >1.2 )) continue;  ///5/5/2010
 621
 622             fPtCut3   ->Fill(trackPt );
 623 //  here the kaons selected by the decay features
 624            fTPCSgnlPa->Fill( track->GetInnerParam()->GetP() ,(track->GetTPCsignal()  ) ) ;
 625 //
 626             //  NO dEdx cut test 9/2/13               if ( nsigma               > 3.5) continue;
 627                      if ( nsigma               > 3.5) continue;
 628 //
 629 //  next plots for the identified kaons by the kink analysis
 630
 631                                      fHistPtKaon->Fill(track->Pt());   //
 632               if(tpcSign >0.)        fHistPtKaoP->Fill( track->Pt()         ) ;   //
 633                if ( tpcSign <0.)    fHistPtKaoN->Fill(  track->Pt()        ) ;   //
 634           fTPCSgnlP->Fill(track->GetInnerParam()->GetP(), (track->GetTPCsignal()  ) ) ;
 635          fRadNclCln->Fill( (kink->GetR()) ,(track->GetTPCclusters(0)  ) ) ;
 636          fRadiusPtcln->Fill( kink->GetR(), trackPt); //
 637            fInvMassMuNuPt ->Fill(invariantMassKmu,  trackPt);
 638
 639       //   fTPCSgnlPtpc->Fill(trMomTPC  , (track->GetTPCsignal()  ) ) ;
 640          //fMothKinkMomSgnl ->Fill(trMomTPCKink  , (track->GetTPCsignal()  ) ) ;
 641          fMothKinkMomSgnl ->Fill(  dEdxKinkDau  , (track->GetTPCsignal()  ) ) ;
 642 //
 643        fTPCSgnlKinkDau->Fill( daughterMKink.Mag(), dEdxKinkDau  ) ;  //  daughter kink
 644 //
 645          fTPCMomNSgnl->Fill(trMomTPC ,nsigmall );
 646          fNSigmTPC   ->Fill(nsigmall );
 647 //
 648                 frad->Fill(kink->GetR());  // kink
 649                fradLK->Fill(lifeKink    );  // kink
 650              fHistEtaK->Fill(trackEta);
 651             frapiKESD ->Fill(rapiditK);  //  rapidityof kaons
 652           fPosiKinKBg->Fill( vposKink[2], kink->GetR() );
 653
 654                      Float_t signPt= tpcSign*trackPt;
 655                   fSignPtNcl->Fill( signPt  ,   tpcNCl   );   ///  28/4/2010
 656                   fSignPtEta->Fill( signPt  , rapiditK  );
 657                   fEtaNcl->Fill( rapiditK, tpcNCl    );
 658                   fSignPt->Fill( signPt );
 659                   fChi2NclTPC->Fill( (track->GetTPCchi2() ) ,  tpcNCl );
 660          fRatChi2Ncl-> Fill (  (track->GetTPCchi2()/track->GetTPCclusters(0)  )) ;
 661  //   fdcatoVxXY->Fill(dcaToVertexXYpos);
 662            //if( dEdxKinkDau> 1.5* (track->GetTPCsignal()   )  )      fKinkMothDau->Fill(track->P(),daughterMKink.Mag());
 663        //    if((dEdxKinkDau>  80. ) && (dEdxKinkDau > 4.*nsigmaPion)   )      fKinkMothDau->Fill(track->P(),daughterMKink.Mag());
 664          if (nsigmaPion>  3.             ) fTPCSgnlPtpc->Fill( daughterMKink.Mag(),  dEdxKinkDau    ) ;
 665          //if (TMath::Abs(dEdxKinkDau -(track->GetTPCsignal() )> 10. )) fTPCSgnlPtpc->Fill( daughterMKink.Mag(),  dEdxKinkDau    ) ;
 666                         flifetime->Fill(( lifeKink*.493667   )  /track->P()   ) ;
 667              fKinkKaon->Fill(track->Pt());
 668             fDCAkink->Fill( Dist2   );
 669
 670                fPtKink->Fill(track->Pt()); ///  K0 bins
 671               if(tpcSign >0.)        fPtKinkPos ->Fill( track->Pt()         ) ;   //K-plus bins K0
 672               if(tpcSign <0.)        fPtKinkNeg ->Fill( track->Pt()         ) ;   //K-minus bins K0
 673              fKinKRbn->Fill(track->Pt());       // TOF
 674
 675               fradPtRpDt->Fill( kink->GetR(), 1./track->Pt(), rapiditK);
 676                fAngMomK->Fill(    track->P(),        kinkAngle);
 677        fPosiKinkK->Fill( vposKink[0], vposKink[1]  );
 678           fPosiKinKXZ->Fill( vposKink[2], vposKink[0]  );
 679           fPosiKinKYZ->Fill( vposKink[2], vposKink[1]  );
 680
 681         }  //  kink selection
 682
 683
 684         }  //End Kink Information
 685
 686
 687   } //track loop
 688 //  } //daughter loop
 689
 690  //     fMultiplMC->Fill(nESDTracK );
 691
 692   PostData(1, fListOfHistos);
 693
 694 }
 695
 696 //________________________________________________________________________
 697 void AliAnalysisKinkESDat::Terminate(Option_t *)
 698 {
 699   // Draw result to the screen
 700   // Called once at the end of the query
 701
 702 }
 703
 704 //____________________________________________________________________//
 705
 706
 707 const AliESDVertex* AliAnalysisKinkESDat::GetEventVertex(const AliESDEvent* esd) const
 708   // Get the vertex from the ESD and returns it if the vertex is valid
 709
 710 {
 711   // Get the vertex
 712
 713    const AliESDVertex* vertex = esd->GetPrimaryVertexTracks();
 714
 715   if((vertex->GetStatus()==kTRUE)) return vertex;
 716   else
 717   {
 718      vertex = esd->GetPrimaryVertexSPD();
 719       if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>0)) return vertex;
 720      else
 721      return 0;
 722   }
 723 }