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