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