PWGPP/pid/AliAnalysisTaskLambdaBayes.cxx

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