PWGLF/SPECTRA/Kinks/AliAnalysisKinkESDatION.cxx

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