PWGPP/pid/AliAnalysisTaskPhiBayes.cxx

   1 #include "AliAnalysisTaskPhiBayes.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 AliAnalysisTaskPhiBayes::fPtPhiMin[AliAnalysisTaskPhiBayes::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 AliAnalysisTaskPhiBayes::fPtPhiMax[AliAnalysisTaskPhiBayes::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(AliAnalysisTaskPhiBayes)
  39 //_____________________________________________________________________________
  40 AliAnalysisTaskPhiBayes::AliAnalysisTaskPhiBayes():
  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   fPtPhi(0.),
  55   fPhiPhi(0.),
  56   fEtaPhi(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   fHmismTOF(0),
  73   fHchannelTOFdistr(0),
  74   fTypeCol(2)
  75 {
  76   // Default constructor (should not be used)
  77   fList->SetName("contPhiBayes1");
  78   fList2->SetName("contPhiBayes2");
  79   fList3->SetName("contPhiBayes3");
  80
  81   fList->SetOwner(kTRUE);
  82   fList2->SetOwner(kTRUE);
  83   fList3->SetOwner(kTRUE);
  84
  85   TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
  86   fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
  87
  88   TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
  89   fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
  90
  91   for(Int_t i=0;i < nCentrBin;i++){
  92     fElTOF[i] = NULL;
  93     fElTPC[i] = NULL;
  94     fPiTOF[i] = NULL;
  95     fPiTPC[i] = NULL;
  96     fKaTOF[i] = NULL;
  97     fKaTPC[i] = NULL;
  98     fPrTOF[i] = NULL;
  99     fPrTPC[i] = NULL;
 100   }
 101   for(Int_t i=0;i < 4;i++){
 102     hMatching[i] = NULL;
 103     hTracking[i] = NULL;
 104   }
 105 }
 106
 107 //______________________________________________________________________________
 108 AliAnalysisTaskPhiBayes::AliAnalysisTaskPhiBayes(const char *name):
 109   AliAnalysisTaskSE(name),
 110   fVtxCut(10.0),  // cut on |vertex| < fVtxCut
 111   fEtaCut(0.8),   // cut on |eta| < fEtaCut
 112   fMinPt(0.15),   // cut on pt > fMinPt
 113   fIsMC(kFALSE),
 114   fQAsw(kFALSE),
 115   fNcluster(70),
 116   fFilterBit(1),
 117   fList(new TList()),
 118   fList2(new TList()),
 119   fList3(new TList()),
 120   fCentrality(-1),
 121   fPsi(0),
 122   fPtPhi(0.),
 123   fPhiPhi(0.),
 124   fEtaPhi(0.),
 125   fMassV0(0.),
 126   fPtKp(0.),
 127   fPhiKp(0.),
 128   fEtaKp(0.),
 129   fPtKn(0.),
 130   fPhiKn(0.),
 131   fEtaKn(0.),
 132   fPidKp(0),
 133   fPidKn(0),
 134   fTOFTPCsignal(0),
 135   fCombsignal(0),
 136   fCutsDaughter(NULL),
 137   fPIDCombined(NULL),
 138   fContPid(NULL),
 139   fContPid2(NULL),
 140   fHmismTOF(0),
 141   fHchannelTOFdistr(0),
 142   fTypeCol(2)
 143 {
 144
 145   DefineOutput(1, TList::Class());
 146   DefineOutput(2, TList::Class());
 147   DefineOutput(3, TList::Class());
 148   DefineOutput(4, TList::Class());
 149
 150   // Output slot #1 writes into a TTree
 151   fList->SetName("contPhiBayes1");
 152   fList2->SetName("contPhiBayes2");
 153   fList3->SetName("contPhiBayes3");
 154
 155   fList->SetOwner(kTRUE);
 156   fList2->SetOwner(kTRUE);
 157   fList3->SetOwner(kTRUE);
 158
 159   TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
 160   fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
 161
 162   TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
 163   fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
 164
 165   for(Int_t i=0;i < nCentrBin;i++){
 166     fElTOF[i] = NULL;
 167     fElTPC[i] = NULL;
 168     fPiTOF[i] = NULL;
 169     fPiTPC[i] = NULL;
 170     fKaTOF[i] = NULL;
 171     fKaTPC[i] = NULL;
 172     fPrTOF[i] = NULL;
 173     fPrTPC[i] = NULL;
 174   }
 175   for(Int_t i=0;i < 4;i++){
 176     hMatching[i] = NULL;
 177     hTracking[i] = NULL;
 178   }
 179 }
 180 //_____________________________________________________________________________
 181 AliAnalysisTaskPhiBayes::~AliAnalysisTaskPhiBayes()
 182 {
 183
 184 }
 185
 186 //______________________________________________________________________________
 187 void AliAnalysisTaskPhiBayes::UserCreateOutputObjects()
 188 {
 189
 190   fPIDCombined=new AliPIDCombined;
 191   fPIDCombined->SetDefaultTPCPriors();
 192   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
 193
 194   Float_t invMmin = 0.985;
 195   Float_t invMmax = 1.045;
 196
 197   const Int_t nBinPid = 14;
 198
 199   Int_t binPid[nBinPid] = {1/*ptPhi*/,8/*EtaKP*/,20/*pt+*/,1/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,1/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,1/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
 200
 201   Int_t binPid2[nBinPid] = {1/*ptPhi*/,8/*EtaKN*/,1/*pt+*/,20/*pt-*/,1/*P+*/,5/*P-*/,1/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,1/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
 202
 203   fContPid = new AliPIDperfContainer("contPID",nBinPid,binPid);
 204   fContPid->SetTitleX("M_{#phi}");
 205   fContPid->SetTitleY("centrality (%)");
 206   fContPid->SetVarName(0,"p_{T}^#phi}");
 207   fContPid->SetVarRange(0,0,10);
 208   fContPid->SetVarName(1,"#eta^{#phi}");
 209   fContPid->SetVarRange(1,-0.8,0.8);
 210   fContPid->SetVarName(2,"p_{T}^{Kp}");
 211   fContPid->SetVarRange(2,0.3,4.3);
 212   fContPid->SetVarName(3,"p_{T}^{Kn}");
 213   fContPid->SetVarRange(3,0.3,4.3);
 214   fContPid->SetVarName(4,"BayesProb^{Kp}");
 215   fContPid->SetVarRange(4,0,1.);
 216   fContPid->SetVarName(5,"BayesProb^{Kn}");
 217   fContPid->SetVarRange(5,0,1.);
 218   fContPid->SetVarName(6,"isTOFmatch^{Kp}");
 219   fContPid->SetVarRange(6,-0.5,1.5);
 220   fContPid->SetVarName(7,"isTOFmatch^{Kn}");
 221   fContPid->SetVarRange(7,-0.5,1.5);
 222   fContPid->SetVarName(8,"isPhiTrue^{Kn}");
 223   fContPid->SetVarRange(8,-0.5,1.5);
 224   fContPid->SetVarName(9,"N#sigma^{Kp}");
 225   fContPid->SetVarRange(9,1.25,6.25);
 226   fContPid->SetVarName(10,"N#sigma^{Kn}");
 227   fContPid->SetVarRange(10,1.25,6.25);
 228   fContPid->SetVarName(11,"#Delta#phi^{Kp}");
 229   fContPid->SetVarRange(11,-TMath::Pi(),TMath::Pi());
 230   fContPid->SetVarName(12,"#Delta#phi^{Kn}");
 231   fContPid->SetVarRange(12,-TMath::Pi(),TMath::Pi());
 232   fContPid->SetVarName(13,"#Psi");
 233   fContPid->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
 234
 235   fContPid2 = new AliPIDperfContainer("contPID2",nBinPid,binPid2);
 236   fContPid2->SetTitleX("M_{#phi}");
 237   fContPid2->SetTitleY("centrality (%)");
 238   fContPid2->SetVarName(0,"p_{T}^{#phi}");
 239   fContPid2->SetVarRange(0,0,10);
 240   fContPid2->SetVarName(1,"#eta^{#phi}");
 241   fContPid2->SetVarRange(1,-0.8,0.8);
 242   fContPid2->SetVarName(2,"p_{T}^{Kp}");
 243   fContPid2->SetVarRange(2,0.3,4.3);
 244   fContPid2->SetVarName(3,"p_{T}^{Kn}");
 245   fContPid2->SetVarRange(3,0.3,4.3);
 246   fContPid2->SetVarName(4,"BayesProb^{Kp}");
 247   fContPid2->SetVarRange(4,0,1.);
 248   fContPid2->SetVarName(5,"BayesProb^{Kn}");
 249   fContPid2->SetVarRange(5,0,1.);
 250   fContPid2->SetVarName(6,"isTOFmatch^{Kp}");
 251   fContPid2->SetVarRange(6,-0.5,1.5);
 252   fContPid2->SetVarName(7,"isTOFmatch^{Kn}");
 253   fContPid2->SetVarRange(7,-0.5,1.5);
 254   fContPid2->SetVarName(8,"isPhiTrue^{Kn}");
 255   fContPid2->SetVarRange(8,-0.5,1.5);
 256   fContPid2->SetVarName(9,"N#sigma^{Kp}");
 257   fContPid2->SetVarRange(9,1.25,6.25);
 258   fContPid2->SetVarName(10,"N#sigma^{Kn}");
 259   fContPid2->SetVarRange(10,1.25,6.25);
 260   fContPid2->SetVarName(11,"#Delta#phi^{Kp}");
 261   fContPid2->SetVarRange(11,-TMath::Pi(),TMath::Pi());
 262   fContPid2->SetVarName(12,"#Delta#phi^{Kn}");
 263   fContPid2->SetVarRange(12,-TMath::Pi(),TMath::Pi());
 264   fContPid2->SetVarName(13,"#Psi");
 265   fContPid2->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
 266
 267
 268
 269   const Int_t nDETsignal = 100; // mass
 270   Double_t binDETsignal[nDETsignal+1];
 271   for(Int_t i=0;i<nDETsignal+1;i++){
 272     binDETsignal[i] = invMmin + i*(invMmax - invMmin) / nDETsignal;
 273   }
 274   const Int_t nDETsignal2 = 10; // centrality
 275   Double_t binDETsignal2[nDETsignal2+1];
 276   for(Int_t i=0;i<nDETsignal2+1;i++){
 277     binDETsignal2[i] = i*100./ nDETsignal2;
 278   }
 279   fContPid->AddSpecies("phi",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
 280   fContPid2->AddSpecies("phi2",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
 281
 282   fList->Add(fContPid);
 283   fList->Add(fContPid2);
 284
 285   hMatching[0] = new TH2F("hMatchAll","TOF matched (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 286   hMatching[1] = new TH2F("hMatchPi","TOF matched (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 287   hMatching[2] = new TH2F("hMatchKa","TOF matched (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 288   hMatching[3] = new TH2F("hMatchPr","TOF matched (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 289
 290   hTracking[0] = new TH2F("hTrackingAll","TPC tracks (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 291   hTracking[1] = new TH2F("hTrackingPi","TPC tracks (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 292   hTracking[2] = new TH2F("hTrackingKa","TPC tracks (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 293   hTracking[3] = new TH2F("hTrackingPr","TPC tracks (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
 294
 295   fList2->Add(hMatching[0]);
 296   fList2->Add(hMatching[1]);
 297   fList2->Add(hMatching[2]);
 298   fList2->Add(hMatching[3]);
 299   fList2->Add(hTracking[0]);
 300   fList2->Add(hTracking[1]);
 301   fList2->Add(hTracking[2]);
 302   fList2->Add(hTracking[3]);
 303
 304   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);
 305   fList2->Add(fTOFTPCsignal);
 306   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);
 307   fList2->Add(fCombsignal);
 308
 309   // QA plots
 310   char name[100],title[400];
 311   for(Int_t i=0;i < nCentrBin;i++){
 312     snprintf(name,100,"hPiTPCQA_%i",i);
 313     snprintf(title,400,"TPC signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 314     fPiTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
 315     fList3->Add(fPiTPC[i]);
 316
 317     snprintf(name,100,"hKaTPCQA_%i",i);
 318     snprintf(title,400,"TPC signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 319     fKaTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
 320     fList3->Add(fKaTPC[i]);
 321
 322     snprintf(name,100,"hPrTPCQA_%i",i);
 323     snprintf(title,400,"TPC signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 324     fPrTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
 325     fList3->Add(fPrTPC[i]);
 326
 327     snprintf(name,100,"hElTPCQA_%i",i);
 328     snprintf(title,400,"TPC signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 329     fElTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
 330     fList3->Add(fElTPC[i]);
 331
 332     snprintf(name,100,"hPiTOFQA_%i",i);
 333     snprintf(title,400,"TOF signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 334     fPiTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
 335     fList3->Add(fPiTOF[i]);
 336
 337     snprintf(name,100,"hKaTOFQA_%i",i);
 338     snprintf(title,400,"TOF signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 339     fKaTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
 340     fList3->Add(fKaTOF[i]);
 341
 342     snprintf(name,100,"hPrTOFQA_%i",i);
 343     snprintf(title,400,"TOF signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 344     fPrTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
 345     fList3->Add(fPrTOF[i]);
 346
 347     snprintf(name,100,"hElTOFQA_%i",i);
 348     snprintf(title,400,"TOF signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
 349     fElTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
 350     fList3->Add(fElTOF[i]);
 351
 352   }
 353
 354   // Post output data.
 355   PostData(1, fList);
 356   PostData(2, fList2);
 357   PostData(3, fList3);
 358 }
 359
 360 //______________________________________________________________________________
 361 void AliAnalysisTaskPhiBayes::UserExec(Option_t *)
 362 {
 363     // Main loop
 364     // Called for each event
 365
 366     fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 367     if(!fOutputAOD){
 368         Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 369         this->Dump();
 370         return;
 371     }
 372
 373     Float_t zvtx = GetVertex(fOutputAOD);
 374
 375
 376
 377     //Get the MC object
 378     if(fIsMC){
 379       AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 380       if (!mcHeader) {
 381         AliError("Could not find MC Header in AOD");
 382         return;
 383       }
 384     }
 385
 386     if (TMath::Abs(zvtx) < fVtxCut) {
 387
 388       //Centrality
 389       Float_t v0Centr  = -10.;
 390       Float_t trkCentr  = -10.;
 391       AliCentrality *centrality = fOutputAOD->GetCentrality();
 392       if (centrality){
 393         trkCentr  = centrality->GetCentralityPercentile("V0M");
 394         v0Centr = centrality->GetCentralityPercentile("TRK");
 395       }
 396
 397       if(!fTypeCol){
 398         v0Centr=100./(fOutputAOD->GetNumberOfTracks()/12.+1);
 399         trkCentr=v0Centr;
 400       }
 401
 402       if((TMath::Abs(v0Centr - trkCentr) < 5.0 || (fTypeCol!=2)) && v0Centr>0){ // consistency cut on centrality selection
 403         fCentrality = v0Centr;
 404         Analyze(fOutputAOD); // Do analysis!!!!
 405
 406       }
 407     }
 408
 409 }
 410
 411 //________________________________________________________________________
 412 void AliAnalysisTaskPhiBayes::Analyze(AliAODEvent* aodEvent)
 413 {
 414   Int_t ntrack = aodEvent->GetNumberOfTracks();
 415
 416   fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
 417   fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
 418   fPidKp=0,fPidKn=0;
 419   fMassV0=-1;
 420
 421   TLorentzVector phiV;
 422
 423   Double_t px,py,pz,E;
 424
 425   Float_t invMmin = 0.985;
 426   Float_t invMmax = 1.045;
 427
 428   Int_t icentr = 8;
 429   if(fCentrality < 0) icentr = 8;
 430   else if(fCentrality < 10) icentr = 0;
 431   else if(fCentrality < 20) icentr = 1;
 432   else if(fCentrality < 30) icentr = 2;
 433   else if(fCentrality < 40) icentr = 3;
 434   else if(fCentrality < 50) icentr = 4;
 435   else if(fCentrality < 60) icentr = 5;
 436   else if(fCentrality < 70) icentr = 6;
 437   else if(fCentrality < 80) icentr = 7;
 438
 439   Float_t addMismatchForMC = 0.005;
 440   if(fCentrality < 50) addMismatchForMC += 0.005;
 441   if(fCentrality < 20) addMismatchForMC += 0.02;
 442
 443   if(fTypeCol == 0 || fTypeCol == 1) addMismatchForMC = 0.005;
 444
 445   fPsi = 0;
 446   /* Compute TPC EP */
 447   Double_t Qx2 = 0, Qy2 = 0;
 448   Double_t Qx3 = 0, Qy3 = 0;
 449   for(Int_t iT = 0; iT < ntrack; iT++) {
 450     AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
 451
 452     if (!aodTrack){
 453       continue;
 454     }
 455
 456     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
 457
 458     if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag)
 459       continue;
 460
 461     Double_t b[2] = {-99., -99.};
 462     Double_t bCov[3] = {-99., -99., -99.};
 463     if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
 464       continue;
 465
 466     if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
 467       continue;
 468
 469     Qx2 += TMath::Cos(2*aodTrack->Phi());
 470     Qy2 += TMath::Sin(2*aodTrack->Phi());
 471     Qx3 += TMath::Cos(3*aodTrack->Phi());
 472     Qy3 += TMath::Sin(3*aodTrack->Phi());
 473
 474   }
 475
 476   fPsi = TMath::ATan2(Qy2, Qx2)/2.;
 477
 478   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 479   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 480   AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
 481
 482 //   PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
 483 //   PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
 484 //   PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
 485 //   PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
 486
 487   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 488
 489   Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 490   Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 491   Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
 492
 493
 494   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 495   TClonesArray *mcArray = NULL;
 496   if (mcHeader)
 497     mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 498
 499   Int_t nmc = 0;
 500   if(mcArray)
 501     nmc = mcArray->GetEntries();
 502
 503   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
 504     AliAODTrack* track = aodEvent->GetTrack(i);
 505
 506     AliAODMCParticle *mcp = NULL;
 507     Int_t pdg = 0;
 508
 509     if (!track){
 510       continue;
 511     }
 512
 513     Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
 514
 515     Int_t label = -1;
 516     if(mcArray){
 517       label = track->GetLabel();
 518       if(label != -1 && label < nmc){
 519         label = TMath::Abs(label);
 520         mcp = (AliAODMCParticle*)mcArray->At(label);
 521         pdg = TMath::Abs(mcp->GetPdgCode());
 522       }
 523       else
 524         label = -1;
 525     }
 526     else{
 527       /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
 528     }
 529
 530     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 531       hTracking[0]->Fill(track->Pt(),fCentrality);
 532       if(pdg == 211)
 533         hTracking[1]->Fill(track->Pt(),fCentrality);
 534       else if(pdg == 321)
 535         hTracking[2]->Fill(track->Pt(),fCentrality);
 536       else if(pdg == 2212)
 537         hTracking[3]->Fill(track->Pt(),fCentrality);
 538       else if(! mcp){ // Fill matching histo with the prob
 539         hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
 540         hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
 541         hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
 542       }
 543     }
 544
 545     if(!tofMatch) continue;
 546
 547     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 548       hMatching[0]->Fill(track->Pt(),fCentrality);
 549       if(pdg == 211)
 550         hMatching[1]->Fill(track->Pt(),fCentrality);
 551       else if(pdg == 321)
 552         hMatching[2]->Fill(track->Pt(),fCentrality);
 553       else if(pdg == 2212)
 554         hMatching[3]->Fill(track->Pt(),fCentrality);
 555       else if(! mcp){ // Fill matching histo with the prob
 556         hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
 557         hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
 558         hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
 559       }
 560     }
 561   }
 562
 563 //   Int_t pdg1 = -1;
 564 //   Int_t pdg2 = -1;
 565
 566
 567   // start analysis phi
 568   for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
 569     AliAODTrack* KpTrack = aodEvent->GetTrack(i);
 570
 571     if (!KpTrack){
 572       continue;
 573     }
 574
 575     if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
 576
 577     nSigmaTOF=5;
 578     nSigmaComb=5;
 579     fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
 580     fPidKp=0;
 581
 582     UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 583
 584     Double_t oldpP[10];
 585     fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
 586
 587     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon);
 588
 589     if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
 590
 591     Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
 592     /*
 593     if(mcArray){
 594       Int_t labelK = TMath::Abs(KpTrack->GetLabel());
 595       AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
 596       pdg1 = TMath::Abs(mcp1->GetPdgCode());
 597     }
 598     */
 599
 600     fPidKp = Int_t(probP[3]*100);
 601
 602     if(tofMatch1){
 603       if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
 604         // remove this tof hit
 605         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 606         detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 607         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 608         fPidKp = Int_t(probP[4]*100);
 609         tofMatch1=0;
 610       }
 611       else{
 612         if(probP[3] > probP[2] && probP[3] > probP[4] && probP[3] > probP[0]) fPidKp += 128; // max prob
 613
 614         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
 615         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton))<1) fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 616         if(TMath::Abs(nSigmaTOF)<1) fPrTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton));
 617         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
 618         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron))<1) fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 619         if(TMath::Abs(nSigmaTOF)<1) fElTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron));
 620         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
 621         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion))<1) fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 622         if(TMath::Abs(nSigmaTOF)<1) fPiTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion));
 623         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
 624         if(TMath::Abs(PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon))<1) fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 625         if(TMath::Abs(nSigmaTOF)<1) fKaTPC[icentr]->Fill(fPtKp,PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon));
 626
 627         if(fIsMC){
 628           Float_t mismAdd = addMismatchForMC;
 629           if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
 630
 631           if(gRandom->Rndm() < mismAdd){
 632             Float_t etaAbs = TMath::Abs(KpTrack->Eta());
 633             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 634             channel = channel % 8736;
 635             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 636
 637             // generate random time
 638             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 639             Double_t times[10];
 640             KpTrack->GetIntegratedTimes(times);
 641             nSigmaTOF = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kKaon);
 642           }
 643         }
 644
 645         if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
 646         nSigmaTOF = TMath::Abs(nSigmaTOF);
 647
 648         if(nSigmaTOF < 2) fPidKp += 256;
 649         else if(nSigmaTOF < 3) fPidKp += 512;
 650       }
 651     }
 652
 653     if(tofMatch1){
 654       nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
 655       if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
 656         fCombsignal->Fill(nSigmaComb);
 657       }
 658     } else {
 659       nSigmaComb = TMath::Abs(nSigmaTPC);
 660     }
 661
 662     // use sigmaTOF instead of sigmaComb
 663     if(tofMatch1) nSigmaComb = nSigmaTOF;
 664
 665     if(nSigmaComb < 2) nSigmaComb = 2;
 666     else if(nSigmaComb < 3) nSigmaComb = 3;
 667     else if(nSigmaComb < 5) nSigmaComb = 4.99;
 668     else nSigmaComb = 6;
 669
 670     if(fPtKp > 4.299) fPtKp = 4.299;
 671
 672     for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
 673       AliAODTrack* KnTrack = aodEvent->GetTrack(j);
 674
 675       if (!KnTrack){
 676         continue;
 677       }
 678
 679       if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
 680
 681       px = KpTrack->Px() + KnTrack->Px();
 682       py = KpTrack->Py() + KnTrack->Py();
 683       pz = KpTrack->Pz() + KnTrack->Pz();
 684       E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 4.93676999999999977e-01*4.93676999999999977e-01);
 685       E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 4.93676999999999977e-01*4.93676999999999977e-01);
 686
 687       phiV.SetPxPyPzE(px,py,pz,E);
 688       fMassV0 = phiV.M();
 689
 690       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
 691
 692       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
 693       fPidKn=0;
 694
 695       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 696       Double_t oldpN[10];
 697       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
 698
 699       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kKaon);
 700
 701       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
 702
 703       nSigmaTOF2=5;
 704       nSigmaComb2=5;
 705
 706       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
 707       /*
 708       if(mcArray){
 709         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
 710         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
 711         pdg2 = TMath::Abs(mcp2->GetPdgCode());
 712       }
 713       */
 714
 715       fPidKn = Int_t(probN[3]*100);
 716
 717       if(tofMatch2){
 718         if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
 719           // remove this tof hit
 720           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 721           detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 722           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 723           fPidKn = Int_t(probN[4]*100);
 724           tofMatch2=0;
 725         }
 726         else{
 727           if(probN[3] > probN[2] && probN[3] > probN[4] && probN[3] > probN[0]) fPidKn += 128; // max prob
 728
 729           nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kKaon);
 730
 731           nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
 732
 733           if(fIsMC){
 734             Float_t mismAdd = addMismatchForMC;
 735             if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
 736
 737             if(gRandom->Rndm() < mismAdd){
 738               Float_t etaAbs = TMath::Abs(KnTrack->Eta());
 739               Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 740               channel = channel % 8736;
 741               Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 742
 743               // generate random time
 744               Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 745               Double_t times[10];
 746               KnTrack->GetIntegratedTimes(times);
 747               nSigmaTOF2 = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kKaon);
 748             }
 749           }
 750
 751           if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
 752
 753           if(nSigmaTOF2 < 2) fPidKn += 256;
 754           else if(nSigmaTOF2 < 3) fPidKn += 512;
 755         }
 756       }
 757
 758       fPtPhi = phiV.Pt();
 759       fEtaPhi = phiV.Eta();
 760       fPhiPhi = phiV.Phi();
 761
 762       if(tofMatch2){
 763         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
 764         if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
 765           fCombsignal->Fill(nSigmaComb2);
 766         }
 767       } else {
 768         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
 769       }
 770
 771       // use sigmaTOF instead of sigmaComb
 772       if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
 773
 774       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
 775       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
 776       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
 777       else nSigmaComb2 = 6;
 778
 779       Bool_t isTrue = kFALSE;
 780
 781       if(mcArray){
 782         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
 783         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
 784
 785         if(labelP > -1 && labelN > -1){
 786           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
 787           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
 788
 789           Int_t mP = partP->GetMother();
 790           Int_t mN = partN->GetMother();
 791           if(mP == mN && mP > -1){
 792             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
 793             Int_t pdgM = partM->GetPdgCode();
 794             if(pdgM == 333) isTrue = kTRUE;
 795           }
 796         }
 797
 798       }
 799
 800       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
 801       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
 802
 803       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
 804       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
 805       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
 806       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
 807
 808       if(fPtKn > 4.299) fPtKn = 4.299;
 809
 810       Float_t xTOfill[] = {fPtPhi,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
 811
 812
 813       Int_t ipt = 0;
 814       while(fPtPhiMin[ipt] < fPtPhi && ipt < nPtBin){
 815         ipt++;
 816       }
 817       ipt--;
 818       if(ipt < 0) ipt = 0; // just to be sure
 819
 820       if(TMath::Abs(fEtaPhi) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
 821         if((fPidKn%128) > 80) fContPid->Fill(0,fMassV0,fCentrality,xTOfill); // use tagging on negative track
 822         xTOfill[1] = KnTrack->Eta();
 823         if((fPidKp%128) > 80) fContPid2->Fill(0,fMassV0,fCentrality,xTOfill);// use tagging on positive track
 824
 825       }
 826
 827
 828     }
 829   } // end analysi phi
 830
 831 }
 832
 833 //_____________________________________________________________________________
 834 Float_t AliAnalysisTaskPhiBayes::GetVertex(AliAODEvent* aod) const
 835 {
 836
 837   Float_t zvtx = -999;
 838
 839   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
 840   if (!vtxAOD)
 841     return zvtx;
 842   if(vtxAOD->GetNContributors()>0)
 843     zvtx = vtxAOD->GetZ();
 844
 845   return zvtx;
 846 }
 847 //_____________________________________________________________________________
 848 void AliAnalysisTaskPhiBayes::Terminate(Option_t *)
 849 {
 850   // Terminate loop
 851   Printf("Terminate()");
 852 }
 853 //-------------------------------------------------------------------------
 854 Int_t AliAnalysisTaskPhiBayes::IsChannelValid(Float_t etaAbs){
 855   if(!fIsMC) return 1; // used only on MC
 856
 857   if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
 858     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 859
 860     if(!(channel%20)) return 0; // 5% additional loss in MC
 861   }
 862
 863   return 1;
 864 }