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((TMath::Abs(v0Centr - trkCentr) < 5.0 || (fTypeCol!=2)) && v0Centr>0){ // consistency cut on centrality selection
 398         fCentrality = v0Centr;
 399         Analyze(fOutputAOD); // Do analysis!!!!
 400
 401       }
 402     }
 403
 404 }
 405
 406 //________________________________________________________________________
 407 void AliAnalysisTaskPhiBayes::Analyze(AliAODEvent* aodEvent)
 408 {
 409   Int_t ntrack = aodEvent->GetNumberOfTracks();
 410
 411   fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
 412   fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
 413   fPidKp=0,fPidKn=0;
 414   fMassV0=-1;
 415
 416   TLorentzVector phiV;
 417
 418   Double_t px,py,pz,E;
 419
 420   Float_t invMmin = 0.985;
 421   Float_t invMmax = 1.045;
 422
 423   Int_t icentr = 8;
 424   if(fCentrality < 0) icentr = 8;
 425   else if(fCentrality < 10) icentr = 0;
 426   else if(fCentrality < 20) icentr = 1;
 427   else if(fCentrality < 30) icentr = 2;
 428   else if(fCentrality < 40) icentr = 3;
 429   else if(fCentrality < 50) icentr = 4;
 430   else if(fCentrality < 60) icentr = 5;
 431   else if(fCentrality < 70) icentr = 6;
 432   else if(fCentrality < 80) icentr = 7;
 433
 434   Float_t addMismatchForMC = 0.005;
 435   if(fCentrality < 50) addMismatchForMC += 0.005;
 436   if(fCentrality < 20) addMismatchForMC += 0.02;
 437
 438   if(fTypeCol == 0 || fTypeCol == 1) addMismatchForMC = 0.005;
 439
 440   fPsi = 0;
 441   /* Compute TPC EP */
 442   Double_t Qx2 = 0, Qy2 = 0;
 443   Double_t Qx3 = 0, Qy3 = 0;
 444   for(Int_t iT = 0; iT < ntrack; iT++) {
 445     AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
 446
 447     if (!aodTrack){
 448       continue;
 449     }
 450
 451     Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
 452
 453     if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag)
 454       continue;
 455
 456     Double_t b[2] = {-99., -99.};
 457     Double_t bCov[3] = {-99., -99., -99.};
 458     if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
 459       continue;
 460
 461     if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
 462       continue;
 463
 464     Qx2 += TMath::Cos(2*aodTrack->Phi());
 465     Qy2 += TMath::Sin(2*aodTrack->Phi());
 466     Qx3 += TMath::Cos(3*aodTrack->Phi());
 467     Qy3 += TMath::Sin(3*aodTrack->Phi());
 468
 469   }
 470
 471   fPsi = TMath::ATan2(Qy2, Qx2)/2.;
 472
 473   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 474   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 475   AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
 476
 477   PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
 478   PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
 479   PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
 480   PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
 481
 482   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 483
 484   Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 485   Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 486   Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
 487
 488
 489   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 490   TClonesArray *mcArray = NULL;
 491   if (mcHeader)
 492     mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 493
 494   Int_t nmc = 0;
 495   if(mcArray)
 496     nmc = mcArray->GetEntries();
 497
 498   for(Int_t i=0;i < ntrack;i++){ // loop on tracks
 499     AliAODTrack* track = aodEvent->GetTrack(i);
 500
 501     AliAODMCParticle *mcp = NULL;
 502     Int_t pdg = 0;
 503
 504     if (!track){
 505       continue;
 506     }
 507
 508     Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
 509
 510     Int_t label = -1;
 511     if(mcArray){
 512       label = track->GetLabel();
 513       if(label != -1 && label < nmc){
 514         label = TMath::Abs(label);
 515         mcp = (AliAODMCParticle*)mcArray->At(label);
 516         pdg = TMath::Abs(mcp->GetPdgCode());
 517       }
 518       else
 519         label = -1;
 520     }
 521     else{
 522       /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
 523     }
 524
 525     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 526       hTracking[0]->Fill(track->Pt(),fCentrality);
 527       if(pdg == 211)
 528         hTracking[1]->Fill(track->Pt(),fCentrality);
 529       else if(pdg == 321)
 530         hTracking[2]->Fill(track->Pt(),fCentrality);
 531       else if(pdg == 2212)
 532         hTracking[3]->Fill(track->Pt(),fCentrality);
 533       else if(! mcp){ // Fill matching histo with the prob
 534         hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
 535         hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
 536         hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
 537       }
 538     }
 539
 540     if(!tofMatch) continue;
 541
 542     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
 543       hMatching[0]->Fill(track->Pt(),fCentrality);
 544       if(pdg == 211)
 545         hMatching[1]->Fill(track->Pt(),fCentrality);
 546       else if(pdg == 321)
 547         hMatching[2]->Fill(track->Pt(),fCentrality);
 548       else if(pdg == 2212)
 549         hMatching[3]->Fill(track->Pt(),fCentrality);
 550       else if(! mcp){ // Fill matching histo with the prob
 551         hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
 552         hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
 553         hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
 554       }
 555     }
 556   }
 557
 558 //   Int_t pdg1 = -1;
 559 //   Int_t pdg2 = -1;
 560
 561
 562   // start analysis phi
 563   for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
 564     AliAODTrack* KpTrack = aodEvent->GetTrack(i);
 565
 566     if (!KpTrack){
 567       continue;
 568     }
 569
 570     if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3  && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
 571
 572     nSigmaTOF=5;
 573     nSigmaComb=5;
 574     fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
 575     fPidKp=0;
 576
 577     UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 578
 579     Double_t oldpP[10];
 580     fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
 581
 582     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton);
 583     fPrTPC[icentr]->Fill(fPtKp,nSigmaTPC);
 584     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron);
 585     fElTPC[icentr]->Fill(fPtKp,nSigmaTPC);
 586     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
 587     fPiTPC[icentr]->Fill(fPtKp,nSigmaTPC);
 588     nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon);
 589     fKaTPC[icentr]->Fill(fPtKp,nSigmaTPC);
 590
 591     if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
 592
 593     Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
 594     /*
 595     if(mcArray){
 596       Int_t labelK = TMath::Abs(KpTrack->GetLabel());
 597       AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
 598       pdg1 = TMath::Abs(mcp1->GetPdgCode());
 599     }
 600     */
 601
 602     fPidKp = Int_t(probP[3]*100);
 603
 604     if(tofMatch1){
 605       if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
 606         // remove this tof hit
 607         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 608         detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
 609         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 610         fPidKp = Int_t(probP[4]*100);
 611         tofMatch1=0;
 612       }
 613       else{
 614         if(probP[3] > probP[2] && probP[3] > probP[4] && probP[3] > probP[0]) fPidKp += 128; // max prob
 615
 616         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
 617         fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 618         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
 619         fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 620         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
 621         fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 622         nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
 623         fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
 624
 625         if(fIsMC){
 626           Float_t mismAdd = addMismatchForMC;
 627           if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
 628
 629           if(gRandom->Rndm() < mismAdd){
 630             Float_t etaAbs = TMath::Abs(KpTrack->Eta());
 631             Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 632             channel = channel % 8736;
 633             Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 634
 635             // generate random time
 636             Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 637             Double_t times[10];
 638             KpTrack->GetIntegratedTimes(times);
 639             nSigmaTOF = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kKaon);
 640           }
 641         }
 642
 643         if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
 644         nSigmaTOF = TMath::Abs(nSigmaTOF);
 645
 646         if(nSigmaTOF < 2) fPidKp += 256;
 647         else if(nSigmaTOF < 3) fPidKp += 512;
 648       }
 649     }
 650
 651     if(tofMatch1){
 652       nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
 653       if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
 654         fCombsignal->Fill(nSigmaComb);
 655       }
 656     } else {
 657       nSigmaComb = TMath::Abs(nSigmaTPC);
 658     }
 659
 660     // use sigmaTOF instead of sigmaComb
 661     if(tofMatch1) nSigmaComb = nSigmaTOF;
 662
 663     if(nSigmaComb < 2) nSigmaComb = 2;
 664     else if(nSigmaComb < 3) nSigmaComb = 3;
 665     else if(nSigmaComb < 5) nSigmaComb = 4.99;
 666     else nSigmaComb = 6;
 667
 668     if(fPtKp > 4.299) fPtKp = 4.299;
 669
 670     for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
 671       AliAODTrack* KnTrack = aodEvent->GetTrack(j);
 672
 673       if (!KnTrack){
 674         continue;
 675       }
 676
 677       if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
 678
 679       px = KpTrack->Px() + KnTrack->Px();
 680       py = KpTrack->Py() + KnTrack->Py();
 681       pz = KpTrack->Pz() + KnTrack->Pz();
 682       E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 4.93676999999999977e-01*4.93676999999999977e-01);
 683       E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 4.93676999999999977e-01*4.93676999999999977e-01);
 684
 685       phiV.SetPxPyPzE(px,py,pz,E);
 686       fMassV0 = phiV.M();
 687
 688       if(fMassV0 <  invMmin || fMassV0 > invMmax) continue;
 689
 690       fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
 691       fPidKn=0;
 692
 693       UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 694       Double_t oldpN[10];
 695       fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
 696
 697       nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kKaon);
 698
 699       if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
 700
 701       nSigmaTOF2=5;
 702       nSigmaComb2=5;
 703
 704       Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
 705       /*
 706       if(mcArray){
 707         Int_t labelK = TMath::Abs(KnTrack->GetLabel());
 708         AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
 709         pdg2 = TMath::Abs(mcp2->GetPdgCode());
 710       }
 711       */
 712
 713       fPidKn = Int_t(probN[3]*100);
 714
 715       if(tofMatch2){
 716         if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
 717           // remove this tof hit
 718           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 719           detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
 720           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 721           fPidKn = Int_t(probN[4]*100);
 722           tofMatch2=0;
 723         }
 724         else{
 725           if(probN[3] > probN[2] && probN[3] > probN[4] && probN[3] > probN[0]) fPidKn += 128; // max prob
 726
 727           nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kKaon);
 728
 729           nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
 730
 731           if(fIsMC){
 732             Float_t mismAdd = addMismatchForMC;
 733             if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
 734
 735             if(gRandom->Rndm() < mismAdd){
 736               Float_t etaAbs = TMath::Abs(KnTrack->Eta());
 737               Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 738               channel = channel % 8736;
 739               Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 740
 741               // generate random time
 742               Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
 743               Double_t times[10];
 744               KnTrack->GetIntegratedTimes(times);
 745               nSigmaTOF2 = TMath::Abs(timeRandom - times[3])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kKaon);
 746             }
 747           }
 748
 749           if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
 750
 751           if(nSigmaTOF2 < 2) fPidKn += 256;
 752           else if(nSigmaTOF2 < 3) fPidKn += 512;
 753         }
 754       }
 755
 756       fPtPhi = phiV.Pt();
 757       fEtaPhi = phiV.Eta();
 758       fPhiPhi = phiV.Phi();
 759
 760       if(tofMatch2){
 761         nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
 762         if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
 763           fCombsignal->Fill(nSigmaComb2);
 764         }
 765       } else {
 766         nSigmaComb2 = TMath::Abs(nSigmaTPC2);
 767       }
 768
 769       // use sigmaTOF instead of sigmaComb
 770       if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
 771
 772       if(nSigmaComb2 < 2) nSigmaComb2 = 2;
 773       else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
 774       else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
 775       else nSigmaComb2 = 6;
 776
 777       Bool_t isTrue = kFALSE;
 778
 779       if(mcArray){
 780         Int_t labelP = TMath::Abs(KpTrack->GetLabel());
 781         Int_t labelN = TMath::Abs(KnTrack->GetLabel());
 782
 783         if(labelP > -1 && labelN > -1){
 784           AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
 785           AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
 786
 787           Int_t mP = partP->GetMother();
 788           Int_t mN = partN->GetMother();
 789           if(mP == mN && mP > -1){
 790             AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
 791             Int_t pdgM = partM->GetPdgCode();
 792             if(pdgM == 333) isTrue = kTRUE;
 793           }
 794         }
 795
 796       }
 797
 798       Double_t deltaphi1 = KpTrack->Phi() - fPsi;
 799       Double_t deltaphi2 = KnTrack->Phi() - fPsi;
 800
 801       while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
 802       while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
 803       while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
 804       while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
 805
 806       if(fPtKn > 4.299) fPtKn = 4.299;
 807
 808       Float_t xTOfill[] = {fPtPhi,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
 809
 810
 811       Int_t ipt = 0;
 812       while(fPtPhiMin[ipt] < fPtPhi && ipt < nPtBin){
 813         ipt++;
 814       }
 815       ipt--;
 816       if(ipt < 0) ipt = 0; // just to be sure
 817
 818       if(TMath::Abs(fEtaPhi) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
 819         if((fPidKn%128) > 80) fContPid->Fill(0,fMassV0,fCentrality,xTOfill); // use tagging on negative track
 820         xTOfill[1] = KnTrack->Eta();
 821         if((fPidKp%128) > 80) fContPid2->Fill(0,fMassV0,fCentrality,xTOfill);// use tagging on positive track
 822
 823       }
 824
 825
 826     }
 827   } // end analysi phi
 828
 829 }
 830
 831 //_____________________________________________________________________________
 832 Float_t AliAnalysisTaskPhiBayes::GetVertex(AliAODEvent* aod) const
 833 {
 834
 835   Float_t zvtx = -999;
 836
 837   const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
 838   if (!vtxAOD)
 839     return zvtx;
 840   if(vtxAOD->GetNContributors()>0)
 841     zvtx = vtxAOD->GetZ();
 842
 843   return zvtx;
 844 }
 845 //_____________________________________________________________________________
 846 void AliAnalysisTaskPhiBayes::Terminate(Option_t *)
 847 {
 848   // Terminate loop
 849   Printf("Terminate()");
 850 }
 851 //-------------------------------------------------------------------------
 852 Int_t AliAnalysisTaskPhiBayes::IsChannelValid(Float_t etaAbs){
 853   if(!fIsMC) return 1; // used only on MC
 854
 855   if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
 856     Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 857
 858     if(!(channel%20)) return 0; // 5% additional loss in MC
 859   }
 860
 861   return 1;
 862 }