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((TMath::Abs(v0Centr - trkCentr) < 5.0 || (fTypeCol!=2)) && v0Centr>0){ // consistency cut on centrality selection
 424         fCentrality = v0Centr;
 425         Analyze(fOutputAOD); // Do analysis!!!!
 426
 427       }
 428     }
 429
 430 }
 431
 432 //________________________________________________________________________
 433 void AliAnalysisTaskK0sBayes::Analyze(AliAODEvent* aodEvent)
 434 {
 435
 436   Int_t ntrack = aodEvent->GetNumberOfTracks();
 437
 438   fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
 439   fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
 440   fPidKp=0,fPidKn=0;
 441   fMassV0=-1;
 442
 443   TLorentzVector ksV;
 444
 445   Double_t px,py,pz,E;
 446
 447   Float_t invMmin = 0.497-0.005*8;
 448   Float_t invMmax = 0.497+0.005*8;
 449
 450   Int_t icentr = 8;
 451   if(fCentrality < 0) icentr = 8;
 452   else if(fCentrality < 10) icentr = 0;
 453   else if(fCentrality < 20) icentr = 1;
 454   else if(fCentrality < 30) icentr = 2;
 455   else if(fCentrality < 40) icentr = 3;
 456   else if(fCentrality < 50) icentr = 4;
 457   else if(fCentrality < 60) icentr = 5;
 458   else if(fCentrality < 70) icentr = 6;
 459   else if(fCentrality < 80) icentr = 7;
 460
 461   Float_t addMismatchForMC = 0.005;
 462   if(fCentrality < 50) addMismatchForMC += 0.01;
 463   if(fCentrality < 20) addMismatchForMC += 0.02;
 464
 465   if(fTypeCol == 0) addMismatchForMC = 0.005;
 466   else if(fTypeCol == 1) addMismatchForMC *= 0.5;
 467
 468   fPsi = 0;
 469   /* Compute TPC EP */
 470   Double_t Qx2 = 0, Qy2 = 0;
 471   Double_t Qx3 = 0, Qy3 = 0;
 472   for(Int_t iT = 0; iT < ntrack; iT++) {
 473     AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
 474
 475     if (!aodTrack){
 476       continue;
 477     }
 478
 479     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
 480
 481     if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag)
 482       continue;
 483
 484     Double_t b[2] = {-99., -99.};
 485     Double_t bCov[3] = {-99., -99., -99.};
 486     if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
 487       continue;
 488
 489     if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
 490       continue;
 491
 492     Qx2 += TMath::Cos(2*aodTrack->Phi());
 493     Qy2 += TMath::Sin(2*aodTrack->Phi());
 494     Qx3 += TMath::Cos(3*aodTrack->Phi());
 495     Qy3 += TMath::Sin(3*aodTrack->Phi());
 496
 497   }
 498
 499   fPsi = TMath::ATan2(Qy2, Qx2)/2.;
 500
 501   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 502   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 503   AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
 504
 505   PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
 506   PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
 507   PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
 508   PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
 509
 510   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 511
 512   Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 513   Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 514   Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
 515
 516
 517   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 518   TClonesArray *mcArray = NULL;
 519   if (mcHeader)
 520     mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 521
 522   Int_t nmc = 0;
 523   if(mcArray)
 524     nmc = mcArray->GetEntries();
 525
 526   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
 527     AliAODTrack* track = aodEvent->GetTrack(i);
 528
 529     AliAODMCParticle *mcp = NULL;
 530     Int_t pdg = 0;
 531
 532     if (!track){
 533       continue;
 534     }
 535
 536     Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
 537
 538     Int_t label = -1;
 539     if(mcArray){
 540       label = track->GetLabel();
 541       if(label != -1 && label < nmc){
 542         label = TMath::Abs(label);
 543         mcp = (AliAODMCParticle*)mcArray->At(label);
 544         pdg = TMath::Abs(mcp->GetPdgCode());
 545       }
 546       else
 547         label = -1;
 548     }
 549     else{
 550       /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
 551     }
 552
 553     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 554       hTracking[0]->Fill(track->Pt(),fCentrality);
 555       if(pdg == 211)
 556         hTracking[1]->Fill(track->Pt(),fCentrality);
 557       else if(pdg == 321)
 558         hTracking[2]->Fill(track->Pt(),fCentrality);
 559       else if(pdg == 2212)
 560         hTracking[3]->Fill(track->Pt(),fCentrality);
 561       else if(! mcp){ // Fill matching histo with the prob
 562         hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
 563         hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
 564         hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
 565       }
 566     }
 567
 568     if(!tofMatch) continue;
 569
 570     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 571       hMatching[0]->Fill(track->Pt(),fCentrality);
 572       if(pdg == 211)
 573         hMatching[1]->Fill(track->Pt(),fCentrality);
 574       else if(pdg == 321)
 575         hMatching[2]->Fill(track->Pt(),fCentrality);
 576       else if(pdg == 2212)
 577         hMatching[3]->Fill(track->Pt(),fCentrality);
 578       else if(! mcp){ // Fill matching histo with the prob
 579         hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
 580         hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
 581         hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
 582       }
 583     }
 584   }
 585
 586 //   Int_t pdg1 = -1;
 587 //   Int_t pdg2 = -1;
 588
 589
 590   // start analysis K0s
 591   for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
 592     AliAODTrack* KpTrack = aodEvent->GetTrack(i);
 593
 594     if (!KpTrack){
 595       continue;
 596     }
 597
 598     if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
 599
 600     nSigmaComb=5;
 601     nSigmaTOF = 5;
 602     fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
 603     fPidKp=0;
 604
 605     UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 606
 607     Double_t oldpP[10];
 608     fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
 609
 610     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon);
 611     fKaTPC[icentr]->Fill(fPtKp,nSigmaTPC);
 612     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton);
 613     fPrTPC[icentr]->Fill(fPtKp,nSigmaTPC);
 614     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron);
 615     fElTPC[icentr]->Fill(fPtKp,nSigmaTPC);
 616
 617     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
 618     fPiTPC[icentr]->Fill(fPtKp,nSigmaTPC);
 619
 620     if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
 621
 622     Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
 623
 624     /*
 625     if(mcArray){
 626       Int_t labelK = TMath::Abs(KpTrack->GetLabel());
 627       AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
 628       pdg1 = TMath::Abs(mcp1->GetPdgCode());
 629     }
 630     */
 631
 632     fPidKp = Int_t(probP[2]*100);
 633
 634     if(tofMatch1){
 635       if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
 636         // remove this tof hit
 637         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 638         detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 639         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 640         fPidKp = Int_t(probP[4]*100);
 641         tofMatch1=0;
 642       }
 643       else{
 644         if(probP[2] > probP[3] && probP[2] > probP[4] && probP[2] > probP[0]) fPidKp += 128; // max prob
 645
 646         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
 647         fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 648         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
 649         fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 650
 651         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
 652         fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 653
 654         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
 655         fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 656
 657         nSigmaTOF = TMath::Abs(nSigmaTOF);
 658
 659         if(fIsMC){
 660           Float_t mismAdd = addMismatchForMC;
 661           if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
 662
 663           if(gRandom->Rndm() < mismAdd){
 664             Float_t etaAbs = TMath::Abs(KpTrack->Eta());
 665             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 666             channel = channel % 8736;
 667             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 668
 669             // generate random time
 670             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 671             Double_t times[10];
 672             KpTrack->GetIntegratedTimes(times);
 673             nSigmaTOF = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kPion);
 674           }
 675         }
 676
 677         if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
 678
 679         if(nSigmaTOF < 2) fPidKp += 256;
 680         else if(nSigmaTOF < 3) fPidKp += 512;
 681       }
 682     }
 683
 684     if(tofMatch1){
 685       nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
 686       if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
 687         fCombsignal->Fill(nSigmaComb);
 688       }
 689     } else {
 690       nSigmaComb = TMath::Abs(nSigmaTPC);
 691     }
 692
 693     // use sigmaTOF instead of sigmaComb
 694     if(tofMatch1) nSigmaComb = nSigmaTOF;
 695
 696     if(nSigmaComb < 2) nSigmaComb = 2;
 697     else if(nSigmaComb < 3) nSigmaComb = 3;
 698     else if(nSigmaComb < 5) nSigmaComb = 4.99;
 699     else nSigmaComb = 6;
 700
 701     Int_t iks=-1;
 702     for(Int_t k=0;k < fNK0s;k++){ // find the K0s which contains the positive track
 703       if(i == fIpiP[k]) iks = k;
 704     }
 705
 706     if(fPtKp > 4.299) fPtKp = 4.299;
 707
 708     if(iks > -1 && fIpiN[iks] > -1){
 709       //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
 710       Int_t j = fIpiN[iks];
 711       AliAODTrack* KnTrack = aodEvent->GetTrack(j);
 712
 713       if (!KnTrack){
 714         continue;
 715       }
 716
 717       if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
 718
 719       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
 720       fPidKn=0;
 721
 722       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 723       Double_t oldpN[10];
 724       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
 725
 726       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
 727
 728       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
 729
 730       nSigmaComb2=5;
 731       nSigmaTOF2=5;
 732
 733       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
 734       /*
 735       if(mcArray){
 736         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
 737         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
 738         pdg2 = TMath::Abs(mcp2->GetPdgCode());
 739       }
 740       */
 741
 742       fPidKn = Int_t(probN[2]*100);
 743
 744       if(tofMatch2){
 745         if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
 746           // remove this tof hit
 747           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 748           detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 749           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 750           fPidKn = Int_t(probN[4]*100);
 751           tofMatch2=0;
 752         }
 753         else{
 754           if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
 755
 756           nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
 757
 758           nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
 759
 760           if(fIsMC){
 761             Float_t mismAdd = addMismatchForMC;
 762             if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
 763
 764             if(gRandom->Rndm() < mismAdd){
 765               Float_t etaAbs = TMath::Abs(KnTrack->Eta());
 766               Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 767               channel = channel % 8736;
 768               Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 769
 770               // generate random time
 771               Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 772               Double_t times[10];
 773               KnTrack->GetIntegratedTimes(times);
 774               nSigmaTOF2 = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kPion);
 775             }
 776           }
 777
 778           if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
 779
 780           if(nSigmaTOF2 < 2) fPidKn += 256;
 781           else if(nSigmaTOF2 < 3) fPidKn += 512;
 782         }
 783       }
 784
 785       px = KpTrack->Px() + KnTrack->Px();
 786       py = KpTrack->Py() + KnTrack->Py();
 787       pz = KpTrack->Pz() + KnTrack->Pz();
 788       E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 1.39e-01*1.39e-01);
 789       E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
 790
 791       ksV.SetPxPyPzE(px,py,pz,E);
 792       fMassV0 = fMassKs[iks];
 793
 794       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
 795
 796
 797       fPtKs = ksV.Pt();
 798       fEtaKs = ksV.Eta();
 799       fPhiKs = ksV.Phi();
 800
 801       if(tofMatch2){
 802         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
 803         if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
 804           fCombsignal->Fill(nSigmaComb2);
 805         }
 806       } else {
 807         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
 808       }
 809
 810       // use sigmaTOF instead of sigmaComb
 811       if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
 812
 813       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
 814       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
 815       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
 816       else nSigmaComb2 = 6;
 817
 818       Bool_t isTrue = kFALSE;
 819
 820       if(mcArray){
 821         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
 822         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
 823
 824         if(labelP > -1 && labelN > -1){
 825           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
 826           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
 827
 828           Int_t mP = partP->GetMother();
 829           Int_t mN = partN->GetMother();
 830           if(mP == mN && mP > -1){
 831             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
 832             Int_t pdgM = partM->GetPdgCode();
 833             if(pdgM == 310) isTrue = kTRUE;
 834           }
 835         }
 836
 837       }
 838
 839       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
 840       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
 841
 842       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
 843       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
 844       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
 845       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
 846
 847       if(fPtKn > 4.299) fPtKn = 4.299;
 848
 849       Float_t xTOfill[] = {fPtKs,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
 850
 851
 852       Int_t ipt = 0;
 853       while(fPtKsMin[ipt] < fPtKs && ipt < nPtBin){
 854         ipt++;
 855       }
 856       ipt--;
 857       if(ipt < 0) ipt = 0; // just to be sure
 858
 859       if(TMath::Abs(fEtaKs) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
 860         fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
 861         xTOfill[1] = KnTrack->Eta();
 862         fContPid2->Fill(0,fMassV0,fCentrality,xTOfill);
 863       }
 864
 865
 866     }
 867   } // end analysi K0s
 868 }
 869
 870 //_____________________________________________________________________________
 871 Float_t AliAnalysisTaskK0sBayes::GetVertex(AliAODEvent* aod) const
 872 {
 873
 874   Float_t zvtx = -999;
 875
 876   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
 877   if (!vtxAOD)
 878     return zvtx;
 879   if(vtxAOD->GetNContributors()>0)
 880     zvtx = vtxAOD->GetZ();
 881
 882   return zvtx;
 883 }
 884 //_____________________________________________________________________________
 885 void AliAnalysisTaskK0sBayes::Terminate(Option_t *)
 886 {
 887   // Terminate loop
 888   Printf("Terminate()");
 889 }
 890 //=======================================================================
 891 void AliAnalysisTaskK0sBayes::SelectK0s(){
 892   fNK0s=0;
 893   fNpiPos=0;
 894   fNpiNeg=0;
 895
 896   Int_t nV0s = fOutputAOD->GetNumberOfV0s();
 897   AliAODv0 *myV0;
 898   Double_t dMASS=0.0;
 899   for (Int_t i=0; i!=nV0s; ++i) {
 900     myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
 901     if(!myV0) continue;
 902     if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
 903     Int_t pass = PassesAODCuts(myV0,fOutputAOD,0); // check for K0s
 904     if(pass) {
 905       dMASS = myV0->MassK0Short();
 906       Float_t massLambda = myV0->MassLambda();
 907       Float_t massAntiLambda = myV0->MassAntiLambda();
 908
 909       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){
 910         AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
 911         AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
 912         if(iT->Charge()<0){
 913           iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
 914           jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
 915         }
 916         fPhiK0s[fNK0s] = myV0->Phi();
 917         fPtK0s[fNK0s] = myV0->Pt();
 918         fIpiP[fNK0s] = FindDaugheterIndex(iT);
 919         fIpiN[fNK0s] = FindDaugheterIndex(jT);
 920         fMassKs[fNK0s] = dMASS;
 921         if(fIpiP[fNK0s] > -1 && fIpiN[fNK0s] > -1)
 922           fNK0s++;
 923       }
 924     }
 925   }
 926
 927   /* My V0 code
 928   // fill pion stacks
 929   Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
 930   for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
 931     AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
 932
 933     if (!aodTrack){
 934       continue;
 935     }
 936
 937     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
 938
 939     if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
 940       continue;
 941     }
 942
 943     Double_t b[2] = {-99., -99.};
 944     Double_t bCov[3] = {-99., -99., -99.};
 945     if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
 946       continue;
 947
 948     if(TMath::Abs(b[0]) < 0.3/aodTrack->Pt()) continue;
 949
 950
 951     Int_t charge = aodTrack->Charge();
 952     if(charge > 0){
 953       fIPiPos[fNpiPos] = iT;
 954       fNpiPos++;
 955     }
 956     else{
 957       fIPiNeg[fNpiNeg] = iT;
 958       fNpiNeg++;
 959     }
 960   }
 961
 962   for(Int_t i=0;i < fNpiPos;i++){
 963     AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
 964     AliESDtrack pipE(pip);
 965
 966     for(Int_t j=0;j < fNpiNeg;j++){
 967       AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
 968       AliESDtrack pinE(pin);
 969
 970       Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
 971
 972       Double_t pPos[3];
 973       Double_t pNeg[3];
 974       pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
 975       pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
 976
 977       Float_t length = (xp+xn)*0.5;
 978
 979       Float_t pxs = pPos[0] + pNeg[0];
 980       Float_t pys = pPos[1] + pNeg[1];
 981       Float_t pzs = pPos[2] + pNeg[2];
 982       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);
 983
 984       Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
 985       Float_t phi = TMath::ATan2(pys,pxs);
 986       Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
 987
 988       //      if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
 989
 990       if(mindist < 0.4&& length > 0.7 && length < 25){
 991         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);
 992         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);
 993
 994         Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
 995         Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
 996
 997         if(TMath::Abs(mass-0.497)/0.005 < 8 && massaL > 1.15 && massaAL > 1.15){
 998           fPhiK0s[fNK0s] = phi;
 999           fPtK0s[fNK0s] = pt;
1000           fIpiP[fNK0s] =fIPiPos[i] ;
1001           fIpiN[fNK0s] = fIPiNeg[j];
1002           fMassKs[fNK0s] = mass;
1003           fNK0s++;
1004         }
1005       }
1006     }
1007   }
1008   */
1009 }
1010
1011 //=======================================================================
1012 Int_t AliAnalysisTaskK0sBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1013 {
1014   if (myV0->GetOnFlyStatus() ) return 0;
1015
1016   //the following is needed in order to evualuate track-quality
1017   AliAODTrack *iT, *jT;
1018   AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1019   Double_t pos[3],cov[6];
1020   vV0s->GetXYZ(pos);
1021   vV0s->GetCovarianceMatrix(cov);
1022   const AliESDVertex vESD(pos,cov,100.,100);
1023
1024   // TESTING CHARGE
1025   int iPos, iNeg;
1026   iT=(AliAODTrack*) myV0->GetDaughter(0);
1027   if(iT->Charge()>0) {
1028     iPos = 0; iNeg = 1;
1029   } else {
1030     iPos = 1; iNeg = 0;
1031   }
1032   // END OF TEST
1033
1034   iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1035
1036   jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1037
1038   Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1039   if(!trkFlag) return 0;
1040   Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1041   if(!trkFlag2) return 0;
1042
1043   Double_t pvertex[3];
1044   pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1045   pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1046   pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1047
1048   Double_t dDL=myV0->DecayLengthV0( pvertex );
1049   if(dDL  < 0.5 || dDL > 25) return 0;
1050
1051   Double_t dDCA=myV0->DcaV0Daughters();
1052   if(dDCA >0.5) return 0;
1053
1054   Double_t dCTP=myV0->CosPointingAngle( pvertex );
1055   if(dCTP < -1) return 0;
1056
1057 //   AliESDtrack ieT( iT );
1058 //   Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1059 //   if(TMath::Abs(dD0P) < 0]) return 0;
1060
1061 //   AliESDtrack jeT( jT );
1062 //   Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1063 //   if(TMath::Abs(dD0M) < 0) return 0;
1064
1065 //   Double_t dD0D0=dD0P*dD0M;
1066 //   if(dD0D0>0) return 0;
1067
1068 //   Double_t dETA=myV0->Eta();
1069 //   if(dETA <-0.8) return 0;
1070 //   if(dETA >0.8) return 0;
1071
1072 //   Double_t dQT=myV0->PtArmV0();
1073 //   if(specie==0) if(dQT<???) return 0;
1074
1075   Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1076   if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1077
1078   if(specie==1 && dALPHA<0) return 2; // antilambda
1079   return 1; // K0s or lambda
1080 }
1081 //-------------------------------------------------------------------------
1082 Int_t AliAnalysisTaskK0sBayes::FindDaugheterIndex(AliAODTrack *trk){
1083   Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1084
1085   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1086     AliAODTrack* track = fOutputAOD->GetTrack(i);
1087     if(track == trk) return i;
1088   }
1089
1090   printf("Daughter for %p not found\n",trk);
1091   return -1;
1092 }
1093 //-------------------------------------------------------------------------
1094 Int_t AliAnalysisTaskK0sBayes::IsChannelValid(Float_t etaAbs){
1095   if(!fIsMC) return 1; // used only on MC
1096
1097   if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
1098     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1099
1100     if(!(channel%20)) return 0; // 5% additional loss in MC
1101   }
1102
1103   return 1;
1104 }