PWGLF/STRANGENESS/LambdaK0PbPb/AliV0CutVariations.C

   1 #define AliV0CutVariations_cxx
   2 // The class definition in AliV0CutVariations.h has been generated automatically
   3 // by the ROOT utility TTree::MakeSelector(). This class is derived
   4 // from the ROOT class TSelector. For more information on the TSelector
   5 // framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
   6
   7 // The following methods are defined in this file:
   8 //    Begin():        called every time a loop on the tree starts,
   9 //                    a convenient place to create your histograms.
  10 //    SlaveBegin():   called after Begin(), when on PROOF called only on the
  11 //                    slave servers.
  12 //    Process():      called for each event, in this function you decide what
  13 //                    to read and fill your histograms.
  14 //    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
  15 //                    called only on the slave servers.
  16 //    Terminate():    called at the end of the loop on the tree,
  17 //                    a convenient place to draw/fit your histograms.
  18 //
  19 // To use this file, try the following session on your Tree T:
  20 //
  21 // Root > T->Process("AliV0CutVariations.C")
  22 // Root > T->Process("AliV0CutVariations.C","some options")
  23 // Root > T->Process("AliV0CutVariations.C+")
  24 //
  25
  26 #include "AliV0CutVariations.h"
  27
  28 #include "Riostream.h"
  29 #include "TCanvas.h"
  30 #include "TPDGCode.h"
  31 #include <TDatabasePDG.h>
  32 #include <TH2.h>
  33 #include <TStyle.h>
  34 #include <TH1F.h>
  35 #include <TH2F.h>
  36 #include <TH3F.h>
  37
  38 AliV0CutVariations::AliV0CutVariations(TTree * /*tree*/ ) :
  39 fChain(0),
  40 fIsMC(kFALSE),
  41 fSelectNonInjected(kFALSE),
  42 fCMin(0.),
  43 fCMax(90.),
  44 fCPA(0.9975),
  45 fDCA(1.0),
  46 fTPCcr(70.),
  47 fTPCcrfd(0.8),
  48 fDCApv(0.1),
  49 fMult(0),
  50 fdEdx(0),
  51 fdEdxPid(0),
  52 fCosPA(0),
  53 fDtrDCA(0),
  54 fTPCrows(0),
  55 fTPCratio(0),
  56 fPrimDCA(0),
  57
  58 fK0sM(0),
  59 fK0sSi(0),
  60 fK0sMC(0),
  61 fK0sAs(0),
  62
  63 fLambdaM(0),
  64 fLambdaSi(0),
  65 fLambdaMC(0),
  66 fLambdaAs(0),
  67
  68 fLambdaEff(0),
  69 fLambdaPt(0),
  70
  71 fLambdaBarM(0),
  72 fLambdaBarSi(0),
  73 fLambdaBarMC(0),
  74 fLambdaBarAs(0),
  75
  76 fLambdaBarEff(0),
  77 fLambdaBarPt(0),
  78
  79 fLambdaFromXi(0),
  80 fXiM(0),
  81 fXiSiP(0),
  82
  83 fLambdaBarFromXiBar(0),
  84 fXiBarM(0),
  85 fXiBarSiP(0)
  86 { }
  87
  88 void AliV0CutVariations::Begin(TTree * /*tree*/)
  89 {
  90    // The Begin() function is called at the start of the query.
  91    // When running with PROOF Begin() is only called on the client.
  92    // The tree argument is deprecated (on PROOF 0 is passed).
  93
  94    TString option = GetOption();
  95
  96 }
  97
  98 void AliV0CutVariations::SlaveBegin(TTree * tree)
  99 {
 100    // The SlaveBegin() function is called after the Begin() function.
 101    // When running with PROOF SlaveBegin() is called on each slave server.
 102    // The tree argument is deprecated (on PROOF 0 is passed).
 103
 104   Init(tree);
 105
 106   //cout<<tree->GetEntries()<<endl;
 107
 108   TString option = GetOption();
 109
 110   Int_t    lbins=100;  // number of bins in lt
 111   Int_t    kbins=33;  // number of bins in pt of Xi
 112   Double_t ltMax=100.;
 113   Double_t ptMax=12.;
 114
 115   const Double_t xBins[]={
 116     0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,
 117     1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,
 118     2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,
 119     4.5,5.0,5.5,6.5,8.0,10.0,12.0
 120   };
 121   const Int_t nBins=sizeof(xBins)/sizeof(Double_t) - 1;
 122
 123   Double_t yBins[lbins+1];
 124   for (Int_t i=0; i<=(lbins+1); i++) yBins[i]=i*ltMax/lbins;
 125   Double_t zBins[kbins+1];
 126   for (Int_t i=0; i<=(kbins+1); i++) zBins[i]=i*(ptMax+2)/kbins;
 127
 128
 129   //fMult=new TH1F("fMult","Multiplicity",1100,0.,3300);
 130   //fMult->GetXaxis()->SetTitle("N tracks");
 131   fMult=new TH1F("fMult","Centrality",100,0.,100);
 132   fMult->GetXaxis()->SetTitle("Centrality (%)");
 133   fOutput->Add(fMult);
 134
 135
 136   fCosPA=new TH1F("fCosPA","Cos(PA) distribution",50,0.9975,1.0005);
 137   fCosPA->GetXaxis()->SetTitle("Cos(PA)");
 138   fOutput->Add(fCosPA);
 139
 140   fDtrDCA=new TH1F("fDtrDCA","DCA between V0 daughters",50,0.0,1.5);
 141   fDtrDCA->GetXaxis()->SetTitle("DCA (rel. u.)");
 142   fOutput->Add(fDtrDCA);
 143
 144   fTPCrows=new TH1F("fTPCrows","TPC crossed pad rows",180,0.,180.);
 145   fTPCrows->GetXaxis()->SetTitle("TPC crossed pad rows");
 146   fOutput->Add(fTPCrows);
 147
 148   fTPCratio=new TH1F("fTPCratio","TPC crossed/findable pad rows",50,0.0,1.5);
 149   fTPCratio->GetXaxis()->SetTitle("TPC crossed/findable pad rows");
 150   fOutput->Add(fTPCratio);
 151
 152   fPrimDCA=new TH1F("fPrimDCA","DCA wrt the primary vertex",50,0.0,1.5);
 153   fPrimDCA->GetXaxis()->SetTitle("DCA wrt the PV (cm)");
 154   fOutput->Add(fPrimDCA);
 155
 156
 157   fdEdx=new TH2F("fdEdx","dE/dx",50,0.2,3,50,0.,6.);
 158   fOutput->Add(fdEdx);
 159
 160   fdEdxPid=new TH2F("fdEdxPid","dE/dx with PID",50,0.2,3,50,0.,6.);
 161   fOutput->Add(fdEdxPid);
 162
 163   fK0sM =
 164   new TH2F("fK0sM", "Mass for K^{0}_{s}", 50,0.448,0.548,nBins,xBins);
 165   fK0sM->GetXaxis()->SetTitle("Mass (GeV/c)");
 166   fOutput->Add(fK0sM);
 167
 168   fK0sSi =
 169   new TH2F("fK0sSi","L_{T} vs p_{T} for K^{0}_{s}, side-band subtracted",
 170   nBins,xBins,lbins,0.,ltMax);
 171   fK0sSi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 172   fK0sSi->GetYaxis()->SetTitle("L_{T} (cm)");
 173   fOutput->Add(fK0sSi);
 174
 175   fK0sMC =
 176   new TH2F("fK0sMC","L_{T} vs p_{T} for K^{0}_{s}, from MC stack",
 177   nBins,xBins,lbins,0.,ltMax);
 178   fK0sMC->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 179   fK0sMC->GetYaxis()->SetTitle("L_{T} (cm)");
 180   fOutput->Add(fK0sMC);
 181
 182   fK0sAs =
 183   new TH2F("fK0sAs", "L_{T} vs p_{T} for K^{0}_{s}, associated",
 184   nBins,xBins,lbins,0.,ltMax);
 185   fK0sAs->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 186   fK0sAs->GetYaxis()->SetTitle("L_{T} (cm)");
 187   fOutput->Add(fK0sAs);
 188
 189   //----------------------
 190
 191   fLambdaM =
 192   new TH2F("fLambdaM","Mass for \\Lambda", 100, 1.065, 1.165,nBins,xBins);
 193   fLambdaM->GetXaxis()->SetTitle("Mass (GeV/c)");
 194   fOutput->Add(fLambdaM);
 195
 196  fLambdaSi =
 197   new TH2F("fLambdaSi","L_{T} vs p_{T} for \\Lambda, side-band subtructed",
 198   nBins,xBins,lbins,0.,ltMax);
 199   fLambdaSi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 200   fLambdaSi->GetYaxis()->SetTitle("L_{T} (cm)");
 201   fOutput->Add(fLambdaSi);
 202
 203   fLambdaMC =
 204   new TH2F("fLambdaMC","L_{T} vs p_{T} for \\Lambda, from MC stack",
 205   nBins,xBins,lbins,0.,ltMax);
 206   fLambdaMC->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 207   fLambdaMC->GetYaxis()->SetTitle("L_{T} (cm)");
 208   fOutput->Add(fLambdaMC);
 209
 210   fLambdaAs =
 211   new TH2F("fLambdaAs","L_{T} vs p_{T} for \\Lambda, associated",
 212   nBins,xBins,lbins,0.,ltMax);
 213   fLambdaAs->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 214   fLambdaAs->GetYaxis()->SetTitle("L_{T} (cm)");
 215   fOutput->Add(fLambdaAs);
 216
 217   //----------------------
 218
 219   fLambdaEff=fLambdaAs->ProjectionX();
 220   fLambdaEff->SetName("fLambdaEff");
 221   fLambdaEff->SetTitle("Efficiency for #Lambda");
 222   fOutput->Add(fLambdaEff);
 223
 224   fLambdaPt=fLambdaAs->ProjectionX();
 225   fLambdaPt->SetName("fLambdaPt");
 226   fLambdaPt->SetTitle("Raw #Lambda pT spectrum");
 227   fOutput->Add(fLambdaPt);
 228   //----------------------
 229
 230   fLambdaBarM =
 231   new TH2F("fLambdaBarM","Mass for anti-\\Lambda",100,1.065,1.165,nBins,xBins);
 232   fLambdaBarM->GetXaxis()->SetTitle("Mass (GeV/c)");
 233   fOutput->Add(fLambdaBarM);
 234
 235   fLambdaBarSi =
 236   new TH2F("fLambdaBarSi","L_{T} vs p_{T} for anti-\\Lambda, side-band subtructed",
 237   nBins,xBins,lbins,0.,ltMax);
 238   fLambdaBarSi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 239   fLambdaBarSi->GetYaxis()->SetTitle("L_{T} (cm)");
 240   fOutput->Add(fLambdaBarSi);
 241
 242   fLambdaBarMC =
 243   new TH2F("fLambdaBarMC","L_{T} vs p_{T} for anti-\\Lambda, from MC stack",
 244   nBins,xBins,lbins,0.,ltMax);
 245   fLambdaBarMC->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 246   fLambdaBarMC->GetYaxis()->SetTitle("L_{T} (cm)");
 247   fOutput->Add(fLambdaBarMC);
 248
 249   fLambdaBarAs =
 250   new TH2F("fLambdaBarAs","L_{T} vs p_{T} for anti-\\Lambda, associated",
 251   nBins,xBins,lbins,0.,ltMax);
 252   fLambdaBarAs->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 253   fLambdaBarAs->GetYaxis()->SetTitle("L_{T} (cm)");
 254   fOutput->Add(fLambdaBarAs);
 255
 256
 257   //----------------------
 258
 259   fLambdaBarEff=fLambdaBarAs->ProjectionX();
 260   fLambdaBarEff->SetName("fLambdaBarEff");
 261   fLambdaBarEff->SetTitle("Efficiency for anti-#Lambda");
 262   fOutput->Add(fLambdaBarEff);
 263
 264   fLambdaBarPt=fLambdaBarAs->ProjectionX();
 265   fLambdaBarPt->SetName("fLambdaBarPt");
 266   fLambdaBarPt->SetTitle("Raw anti-#Lambda pT spectrum");
 267   fOutput->Add(fLambdaBarPt);
 268
 269   //----------------------
 270   fLambdaFromXi=new TH3F("fLambdaFromXi","L_{T} vs p_{T} vs p_{T} of \\Xi for \\Lambda from \\Xi",
 271   nBins,xBins,lbins,yBins,kbins,zBins);
 272   fOutput->Add(fLambdaFromXi);
 273
 274   fXiM  =
 275   new TH2F("fXiM", "\\Xi mass distribution", 50, 1.271, 1.371,33,0.,ptMax+2);
 276   fOutput->Add(fXiM);
 277
 278   fXiSiP  = new TH1F("fXiSiP", "Pt for \\Xi, side-band subracted",
 279   33,0.,ptMax+2);
 280   fOutput->Add(fXiSiP);
 281
 282
 283   fLambdaBarFromXiBar=new TH3F("fLambdaBarFromXiBar","L_{T} vs p_{T} vs p_{T} of anti-\\Xi for anti-\\Lambda from anti-\\Xi",
 284   nBins,xBins,lbins,yBins,kbins,zBins);
 285   fOutput->Add(fLambdaBarFromXiBar);
 286
 287   fXiBarM  =
 288   new TH2F("fXiBarM", "anti-\\Xi mass distribution", 50, 1.271, 1.371,33,0.,ptMax+2);
 289   fOutput->Add(fXiBarM);
 290
 291   fXiBarSiP  = new TH1F("fXiBarSiP", "Pt for anti-\\Xi, side-band subracted",
 292   33,0.,ptMax+2);
 293   fOutput->Add(fXiBarSiP);
 294
 295 }
 296
 297 Bool_t AliV0CutVariations::AcceptTracks() {
 298   //if (!t->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
 299   //if (t->GetKinkIndex(0)>0) return kFALSE;
 300
 301   /*
 302   if (fTreeVariableLeastNbrCrossedRows < fTPCcr) return kFALSE;
 303   if (fTreeVariableLeastRatioCrossedRowsOverFindable < fTPCcrfd) return kFALSE;
 304   */
 305   if (TMath::Abs(fTreeVariableNegEta) > 0.8) return kFALSE;
 306   if (TMath::Abs(fTreeVariablePosEta) > 0.8) return kFALSE;
 307
 308   fTPCrows->Fill(fTreeVariableLeastNbrCrossedRows);
 309   fTPCratio->Fill(fTreeVariableLeastRatioCrossedRowsOverFindable);
 310
 311   return kTRUE;
 312 }
 313
 314 Bool_t AliV0CutVariations::AcceptV0()
 315 {
 316   //if (v0->GetOnFlyStatus()) return kFALSE;
 317
 318   if (fTreeVariableV0CosineOfPointingAngle < fCPA) return kFALSE;
 319   if (fTreeVariableDcaV0Daughters > fDCA) return kFALSE;
 320
 321   if (!AcceptTracks()) return kFALSE;
 322
 323   if (TMath::Abs(fTreeVariableDcaNegToPrimVertex) < fDCApv) return kFALSE;
 324   if (TMath::Abs(fTreeVariableDcaPosToPrimVertex) < fDCApv) return kFALSE;
 325
 326   if (fTreeVariableV0Radius < 0.9) return kFALSE;
 327   if (fTreeVariableV0Radius > 100) return kFALSE;
 328
 329   fCosPA->Fill(fTreeVariableV0CosineOfPointingAngle);
 330   fDtrDCA->Fill(fTreeVariableDcaV0Daughters);
 331   fPrimDCA->Fill(fTreeVariableDcaNegToPrimVertex);
 332   fPrimDCA->Fill(fTreeVariableDcaPosToPrimVertex);
 333
 334   return kTRUE;
 335 }
 336
 337 Bool_t AliV0CutVariations::AcceptPID(Int_t code)
 338 {
 339   //const AliAODPid *pid=ptrack->GetDetPid();
 340   //if (!pid) return kTRUE;
 341   if (code > 0) {
 342      if (fTreeVariablePosTransvMomentum > 1.) return kTRUE;
 343   } else {
 344      if (fTreeVariableNegTransvMomentum > 1.) return kTRUE;
 345   }
 346
 347
 348   if (fIsMC) {
 349     // MC PID
 350
 351     if (fSelectNonInjected && (!fTreeVariableIsNonInjected)) return kFALSE;
 352
 353     if (code > 0) {
 354       if (fTreeVariablePIDPositive == code) return kTRUE;
 355     } else {
 356       if (fTreeVariablePIDNegative == code) return kTRUE;
 357     }
 358   } else {
 359     // Real PID
 360     if (code > 0) {
 361        if (TMath::Abs(fTreeVariableNSigmasPosProton) < 3.) return kTRUE;
 362     } else {
 363        if (TMath::Abs(fTreeVariableNSigmasNegProton) < 3.) return kTRUE;
 364     }
 365   }
 366
 367   return kFALSE;
 368 }
 369
 370
 371 Bool_t AliV0CutVariations::Process(Long64_t entry)
 372 {
 373    // The Process() function is called for each entry in the tree (or possibly
 374    // keyed object in the case of PROOF) to be processed. The entry argument
 375    // specifies which entry in the currently loaded tree is to be processed.
 376    // It can be passed to either AliV0CutVariations::GetEntry() or TBranch::GetEntry()
 377    // to read either all or the required parts of the data. When processing
 378    // keyed objects with PROOF, the object is already loaded and is available
 379    // via the fObject pointer.
 380    //
 381    // This function should contain the "body" of the analysis. It can contain
 382    // simple or elaborate selection criteria, run algorithms on the data
 383    // of the event and typically fill histograms.
 384    //
 385    // The processing can be stopped by calling Abort().
 386    //
 387    // Use fStatus to set the return value of TTree::Process().
 388    //
 389    // The return value is currently not used.
 390
 391    const Double_t yMax=0.5;
 392    const Double_t pMin=0.0;
 393    //const Double_t lMax=0.001;
 394
 395    fChain->GetTree()->GetEntry(entry);
 396
 397    //cout<<entry<<'\r';
 398
 399    if (fTreeVariableMultiplicity<fCMin) return kFALSE;
 400    if (fTreeVariableMultiplicity>fCMax) return kFALSE;
 401
 402    fMult->Fill(fTreeVariableMultiplicity);
 403
 404    Double_t pt=0;
 405    Double_t lt=0;
 406    //+++++++ MC
 407    if (fIsMC) {
 408       if (fSelectNonInjected && (!fTreeVariableIsNonInjected)) goto real;
 409
 410       Int_t code=fTreeVariablePID;
 411
 412       if (code != kK0Short)
 413          if (code != kLambda0)
 414             if (code != kLambda0Bar) goto real;
 415
 416        pt=fTreeVariablePtMC;
 417        if (pt < pMin) goto real;
 418
 419        if (TMath::Abs(fTreeVariableRapMC) > yMax) goto real;
 420
 421        //if (TMath::Abs(fTreeVariableV0CreationRadius) > lMax) goto real;
 422        if (fTreeVariablePrimaryStatus != 1) goto real;
 423
 424        lt=fTreeVariableV0Radius;
 425        switch (code) {
 426        case kK0Short:
 427           fK0sMC->Fill(pt,lt);
 428           break;
 429        case kLambda0:
 430           fLambdaMC->Fill(pt,lt);
 431           break;
 432        case kLambda0Bar:
 433           fLambdaBarMC->Fill(pt,lt);
 434           break;
 435        default: break;
 436        }
 437
 438   }
 439    //+++++++
 440
 441  real:
 442
 443    pt=fTreeVariablePt;
 444    if (pt < pMin) return kFALSE;
 445
 446    if (!AcceptV0()) return kFALSE;
 447
 448    lt=fTreeVariableV0Radius;
 449    if (lt/pt > 3*7.89/1.1157) return kFALSE;
 450
 451    //--- V0 switches
 452    Bool_t isK0s=kTRUE;
 453    Bool_t isLambda=kTRUE;
 454    Bool_t isLambdaBar=kTRUE;
 455
 456    if (0.4977*lt/pt > 3*2.68) isK0s=kFALSE;
 457    if (1.1157*lt/pt > 3*7.89) isLambdaBar=isLambda=kFALSE;
 458
 459    if (fTreeVariablePtArmV0<0.2*TMath::Abs(fTreeVariableAlphaV0)) isK0s=kFALSE;
 460
 461    if (!AcceptPID(kProton)) isLambda=kFALSE;
 462    if (!AcceptPID(kProtonBar)) isLambdaBar=kFALSE;
 463
 464    Double_t yK0s=TMath::Abs(fTreeVariableRapK0Short);
 465    Double_t yLam=TMath::Abs(fTreeVariableRapLambda);
 466    if (yK0s > yMax) isK0s=kFALSE;
 467    if (yLam > yMax) isLambda=isLambdaBar=kFALSE;
 468    //---
 469
 470    Double_t mass=0., m=0., s=0.;
 471    if (isK0s) {
 472       mass=fTreeVariableInvMassK0s;
 473       fK0sM->Fill(mass,pt);
 474
 475       m=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
 476       //s=0.0044 + (0.008-0.0044)/(10-1)*(pt - 1.);
 477       s=0.0029 + 0.00077*pt;
 478       if (TMath::Abs(m-mass) < 3*s) {
 479          fK0sSi->Fill(pt,lt);
 480       } else {
 481          isK0s=kFALSE;
 482       }
 483       if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 484          fK0sSi->Fill(pt,lt,-1);
 485       }
 486       if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 487          fK0sSi->Fill(pt,lt,-1);
 488       }
 489    }
 490
 491    if (isLambda) {
 492       mass=fTreeVariableInvMassLambda;
 493       fLambdaM->Fill(mass,pt);
 494
 495       m=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
 496       //s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
 497       s=0.0021 + 0.00022*pt;
 498       if (TMath::Abs(m-mass) < 3*s) {
 499          fLambdaSi->Fill(pt,lt);
 500       } else {
 501          isLambda=kFALSE;
 502       }
 503       if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 504          fLambdaSi->Fill(pt,lt,-1);
 505       }
 506       if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 507          fLambdaSi->Fill(pt,lt,-1);
 508       }
 509    }
 510
 511    if (isLambdaBar) {
 512       mass=fTreeVariableInvMassAntiLambda;
 513       fLambdaBarM->Fill(mass,pt);
 514
 515       m=TDatabasePDG::Instance()->GetParticle(kLambda0Bar)->Mass();
 516       //s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
 517       s=0.0021 + 0.00022*pt;
 518       if (TMath::Abs(m-mass) < 3*s) {
 519          fLambdaBarSi->Fill(pt,lt);
 520       } else {
 521          isLambdaBar=kFALSE;
 522       }
 523       if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 524          fLambdaBarSi->Fill(pt,lt,-1);
 525       }
 526       if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 527          fLambdaBarSi->Fill(pt,lt,-1);
 528       }
 529    }
 530
 531    if (!fIsMC) return kFALSE;
 532    if (fSelectNonInjected && (!fTreeVariableIsNonInjected)) return kFALSE;
 533
 534
 535    //++++++ MC
 536    if (!isK0s)
 537       if (!isLambda)
 538          if (!isLambdaBar) return kFALSE;//check MC only for the accepted V0s
 539
 540    Int_t code=fTreeVariablePID;
 541    if (code != kK0Short)
 542       if (code != kLambda0)
 543          if (code != kLambda0Bar) return kFALSE;
 544
 545    Double_t ptAs=fTreeVariablePtMC;
 546    if (ptAs < pMin) return kFALSE;
 547
 548    Double_t yAs=fTreeVariableRapMC;
 549    if (TMath::Abs(yAs) > yMax) return kFALSE;
 550
 551    Double_t ltAs=fTreeVariableV0Radius;
 552       if (fTreeVariablePrimaryStatus == 1) {
 553          switch (code) {
 554          case kK0Short:
 555             if (isK0s)       fK0sAs->Fill(ptAs,ltAs);
 556             break;
 557          case kLambda0:
 558             if (isLambda)    fLambdaAs->Fill(ptAs,ltAs);
 559             break;
 560          case kLambda0Bar:
 561             if (isLambdaBar) fLambdaBarAs->Fill(ptAs,ltAs);
 562             break;
 563          default: break;
 564          }
 565       } else {
 566          if (code==kK0Short) return kTRUE;
 567
 568          Int_t xcode=fTreeVariablePIDMother;
 569          Double_t xpt=fTreeVariablePtXiMother;
 570
 571          switch (code) {
 572          case kLambda0:
 573             if (isLambda)
 574             if ((xcode==kXiMinus) || (xcode==3322))
 575                fLambdaFromXi->Fill(ptAs,ltAs,xpt);
 576             break;
 577          case kLambda0Bar:
 578             if (isLambdaBar)
 579             if ((xcode==kXiPlusBar)||(xcode==-3322))
 580                fLambdaBarFromXiBar->Fill(ptAs,ltAs,xpt);
 581             break;
 582          default: break;
 583          }
 584    }
 585
 586    return kTRUE;
 587 }
 588
 589 void AliV0CutVariations::SlaveTerminate()
 590 {
 591    // The SlaveTerminate() function is called after all entries or objects
 592    // have been processed. When running with PROOF SlaveTerminate() is called
 593    // on each slave server.
 594
 595 }
 596
 597 void AliV0CutVariations::Terminate()
 598 {
 599    // The Terminate() function is the last function to be called during
 600    // a query. It always runs on the client, it can be used to present
 601    // the results graphically or save the results to file.
 602
 603   cout<<endl;
 604
 605   if (!fOutput) {
 606      Printf("ERROR: fOutput not available");
 607      return;
 608   }
 609
 610   fMult = dynamic_cast<TH1F*>(fOutput->FindObject("fMult")) ;
 611   if (!fMult) {
 612      Printf("ERROR: fMult not available");
 613      return;
 614   }
 615
 616   fdEdx = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdx")) ;
 617   if (!fdEdx) {
 618      Printf("ERROR: fdEdx not available");
 619      return;
 620   }
 621
 622   fdEdxPid = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdxPid")) ;
 623   if (!fdEdxPid) {
 624      Printf("ERROR: fdEdxPid not available");
 625      return;
 626   }
 627
 628
 629   fK0sMC = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sMC")) ;
 630   if (!fK0sMC) {
 631      Printf("ERROR: fK0sMC not available");
 632      return;
 633   }
 634   TH1D *k0sMcPx=fK0sMC->ProjectionX(); k0sMcPx->Sumw2();
 635   fK0sAs = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sAs")) ;
 636   if (!fK0sAs) {
 637      Printf("ERROR: fK0sAs not available");
 638      return;
 639   }
 640   TH1D *k0sAsPx=fK0sAs->ProjectionX();
 641   k0sAsPx->Sumw2(); //k0sAsPx->Scale(0.69);
 642
 643
 644
 645   // Lambda histograms
 646   fLambdaFromXi = dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaFromXi")) ;
 647   if (!fLambdaFromXi) {
 648      Printf("ERROR: fLambdaFromXi not available");
 649      return;
 650   }
 651   TH1D *lambdaFromXiPx=fLambdaFromXi->ProjectionX(); lambdaFromXiPx->Sumw2();
 652
 653   fLambdaMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaMC")) ;
 654   if (!fLambdaMC) {
 655      Printf("ERROR: fLambdaMC not available");
 656      return;
 657   }
 658   TH1D *lambdaMcPx=fLambdaMC->ProjectionX(); lambdaMcPx->Sumw2();
 659
 660   fLambdaAs = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaAs")) ;
 661   if (!fLambdaAs) {
 662      Printf("ERROR: fLambdaAs not available");
 663      return;
 664   }
 665   TH1D *lambdaAsPx=fLambdaAs->ProjectionX();
 666   lambdaAsPx->Sumw2(); //lambdaAsPx->Scale(0.64);
 667
 668   fLambdaSi = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaSi")) ;
 669   if (!fLambdaSi) {
 670      Printf("ERROR: fLambdaSi not available");
 671      return;
 672   }
 673   TH1D *lambdaSiPx=fLambdaSi->ProjectionX();
 674   lambdaSiPx->SetName("fLambdaPt");
 675   lambdaSiPx->Sumw2();
 676
 677   fLambdaEff = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaEff")) ;
 678   if (!fLambdaEff) {
 679      Printf("ERROR: fLambdaEff not available");
 680      return;
 681   }
 682   fLambdaPt = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaPt")) ;
 683   if (!fLambdaPt) {
 684      Printf("ERROR: fLambdaPt not available");
 685      return;
 686   }
 687
 688
 689   // anti-Lambda histograms
 690   fLambdaBarFromXiBar =
 691   dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaBarFromXiBar")) ;
 692   if (!fLambdaBarFromXiBar) {
 693      Printf("ERROR: fLambdaBarFromXiBar not available");
 694      return;
 695   }
 696   TH1D *lambdaBarFromXiBarPx=
 697   fLambdaBarFromXiBar->ProjectionX("Bar"); lambdaBarFromXiBarPx->Sumw2();
 698
 699   fLambdaBarMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarMC")) ;
 700   if (!fLambdaBarMC) {
 701      Printf("ERROR: fLambdaBarMC not available");
 702      return;
 703   }
 704   TH1D *lambdaBarMcPx=fLambdaBarMC->ProjectionX(); lambdaBarMcPx->Sumw2();
 705
 706   fLambdaBarAs = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarAs")) ;
 707   if (!fLambdaBarAs) {
 708      Printf("ERROR: fLambdaBarAs not available");
 709      return;
 710   }
 711   TH1D *lambdaBarAsPx=fLambdaBarAs->ProjectionX();
 712   lambdaBarAsPx->Sumw2(); //lambdaBarAsPx->Scale(0.64);
 713
 714   fLambdaBarSi = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarSi")) ;
 715   if (!fLambdaBarSi) {
 716      Printf("ERROR: fLambdaBarSi not available");
 717      return;
 718   }
 719   TH1D *lambdaBarSiPx=fLambdaBarSi->ProjectionX();
 720   lambdaBarSiPx->SetName("fLambdaBarPt");
 721   lambdaBarSiPx->Sumw2();
 722
 723   fLambdaBarEff = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaBarEff")) ;
 724   if (!fLambdaBarEff) {
 725      Printf("ERROR: fLambdaBarEff not available");
 726      return;
 727   }
 728   fLambdaBarPt = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaBarPt")) ;
 729   if (!fLambdaBarPt) {
 730      Printf("ERROR: fLambdaBarPt not available");
 731      return;
 732   }
 733
 734
 735   if (!gROOT->IsBatch()) {
 736
 737     TCanvas *c1 = new TCanvas("c1","Mulitplicity");
 738     c1->SetLogy();
 739     fMult->DrawCopy() ;
 740
 741     new TCanvas("c2","dE/dx");
 742     fdEdx->DrawCopy() ;
 743
 744     new TCanvas("c3","dE/dx with PID");
 745     fdEdxPid->DrawCopy() ;
 746
 747     if (fIsMC) {
 748        /*
 749        TH1D effK(*k0sAsPx); effK.SetTitle("Efficiency for K0s");
 750        effK.Divide(k0sAsPx,k0sMcPx,1,1,"b");
 751        new TCanvas("c4","Efficiency for K0s");
 752        effK.DrawCopy("E") ;
 753        */
 754
 755        //+++ Lambdas
 756        fLambdaEff->Divide(lambdaAsPx,lambdaMcPx,1,1,"b");
 757        new TCanvas("c5","Efficiency for #Lambda");
 758        fLambdaEff->DrawCopy("E") ;
 759
 760        lambdaSiPx->Add(lambdaFromXiPx,-1);
 761        lambdaSiPx->Divide(fLambdaEff);
 762
 763        new TCanvas("c6","Corrected #Lambda pt");
 764        lambdaSiPx->SetTitle("Corrected #Lambda pt");
 765       *fLambdaPt = *lambdaSiPx;
 766        fLambdaPt->SetLineColor(2);
 767        fLambdaPt->DrawCopy("E");
 768
 769        lambdaMcPx->DrawCopy("same");
 770
 771
 772        //+++ anti-Lambdas
 773        fLambdaBarEff->Divide(lambdaBarAsPx,lambdaBarMcPx,1,1,"b");
 774        new TCanvas("c7","Efficiency for anti-#Lambda");
 775        fLambdaBarEff->DrawCopy("E") ;
 776
 777        lambdaBarSiPx->Add(lambdaBarFromXiBarPx,-1);
 778        lambdaBarSiPx->Divide(fLambdaBarEff);
 779
 780        new TCanvas("c8","Corrected anti-#Lambda pt");
 781        lambdaBarSiPx->SetTitle("Corrected anti-#Lambda pt");
 782       *fLambdaBarPt = *lambdaBarSiPx;
 783        fLambdaBarPt->SetLineColor(2);
 784        fLambdaBarPt->DrawCopy("E");
 785
 786        lambdaBarMcPx->DrawCopy("same");
 787     } else {
 788        new TCanvas("c6","Raw #Lambda pt");
 789        lambdaSiPx->SetTitle("Raw #Lambda pt");
 790       *fLambdaPt = *lambdaSiPx;
 791        fLambdaPt->SetLineColor(2);
 792        fLambdaPt->DrawCopy("E");
 793
 794
 795        new TCanvas("c7","Raw anti-#Lambda pt");
 796        lambdaBarSiPx->SetTitle("Raw anti-#Lambda pt");
 797       *fLambdaBarPt = *lambdaBarSiPx;
 798        fLambdaBarPt->SetLineColor(2);
 799        fLambdaBarPt->DrawCopy("E");
 800     }
 801   }
 802
 803   TFile *file = 0;
 804   if (fIsMC)
 805     if (fSelectNonInjected)
 806        file=TFile::Open("AliV0CutVariationsMC_nonInj.root", "RECREATE");
 807     else
 808        file=TFile::Open("AliV0CutVariationsMC.root", "RECREATE");
 809   else
 810      file=TFile::Open("AliV0CutVariations.root", "RECREATE");
 811
 812   fOutput->Write() ;
 813   file->Close() ;
 814   delete file ;
 815
 816 }