PWGPP/pid/AliAnalysisTaskK0sBayes.cxx

   1 #include "AliAnalysisTaskK0sBayes.h"
   2
   3 // ROOT includes
   4 #include <TMath.h>
   5
   6 // AliRoot includes
   7 #include "AliInputEventHandler.h"
   8 #include "AliAODEvent.h"
   9 #include "AliAODVertex.h"
  10 #include "AliAODTrack.h"
  11 #include "AliESDtrack.h"
  12 #include "AliCentrality.h"
  13 #include "AliVHeader.h"
  14 #include "AliAODVZERO.h"
  15 #include "TFile.h"
  16 #include "AliOADBContainer.h"
  17 #include "TH2F.h"
  18 #include "TF1.h"
  19 #include "AliGenHijingEventHeader.h"
  20 #include "AliMCEvent.h"
  21 #include "AliAODMCHeader.h"
  22 #include "AliAODMCParticle.h"
  23 #include "TChain.h"
  24 #include "AliESDtrackCuts.h"
  25 #include "AliESDVertex.h"
  26 #include "AliEventplane.h"
  27 #include "AliAnalysisManager.h"
  28 #include "TRandom.h"
  29
  30
  31 Float_t AliAnalysisTaskK0sBayes::fPtKsMin[AliAnalysisTaskK0sBayes::nPtBin] = {0.,0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4.,5.,6.};// ptmin bin
  32 Float_t AliAnalysisTaskK0sBayes::fPtKsMax[AliAnalysisTaskK0sBayes::nPtBin] = {0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4.,5.,6.,10.};// ptmax bin
  33
  34 // STL includes
  35 //#include <iostream>
  36 //using namespace std;
  37
  38 ClassImp(AliAnalysisTaskK0sBayes)
  39 //_____________________________________________________________________________
  40 AliAnalysisTaskK0sBayes::AliAnalysisTaskK0sBayes():
  41   AliAnalysisTaskSE(),
  42   fVtxCut(10.0),  // cut on |vertex| < fVtxCut
  43   fEtaCut(0.8),   // cut on |eta| < fEtaCut
  44   fMinPt(0.15),   // cut on pt > fMinPt
  45   fIsMC(kFALSE),
  46   fQAsw(kFALSE),
  47   fNcluster(70),
  48   fFilterBit(1),
  49   fList(new TList()),
  50   fList2(new TList()),
  51   fList3(new TList()),
  52   fCentrality(-1),
  53   fPsi(0),
  54   fPtKs(0.),
  55   fPhiKs(0.),
  56   fEtaKs(0.),
  57   fMassV0(0.),
  58   fPtKp(0.),
  59   fPhiKp(0.),
  60   fEtaKp(0.),
  61   fPtKn(0.),
  62   fPhiKn(0.),
  63   fEtaKn(0.),
  64   fPidKp(0),
  65   fPidKn(0),
  66   fTOFTPCsignal(0),
  67   fCombsignal(0),
  68   fCutsDaughter(NULL),
  69   fPIDCombined(NULL),
  70   fContPid(NULL),
  71   fContPid2(NULL),
  72   fNK0s(0),
  73   fNpiPos(0),
  74   fNpiNeg(0),
  75   fHmismTOF(0),
  76   fHchannelTOFdistr(0),
  77   fTypeCol(2)
  78 {
  79   // Default constructor (should not be used)
  80   fList->SetName("contKsBayes1");
  81   fList2->SetName("contKsBayes2");
  82   fList3->SetName("contKsBayes3");
  83
  84   fList->SetOwner(kTRUE);
  85   fList2->SetOwner(kTRUE);
  86   fList3->SetOwner(kTRUE);
  87
  88   TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
  89   fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
  90
  91   TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
  92   fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
  93
  94   for(Int_t i=0;i < nCentrBin;i++){
  95     fElTOF[i] = NULL;
  96     fElTPC[i] = NULL;
  97     fPiTOF[i] = NULL;
  98     fPiTPC[i] = NULL;
  99     fKaTOF[i] = NULL;
 100     fKaTPC[i] = NULL;
 101     fPrTOF[i] = NULL;
 102     fPrTPC[i] = NULL;
 103   }
 104   for(Int_t i=0;i < 4;i++){
 105     hMatching[i] = NULL;
 106     hTracking[i] = NULL;
 107   }
 108   for(Int_t i=0;i < 1000;i++){
 109     fPhiK0s[i] = 0.0;
 110     fPtK0s[i] = 0.0;
 111     fIPiPos[i] = 0;
 112     fIPiNeg[i] = 0;
 113     fIpiP[i] = 0;
 114     fIpiN[i] = 0;
 115     fMassKs[i] = 0.0;
 116   }
 117 }
 118
 119 //______________________________________________________________________________
 120 AliAnalysisTaskK0sBayes::AliAnalysisTaskK0sBayes(const char *name):
 121   AliAnalysisTaskSE(name),
 122   fVtxCut(10.0),  // cut on |vertex| < fVtxCut
 123   fEtaCut(0.8),   // cut on |eta| < fEtaCut
 124   fMinPt(0.15),   // cut on pt > fMinPt
 125   fIsMC(kFALSE),
 126   fQAsw(kFALSE),
 127   fNcluster(70),
 128   fFilterBit(1),
 129   fList(new TList()),
 130   fList2(new TList()),
 131   fList3(new TList()),
 132   fCentrality(-1),
 133   fPsi(0),
 134   fPtKs(0.),
 135   fPhiKs(0.),
 136   fEtaKs(0.),
 137   fMassV0(0.),
 138   fPtKp(0.),
 139   fPhiKp(0.),
 140   fEtaKp(0.),
 141   fPtKn(0.),
 142   fPhiKn(0.),
 143   fEtaKn(0.),
 144   fPidKp(0),
 145   fPidKn(0),
 146   fTOFTPCsignal(0),
 147   fCombsignal(0),
 148   fCutsDaughter(NULL),
 149   fPIDCombined(NULL),
 150   fContPid(NULL),
 151   fContPid2(NULL),
 152   fNK0s(0),
 153   fNpiPos(0),
 154   fNpiNeg(0),
 155   fHmismTOF(0),
 156   fHchannelTOFdistr(0),
 157   fTypeCol(2)
 158 {
 159
 160   DefineOutput(1, TList::Class());
 161   DefineOutput(2, TList::Class());
 162   DefineOutput(3, TList::Class());
 163   DefineOutput(4, TList::Class());
 164
 165   // Output slot #1 writes into a TTree
 166   fList->SetName("contKsBayes1");
 167   fList2->SetName("contKsBayes2");
 168   fList3->SetName("contKsBayes3");
 169
 170   fList->SetOwner(kTRUE);
 171   fList2->SetOwner(kTRUE);
 172   fList3->SetOwner(kTRUE);
 173
 174   TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
 175   fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
 176
 177   TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
 178   fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
 179
 180   for(Int_t i=0;i < nCentrBin;i++){
 181     fElTOF[i] = NULL;
 182     fElTPC[i] = NULL;
 183     fPiTOF[i] = NULL;
 184     fPiTPC[i] = NULL;
 185     fKaTOF[i] = NULL;
 186     fKaTPC[i] = NULL;
 187     fPrTOF[i] = NULL;
 188     fPrTPC[i] = NULL;
 189   }
 190   for(Int_t i=0;i < 4;i++){
 191     hMatching[i] = NULL;
 192     hTracking[i] = NULL;
 193   }
 194   for(Int_t i=0;i < 1000;i++){
 195     fPhiK0s[i] = 0.0;
 196     fPtK0s[i] = 0.0;
 197     fIPiPos[i] = 0;
 198     fIPiNeg[i] = 0;
 199     fIpiP[i] = 0;
 200     fIpiN[i] = 0;
 201     fMassKs[i] = 0.0;
 202   }
 203 }
 204 //_____________________________________________________________________________
 205 AliAnalysisTaskK0sBayes::~AliAnalysisTaskK0sBayes()
 206 {
 207
 208 }
 209
 210 //______________________________________________________________________________
 211 void AliAnalysisTaskK0sBayes::UserCreateOutputObjects()
 212 {
 213
 214   fPIDCombined=new AliPIDCombined;
 215   fPIDCombined->SetDefaultTPCPriors();
 216   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
 217
 218   Float_t invMmin = 0.497-0.005*8;
 219   Float_t invMmax = 0.497+0.005*8;
 220
 221   const Int_t nBinPid = 14;
 222
 223   Int_t binPid[nBinPid] = {1/*ptKs*/,8/*EtaPiP*/,20/*pt+*/,1/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,1/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,1/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
 224
 225   Int_t binPid2[nBinPid] = {1/*ptKs*/,8/*EtaPiN*/,1/*pt+*/,20/*pt-*/,1/*P+*/,5/*P-*/,1/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,1/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
 226
 227   fContPid = new AliPIDperfContainer("contPID",nBinPid,binPid);
 228   fContPid->SetTitleX("M_{K^{0}_{s}}");
 229   fContPid->SetTitleY("centrality (%)");
 230   fContPid->SetVarName(0,"p_{T}^{K^{0}_{s}}");
 231   fContPid->SetVarRange(0,0,10);
 232   fContPid->SetVarName(1,"#eta^{K^{0}_{s}}");
 233   fContPid->SetVarRange(1,-0.8,0.8);
 234   fContPid->SetVarName(2,"p_{T}^{Kp}");
 235   fContPid->SetVarRange(2,0.3,4.3);
 236   fContPid->SetVarName(3,"p_{T}^{Kn}");
 237   fContPid->SetVarRange(3,0.3,4.3);
 238   fContPid->SetVarName(4,"BayesProb^{Kp}");
 239   fContPid->SetVarRange(4,0,1.);
 240   fContPid->SetVarName(5,"BayesProb^{Kn}");
 241   fContPid->SetVarRange(5,0,1.);
 242   fContPid->SetVarName(6,"isTOFmatch^{Kp}");
 243   fContPid->SetVarRange(6,-0.5,1.5);
 244   fContPid->SetVarName(7,"isTOFmatch^{Kn}");
 245   fContPid->SetVarRange(7,-0.5,1.5);
 246   fContPid->SetVarName(8,"isKsTrue^{Kn}");
 247   fContPid->SetVarRange(8,-0.5,1.5);
 248   fContPid->SetVarName(9,"N#sigma^{Kp}");
 249   fContPid->SetVarRange(9,1.25,6.25);
 250   fContPid->SetVarName(10,"N#sigma^{Kn}");
 251   fContPid->SetVarRange(10,1.25,6.25);
 252   fContPid->SetVarName(11,"#Delta#phi^{Kp}");
 253   fContPid->SetVarRange(11,-TMath::Pi(),TMath::Pi());
 254   fContPid->SetVarName(12,"#Delta#phi^{Kn}");
 255   fContPid->SetVarRange(12,-TMath::Pi(),TMath::Pi());
 256   fContPid->SetVarName(13,"#Psi");
 257   fContPid->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
 258
 259   fContPid2 = new AliPIDperfContainer("contPID2",nBinPid,binPid2);
 260   fContPid2->SetTitleX("M_{K^{0}_{s}}");
 261   fContPid2->SetTitleY("centrality (%)");
 262   fContPid2->SetVarName(0,"p_{T}^{K^{0}_{s}}");
 263   fContPid2->SetVarRange(0,0,10);
 264   fContPid2->SetVarName(1,"#eta^{K^{0}_{s}}");
 265   fContPid2->SetVarRange(1,-0.8,0.8);
 266   fContPid2->SetVarName(2,"p_{T}^{Kp}");
 267   fContPid2->SetVarRange(2,0.3,4.3);
 268   fContPid2->SetVarName(3,"p_{T}^{Kn}");
 269   fContPid2->SetVarRange(3,0.3,4.3);
 270   fContPid2->SetVarName(4,"BayesProb^{Kp}");
 271   fContPid2->SetVarRange(4,0,1.);
 272   fContPid2->SetVarName(5,"BayesProb^{Kn}");
 273   fContPid2->SetVarRange(5,0,1.);
 274   fContPid2->SetVarName(6,"isTOFmatch^{Kp}");
 275   fContPid2->SetVarRange(6,-0.5,1.5);
 276   fContPid2->SetVarName(7,"isTOFmatch^{Kn}");
 277   fContPid2->SetVarRange(7,-0.5,1.5);
 278   fContPid2->SetVarName(8,"isKsTrue^{Kn}");
 279   fContPid2->SetVarRange(8,-0.5,1.5);
 280   fContPid2->SetVarName(9,"N#sigma^{Kp}");
 281   fContPid2->SetVarRange(9,1.25,6.25);
 282   fContPid2->SetVarName(10,"N#sigma^{Kn}");
 283   fContPid2->SetVarRange(10,1.25,6.25);
 284   fContPid2->SetVarName(11,"#Delta#phi^{Kp}");
 285   fContPid2->SetVarRange(11,-TMath::Pi(),TMath::Pi());
 286   fContPid2->SetVarName(12,"#Delta#phi^{Kn}");
 287   fContPid2->SetVarRange(12,-TMath::Pi(),TMath::Pi());
 288   fContPid2->SetVarName(13,"#Psi");
 289   fContPid2->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
 290
 291
 292
 293   const Int_t nDETsignal = 100; // mass
 294   Double_t binDETsignal[nDETsignal+1];
 295   for(Int_t i=0;i<nDETsignal+1;i++){
 296     binDETsignal[i] = invMmin + i*(invMmax - invMmin) / nDETsignal;
 297   }
 298   const Int_t nDETsignal2 = 10; // centrality
 299   Double_t binDETsignal2[nDETsignal2+1];
 300   for(Int_t i=0;i<nDETsignal2+1;i++){
 301     binDETsignal2[i] = i*100./ nDETsignal2;
 302   }
 303   fContPid->AddSpecies("K0s",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
 304   fContPid2->AddSpecies("K0s2",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
 305
 306   fList->Add(fContPid);
 307   fList->Add(fContPid2);
 308
 309   hMatching[0] = new TH2F("hMatchAll","TOF matched (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 310   hMatching[1] = new TH2F("hMatchPi","TOF matched (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 311   hMatching[2] = new TH2F("hMatchKa","TOF matched (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 312   hMatching[3] = new TH2F("hMatchPr","TOF matched (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 313
 314   hTracking[0] = new TH2F("hTrackingAll","TPC tracks (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 315   hTracking[1] = new TH2F("hTrackingPi","TPC tracks (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 316   hTracking[2] = new TH2F("hTrackingKa","TPC tracks (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 317   hTracking[3] = new TH2F("hTrackingPr","TPC tracks (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 318
 319   fList2->Add(hMatching[0]);
 320   fList2->Add(hMatching[1]);
 321   fList2->Add(hMatching[2]);
 322   fList2->Add(hMatching[3]);
 323   fList2->Add(hTracking[0]);
 324   fList2->Add(hTracking[1]);
 325   fList2->Add(hTracking[2]);
 326   fList2->Add(hTracking[3]);
 327
 328   fTOFTPCsignal = new TH2F("hTOFTPCsignal","TOF-TPC signal for pions (0.8 < p_{T} < 0.9 GeV/#it{c}, cent. = 0-20);N#sigma_{TOF};N#sigma_{TPC}",100,-5,5,100,-5,5);
 329   fList2->Add(fTOFTPCsignal);
 330   fCombsignal = new TH1F("hCombsignal","Combined signal for pions (0.8 < p_{T} < 0.9 GeV/#it{c}, cent. = 0-20);N#sigma",100,0,5);
 331   fList2->Add(fCombsignal);
 332
 333   // QA plots
 334   char name[100],title[400];
 335   for(Int_t i=0;i < nCentrBin;i++){
 336     snprintf(name,100,"hPiTPCQA_%i",i);
 337     snprintf(title,400,"TPC signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 338     fPiTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
 339     fList3->Add(fPiTPC[i]);
 340
 341     snprintf(name,100,"hKaTPCQA_%i",i);
 342     snprintf(title,400,"TPC signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 343     fKaTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
 344     fList3->Add(fKaTPC[i]);
 345
 346     snprintf(name,100,"hPrTPCQA_%i",i);
 347     snprintf(title,400,"TPC signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 348     fPrTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
 349     fList3->Add(fPrTPC[i]);
 350
 351     snprintf(name,100,"hElTPCQA_%i",i);
 352     snprintf(title,400,"TPC signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 353     fElTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
 354     fList3->Add(fElTPC[i]);
 355
 356     snprintf(name,100,"hPiTOFQA_%i",i);
 357     snprintf(title,400,"TOF signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 358     fPiTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
 359     fList3->Add(fPiTOF[i]);
 360
 361     snprintf(name,100,"hKaTOFQA_%i",i);
 362     snprintf(title,400,"TOF signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 363     fKaTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
 364     fList3->Add(fKaTOF[i]);
 365
 366     snprintf(name,100,"hPrTOFQA_%i",i);
 367     snprintf(title,400,"TOF signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 368     fPrTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
 369     fList3->Add(fPrTOF[i]);
 370
 371     snprintf(name,100,"hElTOFQA_%i",i);
 372     snprintf(title,400,"TOF signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 373     fElTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
 374     fList3->Add(fElTOF[i]);
 375
 376   }
 377
 378   // Post output data.
 379   PostData(1, fList);
 380   PostData(2, fList2);
 381   PostData(3, fList3);
 382 }
 383
 384 //______________________________________________________________________________
 385 void AliAnalysisTaskK0sBayes::UserExec(Option_t *)
 386 {
 387     // Main loop
 388     // Called for each event
 389
 390     fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 391     if(!fOutputAOD){
 392         Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 393         this->Dump();
 394         return;
 395     }
 396
 397     Float_t zvtx = GetVertex(fOutputAOD);
 398
 399
 400
 401     //Get the MC object
 402     if(fIsMC){
 403       AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 404       if (!mcHeader) {
 405         AliError("Could not find MC Header in AOD");
 406         return;
 407       }
 408     }
 409
 410     if (TMath::Abs(zvtx) < fVtxCut) {
 411
 412       SelectK0s();
 413
 414       //Centrality
 415       Float_t v0Centr  = -10.;
 416       Float_t trkCentr  = -10.;
 417       AliCentrality *centrality = fOutputAOD->GetCentrality();
 418       if (centrality){
 419         trkCentr  = centrality->GetCentralityPercentile("V0M");
 420         v0Centr = centrality->GetCentralityPercentile("TRK");
 421       }
 422
 423       if(!fTypeCol){
 424         v0Centr=100./(fOutputAOD->GetNumberOfTracks()/12.+1);
 425         trkCentr=v0Centr;
 426       }
 427
 428       if((TMath::Abs(v0Centr - trkCentr) < 5.0 || (fTypeCol!=2)) && v0Centr>0){ // consistency cut on centrality selection
 429         fCentrality = v0Centr;
 430         Analyze(fOutputAOD); // Do analysis!!!!
 431
 432       }
 433     }
 434
 435 }
 436
 437 //________________________________________________________________________
 438 void AliAnalysisTaskK0sBayes::Analyze(AliAODEvent* aodEvent)
 439 {
 440
 441   Int_t ntrack = aodEvent->GetNumberOfTracks();
 442
 443   fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
 444   fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
 445   fPidKp=0,fPidKn=0;
 446   fMassV0=-1;
 447
 448   TLorentzVector ksV;
 449
 450   Double_t px,py,pz,E;
 451
 452   Float_t invMmin = 0.497-0.005*8;
 453   Float_t invMmax = 0.497+0.005*8;
 454
 455   Int_t icentr = 8;
 456   if(fCentrality < 0) icentr = 8;
 457   else if(fCentrality < 10) icentr = 0;
 458   else if(fCentrality < 20) icentr = 1;
 459   else if(fCentrality < 30) icentr = 2;
 460   else if(fCentrality < 40) icentr = 3;
 461   else if(fCentrality < 50) icentr = 4;
 462   else if(fCentrality < 60) icentr = 5;
 463   else if(fCentrality < 70) icentr = 6;
 464   else if(fCentrality < 80) icentr = 7;
 465
 466   Float_t addMismatchForMC = 0.005;
 467   if(fCentrality < 50) addMismatchForMC += 0.005;
 468   if(fCentrality < 20) addMismatchForMC += 0.02;
 469
 470   if(fTypeCol == 0 || fTypeCol == 1) addMismatchForMC = 0.005;
 471
 472   fPsi = 0;
 473   /* Compute TPC EP */
 474   Double_t Qx2 = 0, Qy2 = 0;
 475   Double_t Qx3 = 0, Qy3 = 0;
 476   for(Int_t iT = 0; iT < ntrack; iT++) {
 477     AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
 478
 479     if (!aodTrack){
 480       continue;
 481     }
 482
 483     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
 484
 485     if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag)
 486       continue;
 487
 488     Double_t b[2] = {-99., -99.};
 489     Double_t bCov[3] = {-99., -99., -99.};
 490     if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
 491       continue;
 492
 493     if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
 494       continue;
 495
 496     Qx2 += TMath::Cos(2*aodTrack->Phi());
 497     Qy2 += TMath::Sin(2*aodTrack->Phi());
 498     Qx3 += TMath::Cos(3*aodTrack->Phi());
 499     Qy3 += TMath::Sin(3*aodTrack->Phi());
 500
 501   }
 502
 503   fPsi = TMath::ATan2(Qy2, Qx2)/2.;
 504
 505   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 506   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 507   AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
 508
 509   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 510
 511   Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 512   Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 513   Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
 514
 515
 516   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 517   TClonesArray *mcArray = NULL;
 518   if (mcHeader)
 519     mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 520
 521   Int_t nmc = 0;
 522   if(mcArray)
 523     nmc = mcArray->GetEntries();
 524
 525   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
 526     AliAODTrack* track = aodEvent->GetTrack(i);
 527
 528     AliAODMCParticle *mcp = NULL;
 529     Int_t pdg = 0;
 530
 531     if (!track){
 532       continue;
 533     }
 534
 535     Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
 536
 537     Int_t label = -1;
 538     if(mcArray){
 539       label = track->GetLabel();
 540       if(label != -1 && label < nmc){
 541         label = TMath::Abs(label);
 542         mcp = (AliAODMCParticle*)mcArray->At(label);
 543         pdg = TMath::Abs(mcp->GetPdgCode());
 544       }
 545       else
 546         label = -1;
 547     }
 548     else{
 549       /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
 550     }
 551
 552     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 553       hTracking[0]->Fill(track->Pt(),fCentrality);
 554       if(pdg == 211)
 555         hTracking[1]->Fill(track->Pt(),fCentrality);
 556       else if(pdg == 321)
 557         hTracking[2]->Fill(track->Pt(),fCentrality);
 558       else if(pdg == 2212)
 559         hTracking[3]->Fill(track->Pt(),fCentrality);
 560       else if(! mcp){ // Fill matching histo with the prob
 561         hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
 562         hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
 563         hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
 564       }
 565     }
 566
 567     if(!tofMatch) continue;
 568
 569     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 570       hMatching[0]->Fill(track->Pt(),fCentrality);
 571       if(pdg == 211)
 572         hMatching[1]->Fill(track->Pt(),fCentrality);
 573       else if(pdg == 321)
 574         hMatching[2]->Fill(track->Pt(),fCentrality);
 575       else if(pdg == 2212)
 576         hMatching[3]->Fill(track->Pt(),fCentrality);
 577       else if(! mcp){ // Fill matching histo with the prob
 578         hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
 579         hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
 580         hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
 581       }
 582     }
 583   }
 584
 585 //   Int_t pdg1 = -1;
 586 //   Int_t pdg2 = -1;
 587
 588
 589   // start analysis K0s
 590   for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
 591     AliAODTrack* KpTrack = aodEvent->GetTrack(i);
 592
 593     if (!KpTrack){
 594       continue;
 595     }
 596
 597     if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
 598
 599     nSigmaComb=5;
 600     nSigmaTOF = 5;
 601     fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
 602     fPidKp=0;
 603
 604     UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 605
 606     Double_t oldpP[10];
 607     fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
 608
 609     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
 610
 611     if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
 612
 613     Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
 614
 615     /*
 616     if(mcArray){
 617       Int_t labelK = TMath::Abs(KpTrack->GetLabel());
 618       AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
 619       pdg1 = TMath::Abs(mcp1->GetPdgCode());
 620     }
 621     */
 622
 623     fPidKp = Int_t(probP[2]*100);
 624
 625     if(tofMatch1){
 626       if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
 627         // remove this tof hit
 628         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 629         detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 630         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 631         fPidKp = Int_t(probP[4]*100);
 632         tofMatch1=0;
 633       }
 634       else{
 635         if(probP[2] > probP[3] && probP[2] > probP[4] && probP[2] > probP[0]) fPidKp += 128; // max prob
 636
 637         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
 638         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton))<1) fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 639         if(TMath::Abs(nSigmaTOF)<1) fPrTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton));
 640         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
 641         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron))<1) fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 642         if(TMath::Abs(nSigmaTOF)<1) fElTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron));
 643         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
 644         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon))<1) fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 645         if(TMath::Abs(nSigmaTOF)<1) fKaTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon));
 646         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
 647         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion))<1) fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 648         if(TMath::Abs(nSigmaTOF)<1) fPiTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion));
 649
 650         if(fIsMC){
 651           Float_t mismAdd = addMismatchForMC;
 652           if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
 653
 654           if(gRandom->Rndm() < mismAdd){
 655             Float_t etaAbs = TMath::Abs(KpTrack->Eta());
 656             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 657             channel = channel % 8736;
 658             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 659
 660             // generate random time
 661             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 662             Double_t times[10];
 663             KpTrack->GetIntegratedTimes(times);
 664             nSigmaTOF = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kPion);
 665           }
 666         }
 667
 668         if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
 669         nSigmaTOF = TMath::Abs(nSigmaTOF);
 670
 671         if(nSigmaTOF < 2) fPidKp += 256;
 672         else if(nSigmaTOF < 3) fPidKp += 512;
 673       }
 674     }
 675
 676     if(tofMatch1){
 677       nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
 678       if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
 679         fCombsignal->Fill(nSigmaComb);
 680       }
 681     } else {
 682       nSigmaComb = TMath::Abs(nSigmaTPC);
 683     }
 684
 685     // use sigmaTOF instead of sigmaComb
 686     if(tofMatch1) nSigmaComb = nSigmaTOF;
 687
 688     if(nSigmaComb < 2) nSigmaComb = 2;
 689     else if(nSigmaComb < 3) nSigmaComb = 3;
 690     else if(nSigmaComb < 5) nSigmaComb = 4.99;
 691     else nSigmaComb = 6;
 692
 693     Int_t iks=-1;
 694     for(Int_t k=0;k < fNK0s;k++){ // find the K0s which contains the positive track
 695       if(i == fIpiP[k]) iks = k;
 696     }
 697
 698     if(fPtKp > 4.299) fPtKp = 4.299;
 699
 700     if(iks > -1 && fIpiN[iks] > -1){
 701       //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
 702       Int_t j = fIpiN[iks];
 703       AliAODTrack* KnTrack = aodEvent->GetTrack(j);
 704
 705       if (!KnTrack){
 706         continue;
 707       }
 708
 709       if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
 710
 711       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
 712       fPidKn=0;
 713
 714       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 715       Double_t oldpN[10];
 716       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
 717
 718       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
 719
 720       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
 721
 722       nSigmaComb2=5;
 723       nSigmaTOF2=5;
 724
 725       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
 726       /*
 727       if(mcArray){
 728         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
 729         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
 730         pdg2 = TMath::Abs(mcp2->GetPdgCode());
 731       }
 732       */
 733
 734       fPidKn = Int_t(probN[2]*100);
 735
 736       if(tofMatch2){
 737         if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
 738           // remove this tof hit
 739           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 740           detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 741           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 742           fPidKn = Int_t(probN[4]*100);
 743           tofMatch2=0;
 744         }
 745         else{
 746           if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
 747
 748           nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
 749
 750           nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
 751
 752           if(fIsMC){
 753             Float_t mismAdd = addMismatchForMC;
 754             if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
 755
 756             if(gRandom->Rndm() < mismAdd){
 757               Float_t etaAbs = TMath::Abs(KnTrack->Eta());
 758               Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 759               channel = channel % 8736;
 760               Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 761
 762               // generate random time
 763               Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 764               Double_t times[10];
 765               KnTrack->GetIntegratedTimes(times);
 766               nSigmaTOF2 = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kPion);
 767             }
 768           }
 769
 770           if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
 771
 772           if(nSigmaTOF2 < 2) fPidKn += 256;
 773           else if(nSigmaTOF2 < 3) fPidKn += 512;
 774         }
 775       }
 776
 777       px = KpTrack->Px() + KnTrack->Px();
 778       py = KpTrack->Py() + KnTrack->Py();
 779       pz = KpTrack->Pz() + KnTrack->Pz();
 780       E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 1.39e-01*1.39e-01);
 781       E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
 782
 783       ksV.SetPxPyPzE(px,py,pz,E);
 784       fMassV0 = fMassKs[iks];
 785
 786       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
 787
 788
 789       fPtKs = ksV.Pt();
 790       fEtaKs = ksV.Eta();
 791       fPhiKs = ksV.Phi();
 792
 793       if(tofMatch2){
 794         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
 795         if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
 796           fCombsignal->Fill(nSigmaComb2);
 797         }
 798       } else {
 799         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
 800       }
 801
 802       // use sigmaTOF instead of sigmaComb
 803       if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
 804
 805       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
 806       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
 807       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
 808       else nSigmaComb2 = 6;
 809
 810       Bool_t isTrue = kFALSE;
 811
 812       if(mcArray){
 813         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
 814         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
 815
 816         if(labelP > -1 && labelN > -1){
 817           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
 818           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
 819
 820           Int_t mP = partP->GetMother();
 821           Int_t mN = partN->GetMother();
 822           if(mP == mN && mP > -1){
 823             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
 824             Int_t pdgM = partM->GetPdgCode();
 825             if(pdgM == 310) isTrue = kTRUE;
 826           }
 827         }
 828
 829       }
 830
 831       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
 832       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
 833
 834       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
 835       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
 836       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
 837       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
 838
 839       if(fPtKn > 4.299) fPtKn = 4.299;
 840
 841       Float_t xTOfill[] = {fPtKs,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
 842
 843
 844       Int_t ipt = 0;
 845       while(fPtKsMin[ipt] < fPtKs && ipt < nPtBin){
 846         ipt++;
 847       }
 848       ipt--;
 849       if(ipt < 0) ipt = 0; // just to be sure
 850
 851       if(TMath::Abs(fEtaKs) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
 852         fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
 853         xTOfill[1] = KnTrack->Eta();
 854         fContPid2->Fill(0,fMassV0,fCentrality,xTOfill);
 855       }
 856
 857
 858     }
 859   } // end analysi K0s
 860 }
 861
 862 //_____________________________________________________________________________
 863 Float_t AliAnalysisTaskK0sBayes::GetVertex(AliAODEvent* aod) const
 864 {
 865
 866   Float_t zvtx = -999;
 867
 868   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
 869   if (!vtxAOD)
 870     return zvtx;
 871   if(vtxAOD->GetNContributors()>0)
 872     zvtx = vtxAOD->GetZ();
 873
 874   return zvtx;
 875 }
 876 //_____________________________________________________________________________
 877 void AliAnalysisTaskK0sBayes::Terminate(Option_t *)
 878 {
 879   // Terminate loop
 880   Printf("Terminate()");
 881 }
 882 //=======================================================================
 883 void AliAnalysisTaskK0sBayes::SelectK0s(){
 884   fNK0s=0;
 885   fNpiPos=0;
 886   fNpiNeg=0;
 887
 888   Int_t nV0s = fOutputAOD->GetNumberOfV0s();
 889   AliAODv0 *myV0;
 890   Double_t dMASS=0.0;
 891   for (Int_t i=0; i!=nV0s; ++i) {
 892     myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
 893     if(!myV0) continue;
 894     if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
 895     Int_t pass = PassesAODCuts(myV0,fOutputAOD,0); // check for K0s
 896     if(pass) {
 897       dMASS = myV0->MassK0Short();
 898       Float_t massLambda = myV0->MassLambda();
 899       Float_t massAntiLambda = myV0->MassAntiLambda();
 900
 901       if(TMath::Abs(dMASS-0.497)/0.005 < 8 && TMath::Abs(massLambda-1.115)/0.005 > 8 && TMath::Abs(massAntiLambda-1.115)/0.005 > 8){
 902         AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
 903         AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
 904         if(iT->Charge()<0){
 905           iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
 906           jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
 907         }
 908         fPhiK0s[fNK0s] = myV0->Phi();
 909         fPtK0s[fNK0s] = myV0->Pt();
 910         fIpiP[fNK0s] = FindDaugheterIndex(iT);
 911         fIpiN[fNK0s] = FindDaugheterIndex(jT);
 912         fMassKs[fNK0s] = dMASS;
 913         if(fIpiP[fNK0s] > -1 && fIpiN[fNK0s] > -1)
 914           fNK0s++;
 915       }
 916     }
 917   }
 918
 919   /* My V0 code
 920   // fill pion stacks
 921   Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
 922   for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
 923     AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
 924
 925     if (!aodTrack){
 926       continue;
 927     }
 928
 929     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
 930
 931     if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
 932       continue;
 933     }
 934
 935     Double_t b[2] = {-99., -99.};
 936     Double_t bCov[3] = {-99., -99., -99.};
 937     if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
 938       continue;
 939
 940     if(TMath::Abs(b[0]) < 0.3/aodTrack->Pt()) continue;
 941
 942
 943     Int_t charge = aodTrack->Charge();
 944     if(charge > 0){
 945       fIPiPos[fNpiPos] = iT;
 946       fNpiPos++;
 947     }
 948     else{
 949       fIPiNeg[fNpiNeg] = iT;
 950       fNpiNeg++;
 951     }
 952   }
 953
 954   for(Int_t i=0;i < fNpiPos;i++){
 955     AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
 956     AliESDtrack pipE(pip);
 957
 958     for(Int_t j=0;j < fNpiNeg;j++){
 959       AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
 960       AliESDtrack pinE(pin);
 961
 962       Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
 963
 964       Double_t pPos[3];
 965       Double_t pNeg[3];
 966       pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
 967       pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
 968
 969       Float_t length = (xp+xn)*0.5;
 970
 971       Float_t pxs = pPos[0] + pNeg[0];
 972       Float_t pys = pPos[1] + pNeg[1];
 973       Float_t pzs = pPos[2] + pNeg[2];
 974       Float_t es = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.13957*0.13957) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.13957*0.13957);
 975
 976       Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
 977       Float_t phi = TMath::ATan2(pys,pxs);
 978       Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
 979
 980       //      if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
 981
 982       if(mindist < 0.4&& length > 0.7 && length < 25){
 983         Float_t esL = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.938*0.938) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.13957*0.13957);
 984         Float_t esAL = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.13957*0.13957) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.938*0.938);
 985
 986         Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
 987         Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
 988
 989         if(TMath::Abs(mass-0.497)/0.005 < 8 && massaL > 1.15 && massaAL > 1.15){
 990           fPhiK0s[fNK0s] = phi;
 991           fPtK0s[fNK0s] = pt;
 992           fIpiP[fNK0s] =fIPiPos[i] ;
 993           fIpiN[fNK0s] = fIPiNeg[j];
 994           fMassKs[fNK0s] = mass;
 995           fNK0s++;
 996         }
 997       }
 998     }
 999   }
1000   */
1001 }
1002
1003 //=======================================================================
1004 Int_t AliAnalysisTaskK0sBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1005 {
1006   if (myV0->GetOnFlyStatus() ) return 0;
1007
1008   //the following is needed in order to evualuate track-quality
1009   AliAODTrack *iT, *jT;
1010   AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1011   Double_t pos[3],cov[6];
1012   vV0s->GetXYZ(pos);
1013   vV0s->GetCovarianceMatrix(cov);
1014   const AliESDVertex vESD(pos,cov,100.,100);
1015
1016   // TESTING CHARGE
1017   int iPos, iNeg;
1018   iT=(AliAODTrack*) myV0->GetDaughter(0);
1019   if(iT->Charge()>0) {
1020     iPos = 0; iNeg = 1;
1021   } else {
1022     iPos = 1; iNeg = 0;
1023   }
1024   // END OF TEST
1025
1026   iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1027
1028   jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1029
1030   Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1031   if(!trkFlag) return 0;
1032   Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1033   if(!trkFlag2) return 0;
1034
1035   Double_t pvertex[3];
1036   pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1037   pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1038   pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1039
1040   Double_t dDL=myV0->DecayLengthV0( pvertex );
1041   if(dDL  < 0.5 || dDL > 25) return 0;
1042
1043   Double_t dDCA=myV0->DcaV0Daughters();
1044   if(dDCA >0.5) return 0;
1045
1046   Double_t dCTP=myV0->CosPointingAngle( pvertex );
1047   if(dCTP < -1) return 0;
1048
1049 //   AliESDtrack ieT( iT );
1050 //   Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1051 //   if(TMath::Abs(dD0P) < 0]) return 0;
1052
1053 //   AliESDtrack jeT( jT );
1054 //   Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1055 //   if(TMath::Abs(dD0M) < 0) return 0;
1056
1057 //   Double_t dD0D0=dD0P*dD0M;
1058 //   if(dD0D0>0) return 0;
1059
1060 //   Double_t dETA=myV0->Eta();
1061 //   if(dETA <-0.8) return 0;
1062 //   if(dETA >0.8) return 0;
1063
1064 //   Double_t dQT=myV0->PtArmV0();
1065 //   if(specie==0) if(dQT<???) return 0;
1066
1067   Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1068   if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1069
1070   if(specie==1 && dALPHA<0) return 2; // antilambda
1071   return 1; // K0s or lambda
1072 }
1073 //-------------------------------------------------------------------------
1074 Int_t AliAnalysisTaskK0sBayes::FindDaugheterIndex(AliAODTrack *trk){
1075   Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1076
1077   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1078     AliAODTrack* track = fOutputAOD->GetTrack(i);
1079     if(track == trk) return i;
1080   }
1081
1082   printf("Daughter for %p not found\n",trk);
1083   return -1;
1084 }
1085 //-------------------------------------------------------------------------
1086 Int_t AliAnalysisTaskK0sBayes::IsChannelValid(Float_t etaAbs){
1087   if(!fIsMC) return 1; // used only on MC
1088
1089   if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
1090     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1091
1092     if(!(channel%20)) return 0; // 5% additional loss in MC
1093   }
1094
1095   return 1;
1096 }