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