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   if (fTreeVariableLeastNbrCrossedRows < fTPCcr) return kFALSE;
 302   if (fTreeVariableLeastRatioCrossedRowsOverFindable < fTPCcrfd) return kFALSE;
 303
 304   if (TMath::Abs(fTreeVariableNegEta) > 0.8) return kFALSE;
 305   if (TMath::Abs(fTreeVariablePosEta) > 0.8) return kFALSE;
 306
 307   fTPCrows->Fill(fTreeVariableLeastNbrCrossedRows);
 308   fTPCratio->Fill(fTreeVariableLeastRatioCrossedRowsOverFindable);
 309
 310   return kTRUE;
 311 }
 312
 313 Bool_t AliV0CutVariations::AcceptV0()
 314 {
 315   //if (v0->GetOnFlyStatus()) return kFALSE;
 316
 317   if (fTreeVariableV0CosineOfPointingAngle < fCPA) return kFALSE;
 318   if (fTreeVariableDcaV0Daughters > fDCA) return kFALSE;
 319
 320   if (!AcceptTracks()) return kFALSE;
 321
 322   if (TMath::Abs(fTreeVariableDcaNegToPrimVertex) < fDCApv) return kFALSE;
 323   if (TMath::Abs(fTreeVariableDcaPosToPrimVertex) < fDCApv) return kFALSE;
 324
 325   if (fTreeVariableV0Radius < 0.9) return kFALSE;
 326   if (fTreeVariableV0Radius > 100) return kFALSE;
 327
 328   fCosPA->Fill(fTreeVariableV0CosineOfPointingAngle);
 329   fDtrDCA->Fill(fTreeVariableDcaV0Daughters);
 330   fPrimDCA->Fill(fTreeVariableDcaNegToPrimVertex);
 331   fPrimDCA->Fill(fTreeVariableDcaPosToPrimVertex);
 332
 333   return kTRUE;
 334 }
 335
 336 Bool_t AliV0CutVariations::AcceptPID(Int_t code)
 337 {
 338   //const AliAODPid *pid=ptrack->GetDetPid();
 339   //if (!pid) return kTRUE;
 340   if (code > 0) {
 341      if (fTreeVariablePosTransvMomentum > 1.) return kTRUE;
 342   } else {
 343      if (fTreeVariableNegTransvMomentum > 1.) return kTRUE;
 344   }
 345
 346
 347   if (fIsMC) {
 348     // MC PID
 349
 350     if (fSelectNonInjected && (!fTreeVariableIsNonInjected)) return kFALSE;
 351
 352     if (code > 0) {
 353       if (fTreeVariablePIDPositive == code) return kTRUE;
 354     } else {
 355       if (fTreeVariablePIDNegative == code) return kTRUE;
 356     }
 357   } else {
 358     // Real PID
 359     if (code > 0) {
 360        if (TMath::Abs(fTreeVariableNSigmasPosProton) < 3.) return kTRUE;
 361     } else {
 362        if (TMath::Abs(fTreeVariableNSigmasNegProton) < 3.) return kTRUE;
 363     }
 364   }
 365
 366   return kFALSE;
 367 }
 368
 369
 370 Bool_t AliV0CutVariations::Process(Long64_t entry)
 371 {
 372    // The Process() function is called for each entry in the tree (or possibly
 373    // keyed object in the case of PROOF) to be processed. The entry argument
 374    // specifies which entry in the currently loaded tree is to be processed.
 375    // It can be passed to either AliV0CutVariations::GetEntry() or TBranch::GetEntry()
 376    // to read either all or the required parts of the data. When processing
 377    // keyed objects with PROOF, the object is already loaded and is available
 378    // via the fObject pointer.
 379    //
 380    // This function should contain the "body" of the analysis. It can contain
 381    // simple or elaborate selection criteria, run algorithms on the data
 382    // of the event and typically fill histograms.
 383    //
 384    // The processing can be stopped by calling Abort().
 385    //
 386    // Use fStatus to set the return value of TTree::Process().
 387    //
 388    // The return value is currently not used.
 389
 390    const Double_t yMax=0.5;
 391    const Double_t pMin=0.0;
 392    //const Double_t lMax=0.001;
 393
 394    fChain->GetTree()->GetEntry(entry);
 395
 396    //cout<<entry<<'\r';
 397
 398    if (fTreeVariableMultiplicity<fCMin) return kFALSE;
 399    if (fTreeVariableMultiplicity>fCMax) return kFALSE;
 400
 401    fMult->Fill(fTreeVariableMultiplicity);
 402
 403    Double_t pt=0;
 404    Double_t lt=0;
 405    //+++++++ MC
 406    if (fIsMC) {
 407       if (fSelectNonInjected && (!fTreeVariableIsNonInjected)) goto real;
 408
 409       Int_t code=fTreeVariablePID;
 410
 411       if (code != kK0Short)
 412          if (code != kLambda0)
 413             if (code != kLambda0Bar) goto real;
 414
 415        pt=fTreeVariablePtMC;
 416        if (pt < pMin) goto real;
 417
 418        if (TMath::Abs(fTreeVariableRapMC) > yMax) goto real;
 419
 420        //if (TMath::Abs(fTreeVariableV0CreationRadius) > lMax) goto real;
 421        if (fTreeVariablePrimaryStatus != 1) goto real;
 422
 423        lt=fTreeVariableV0Radius;
 424        switch (code) {
 425        case kK0Short:
 426           fK0sMC->Fill(pt,lt);
 427           break;
 428        case kLambda0:
 429           fLambdaMC->Fill(pt,lt);
 430           break;
 431        case kLambda0Bar:
 432           fLambdaBarMC->Fill(pt,lt);
 433           break;
 434        default: break;
 435        }
 436
 437   }
 438    //+++++++
 439
 440  real:
 441
 442    pt=fTreeVariablePt;
 443    if (pt < pMin) return kFALSE;
 444
 445    if (!AcceptV0()) return kFALSE;
 446
 447    lt=fTreeVariableV0Radius;
 448    if (lt/pt > 3*7.89/1.1157) return kFALSE;
 449
 450    //--- V0 switches
 451    Bool_t isK0s=kTRUE;
 452    Bool_t isLambda=kTRUE;
 453    Bool_t isLambdaBar=kTRUE;
 454
 455    if (0.4977*lt/pt > 3*2.68) isK0s=kFALSE;
 456    if (1.1157*lt/pt > 3*7.89) isLambdaBar=isLambda=kFALSE;
 457
 458    if (!AcceptPID(kProton)) isLambda=kFALSE;
 459    if (!AcceptPID(kProtonBar)) isLambdaBar=kFALSE;
 460
 461    Double_t yK0s=fTreeVariableRapK0Short;
 462    Double_t yLam=fTreeVariableRapLambda;
 463    if (yK0s > yMax) isK0s=kFALSE;
 464    if (yLam > yMax) isLambda=isLambdaBar=kFALSE;
 465    //---
 466
 467    Double_t mass=0., m=0., s=0.;
 468    if (isK0s) {
 469       mass=fTreeVariableInvMassK0s;
 470       fK0sM->Fill(mass,pt);
 471
 472       m=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
 473       s=0.0044 + (0.008-0.0044)/(10-1)*(pt - 1.);
 474       if (TMath::Abs(m-mass) < 3*s) {
 475          fK0sSi->Fill(pt,lt);
 476       } else {
 477          isK0s=kFALSE;
 478       }
 479       if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 480          fK0sSi->Fill(pt,lt,-1);
 481       }
 482       if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 483          fK0sSi->Fill(pt,lt,-1);
 484       }
 485    }
 486
 487    if (isLambda) {
 488       mass=fTreeVariableInvMassLambda;
 489       fLambdaM->Fill(mass,pt);
 490
 491       m=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
 492       //s=0.0027 + (0.004-0.0027)/(10-1)*(pt-1);
 493       //s=0.0015 + (0.002-0.0015)/(2.6-1)*(pt-1);
 494       s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
 495       if (TMath::Abs(m-mass) < 3*s) {
 496          fLambdaSi->Fill(pt,lt);
 497       } else {
 498          isLambda=kFALSE;
 499       }
 500       if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 501          fLambdaSi->Fill(pt,lt,-1);
 502       }
 503       if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 504          fLambdaSi->Fill(pt,lt,-1);
 505       }
 506    }
 507
 508    if (isLambdaBar) {
 509       mass=fTreeVariableInvMassAntiLambda;
 510       fLambdaBarM->Fill(mass,pt);
 511
 512       m=TDatabasePDG::Instance()->GetParticle(kLambda0Bar)->Mass();
 513       //s=0.0027 + (0.004-0.0027)/(10-1)*(pt-1);
 514       //s=0.0015 + (0.002-0.0015)/(2.6-1)*(pt-1);
 515       s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
 516       if (TMath::Abs(m-mass) < 3*s) {
 517          fLambdaBarSi->Fill(pt,lt);
 518       } else {
 519          isLambdaBar=kFALSE;
 520       }
 521       if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 522          fLambdaBarSi->Fill(pt,lt,-1);
 523       }
 524       if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 525          fLambdaBarSi->Fill(pt,lt,-1);
 526       }
 527    }
 528
 529    if (!fIsMC) return kFALSE;
 530    if (fSelectNonInjected && (!fTreeVariableIsNonInjected)) return kFALSE;
 531
 532
 533    //++++++ MC
 534    if (!isK0s)
 535       if (!isLambda)
 536          if (!isLambdaBar) return kFALSE;//check MC only for the accepted V0s
 537
 538    Int_t code=fTreeVariablePID;
 539    if (code != kK0Short)
 540       if (code != kLambda0)
 541          if (code != kLambda0Bar) return kFALSE;
 542
 543    Double_t ptAs=fTreeVariablePtMC;
 544    if (ptAs < pMin) return kFALSE;
 545
 546    Double_t yAs=fTreeVariableRapMC;
 547    if (TMath::Abs(yAs) > yMax) return kFALSE;
 548
 549    Double_t ltAs=fTreeVariableV0Radius;
 550       if (fTreeVariablePrimaryStatus == 1) {
 551          switch (code) {
 552          case kK0Short:
 553             if (isK0s)       fK0sAs->Fill(ptAs,ltAs);
 554             break;
 555          case kLambda0:
 556             if (isLambda)    fLambdaAs->Fill(ptAs,ltAs);
 557             break;
 558          case kLambda0Bar:
 559             if (isLambdaBar) fLambdaBarAs->Fill(ptAs,ltAs);
 560             break;
 561          default: break;
 562          }
 563       } else {
 564          if (code==kK0Short) return kTRUE;
 565
 566          Int_t xcode=fTreeVariablePIDMother;
 567          Double_t xpt=fTreeVariablePtXiMother;
 568
 569          switch (code) {
 570          case kLambda0:
 571             if (isLambda)
 572             if ((xcode==kXiMinus) || (xcode==3322))
 573                fLambdaFromXi->Fill(ptAs,ltAs,xpt);
 574             break;
 575          case kLambda0Bar:
 576             if (isLambdaBar)
 577             if ((xcode==kXiPlusBar)||(xcode==-3322))
 578                fLambdaBarFromXiBar->Fill(ptAs,ltAs,xpt);
 579             break;
 580          default: break;
 581          }
 582    }
 583
 584    return kTRUE;
 585 }
 586
 587 void AliV0CutVariations::SlaveTerminate()
 588 {
 589    // The SlaveTerminate() function is called after all entries or objects
 590    // have been processed. When running with PROOF SlaveTerminate() is called
 591    // on each slave server.
 592
 593 }
 594
 595 void AliV0CutVariations::Terminate()
 596 {
 597    // The Terminate() function is the last function to be called during
 598    // a query. It always runs on the client, it can be used to present
 599    // the results graphically or save the results to file.
 600
 601   cout<<endl;
 602
 603   if (!fOutput) {
 604      Printf("ERROR: fOutput not available");
 605      return;
 606   }
 607
 608   fMult = dynamic_cast<TH1F*>(fOutput->FindObject("fMult")) ;
 609   if (!fMult) {
 610      Printf("ERROR: fMult not available");
 611      return;
 612   }
 613
 614   fdEdx = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdx")) ;
 615   if (!fdEdx) {
 616      Printf("ERROR: fdEdx not available");
 617      return;
 618   }
 619
 620   fdEdxPid = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdxPid")) ;
 621   if (!fdEdxPid) {
 622      Printf("ERROR: fdEdxPid not available");
 623      return;
 624   }
 625
 626
 627   fK0sMC = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sMC")) ;
 628   if (!fK0sMC) {
 629      Printf("ERROR: fK0sMC not available");
 630      return;
 631   }
 632   TH1D *k0sMcPx=fK0sMC->ProjectionX(); k0sMcPx->Sumw2();
 633   fK0sAs = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sAs")) ;
 634   if (!fK0sAs) {
 635      Printf("ERROR: fK0sAs not available");
 636      return;
 637   }
 638   TH1D *k0sAsPx=fK0sAs->ProjectionX();
 639   k0sAsPx->Sumw2(); //k0sAsPx->Scale(0.69);
 640
 641
 642
 643   // Lambda histograms
 644   fLambdaFromXi = dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaFromXi")) ;
 645   if (!fLambdaFromXi) {
 646      Printf("ERROR: fLambdaFromXi not available");
 647      return;
 648   }
 649   TH1D *lambdaFromXiPx=fLambdaFromXi->ProjectionX(); lambdaFromXiPx->Sumw2();
 650
 651   fLambdaMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaMC")) ;
 652   if (!fLambdaMC) {
 653      Printf("ERROR: fLambdaMC not available");
 654      return;
 655   }
 656   TH1D *lambdaMcPx=fLambdaMC->ProjectionX(); lambdaMcPx->Sumw2();
 657
 658   fLambdaAs = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaAs")) ;
 659   if (!fLambdaAs) {
 660      Printf("ERROR: fLambdaAs not available");
 661      return;
 662   }
 663   TH1D *lambdaAsPx=fLambdaAs->ProjectionX();
 664   lambdaAsPx->Sumw2(); //lambdaAsPx->Scale(0.64);
 665
 666   fLambdaSi = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaSi")) ;
 667   if (!fLambdaSi) {
 668      Printf("ERROR: fLambdaSi not available");
 669      return;
 670   }
 671   TH1D *lambdaSiPx=fLambdaSi->ProjectionX();
 672   lambdaSiPx->SetName("fLambdaPt");
 673   lambdaSiPx->Sumw2();
 674
 675   fLambdaEff = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaEff")) ;
 676   if (!fLambdaEff) {
 677      Printf("ERROR: fLambdaEff not available");
 678      return;
 679   }
 680   fLambdaPt = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaPt")) ;
 681   if (!fLambdaPt) {
 682      Printf("ERROR: fLambdaPt not available");
 683      return;
 684   }
 685
 686
 687   // anti-Lambda histograms
 688   fLambdaBarFromXiBar =
 689   dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaBarFromXiBar")) ;
 690   if (!fLambdaBarFromXiBar) {
 691      Printf("ERROR: fLambdaBarFromXiBar not available");
 692      return;
 693   }
 694   TH1D *lambdaBarFromXiBarPx=
 695   fLambdaBarFromXiBar->ProjectionX("Bar"); lambdaBarFromXiBarPx->Sumw2();
 696
 697   fLambdaBarMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarMC")) ;
 698   if (!fLambdaBarMC) {
 699      Printf("ERROR: fLambdaBarMC not available");
 700      return;
 701   }
 702   TH1D *lambdaBarMcPx=fLambdaBarMC->ProjectionX(); lambdaBarMcPx->Sumw2();
 703
 704   fLambdaBarAs = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarAs")) ;
 705   if (!fLambdaBarAs) {
 706      Printf("ERROR: fLambdaBarAs not available");
 707      return;
 708   }
 709   TH1D *lambdaBarAsPx=fLambdaBarAs->ProjectionX();
 710   lambdaBarAsPx->Sumw2(); //lambdaBarAsPx->Scale(0.64);
 711
 712   fLambdaBarSi = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarSi")) ;
 713   if (!fLambdaBarSi) {
 714      Printf("ERROR: fLambdaBarSi not available");
 715      return;
 716   }
 717   TH1D *lambdaBarSiPx=fLambdaBarSi->ProjectionX();
 718   lambdaBarSiPx->SetName("fLambdaBarPt");
 719   lambdaBarSiPx->Sumw2();
 720
 721   fLambdaBarEff = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaBarEff")) ;
 722   if (!fLambdaBarEff) {
 723      Printf("ERROR: fLambdaBarEff not available");
 724      return;
 725   }
 726   fLambdaBarPt = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaBarPt")) ;
 727   if (!fLambdaBarPt) {
 728      Printf("ERROR: fLambdaBarPt not available");
 729      return;
 730   }
 731
 732
 733   if (!gROOT->IsBatch()) {
 734
 735     TCanvas *c1 = new TCanvas("c1","Mulitplicity");
 736     c1->SetLogy();
 737     fMult->DrawCopy() ;
 738
 739     new TCanvas("c2","dE/dx");
 740     fdEdx->DrawCopy() ;
 741
 742     new TCanvas("c3","dE/dx with PID");
 743     fdEdxPid->DrawCopy() ;
 744
 745     if (fIsMC) {
 746        /*
 747        TH1D effK(*k0sAsPx); effK.SetTitle("Efficiency for K0s");
 748        effK.Divide(k0sAsPx,k0sMcPx,1,1,"b");
 749        new TCanvas("c4","Efficiency for K0s");
 750        effK.DrawCopy("E") ;
 751        */
 752
 753        //+++ Lambdas
 754        fLambdaEff->Divide(lambdaAsPx,lambdaMcPx,1,1,"b");
 755        new TCanvas("c5","Efficiency for #Lambda");
 756        fLambdaEff->DrawCopy("E") ;
 757
 758        lambdaSiPx->Add(lambdaFromXiPx,-1);
 759        lambdaSiPx->Divide(fLambdaEff);
 760
 761        new TCanvas("c6","Corrected #Lambda pt");
 762        lambdaSiPx->SetTitle("Corrected #Lambda pt");
 763       *fLambdaPt = *lambdaSiPx;
 764        fLambdaPt->SetLineColor(2);
 765        fLambdaPt->DrawCopy("E");
 766
 767        lambdaMcPx->DrawCopy("same");
 768
 769
 770        //+++ anti-Lambdas
 771        fLambdaBarEff->Divide(lambdaBarAsPx,lambdaBarMcPx,1,1,"b");
 772        new TCanvas("c7","Efficiency for anti-#Lambda");
 773        fLambdaBarEff->DrawCopy("E") ;
 774
 775        lambdaBarSiPx->Add(lambdaBarFromXiBarPx,-1);
 776        lambdaBarSiPx->Divide(fLambdaBarEff);
 777
 778        new TCanvas("c8","Corrected anti-#Lambda pt");
 779        lambdaBarSiPx->SetTitle("Corrected anti-#Lambda pt");
 780       *fLambdaBarPt = *lambdaBarSiPx;
 781        fLambdaBarPt->SetLineColor(2);
 782        fLambdaBarPt->DrawCopy("E");
 783
 784        lambdaBarMcPx->DrawCopy("same");
 785     } else {
 786        new TCanvas("c6","Raw #Lambda pt");
 787        lambdaSiPx->SetTitle("Raw #Lambda pt");
 788       *fLambdaPt = *lambdaSiPx;
 789        fLambdaPt->SetLineColor(2);
 790        fLambdaPt->DrawCopy("E");
 791
 792
 793        new TCanvas("c7","Raw anti-#Lambda pt");
 794        lambdaBarSiPx->SetTitle("Raw anti-#Lambda pt");
 795       *fLambdaBarPt = *lambdaBarSiPx;
 796        fLambdaBarPt->SetLineColor(2);
 797        fLambdaBarPt->DrawCopy("E");
 798     }
 799   }
 800
 801   TFile *file = 0;
 802   if (fIsMC)
 803      file=TFile::Open("AliV0CutVariationsMC.root", "RECREATE");
 804   else
 805      file=TFile::Open("AliV0CutVariations.root", "RECREATE");
 806
 807   fOutput->Write() ;
 808   file->Close() ;
 809   delete file ;
 810
 811 }