PWG2/SPECTRA/aodV0.C

   1 #define aodV0_cxx
   2 // The class definition in aodV0.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 everytime 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("aodV0.C")
  22 // Root > T->Process("aodV0.C","some options")
  23 // Root > T->Process("aodV0.C+")
  24 //
  25
  26 #include "aodV0.h"
  27 #include <TPDGCode.h>
  28 #include <TSystem.h>
  29 #include <TStyle.h>
  30 #include <TCanvas.h>
  31 #include <TPaveLabel.h>
  32 #include <TPostScript.h>
  33
  34 #if !defined(__CINT__) || defined(__MAKECINT__)
  35 #include "AliESD.h"
  36 #include "AliPID.h"
  37 #include "AliAODv0.h"
  38 #endif
  39
  40
  41
  42 aodV0::aodV0(TTree *) :
  43   TSelector(),
  44   fChain(0),
  45   fESD(0),
  46   fHistV0PerEvent(0),
  47   fHistMassK0(0),
  48   fHistMassLambda(0),
  49   fHistMassAntiLambda(0),
  50   fHistMassLambdaVsProb(0),
  51   fHistMassLambdaCut(0),
  52   fHistMassAntiLambdaCut(0),
  53   fHistPtVsRapK0Short(0),
  54   fHistPtVsRapLambda(0),
  55   fHistArmenterosPodolanski(0)
  56 {
  57   // Constructor. Initialization of pointers
  58 }
  59
  60 aodV0::~aodV0() {
  61   // Remove all pointers
  62
  63   delete fESD;
  64
  65   // histograms are in the output list and deleted when the output
  66   // list is deleted by the TSelector dtor
  67 }
  68
  69 void aodV0::Begin(TTree *)
  70 {
  71   // The Begin() function is called at the start of the query.
  72   // When running with PROOF Begin() is only called on the client.
  73   // The tree argument is deprecated (on PROOF 0 is passed).
  74
  75   TString option = GetOption();
  76 }
  77
  78 void aodV0::SlaveBegin(TTree * tree)
  79 {
  80   // The SlaveBegin() function is called after the Begin() function.
  81   // When running with PROOF SlaveBegin() is called on each slave server.
  82   // The tree argument is deprecated (on PROOF 0 is passed).
  83
  84   Init(tree);
  85
  86   TString option = GetOption();
  87
  88   // Create histograms on each slave server
  89   fHistV0PerEvent = new TH1F("h1V0PerEvent", "V^{0} candidates per event", 50, 0, 50);
  90   fHistV0PerEvent->GetXaxis()->SetTitle("Number of V^{0}");
  91   fHistV0PerEvent->GetYaxis()->SetTitle("Events");
  92   fHistV0PerEvent->SetOption("HE");
  93   fHistV0PerEvent->SetStats(11);
  94
  95   fHistMassK0 = new TH1F("h1MassK0", "K^{0} candidates", 100, 0.4, 0.6);
  96   fHistMassK0->GetXaxis()->SetTitle("Invariant Mass_{#pi^{+}#pi^{-}}  (GeV/c^{2})");
  97   fHistMassK0->SetOption("HE");
  98   fHistMassK0->SetStats(11);
  99
 100   fHistMassLambda = new TH1F("h1MassLambda", "#Lambda^{0} candidates", 75, 1.05, 1.2);
 101   fHistMassLambda->GetXaxis()->SetTitle("Invariant Mass_{p#pi^{-}} (GeV/c^{2})");
 102   fHistMassLambda->SetOption("HE");
 103   fHistMassLambda->SetLineStyle(1);
 104   fHistMassLambda->SetStats(kFALSE);
 105
 106   fHistMassAntiLambda = new TH1F("h1MassAntiLambda", "#bar{#Lambda}^{0} candidates", 75, 1.05, 1.2);
 107   fHistMassAntiLambda->GetXaxis()->SetTitle("Invariant Mass_{#bar{p}#pi^{+}} (GeV/c^{2})");
 108   fHistMassAntiLambda->SetOption("HE");
 109   fHistMassAntiLambda->SetLineStyle(2);
 110   fHistMassAntiLambda->SetStats(kFALSE);
 111
 112   fHistMassLambdaVsProb = new TH2F("h2MassLambdaVsProb", "#Lambda^{0} and #bar{#Lambda}^{0} candidates", 21, -2.5, 102.5, 75, 1.05, 1.2);
 113   fHistMassLambdaVsProb->GetXaxis()->SetTitle("prob(p) (%)");
 114   fHistMassLambdaVsProb->GetYaxis()->SetTitle("Invariant Mass_{p#pi} (GeV/c^{2})");
 115   fHistMassLambdaVsProb->SetOption("BOX");
 116   fHistMassLambdaVsProb->SetFillStyle(0);
 117   fHistMassLambdaVsProb->SetStats(kFALSE);
 118
 119   fHistMassLambdaCut = new TH1F("h1MassLambdaCut", "#Lambda^{0} candidates", 75, 1.05, 1.2);
 120   fHistMassLambdaCut->GetXaxis()->SetTitle("Invariant Mass_{p#pi^{-}} (GeV/c^{2})");
 121   fHistMassLambdaCut->SetOption("E");
 122   fHistMassLambdaCut->SetMarkerStyle(kFullCircle);
 123   fHistMassLambdaCut->SetMarkerColor(kRed);
 124   fHistMassLambdaCut->SetStats(kFALSE);
 125
 126   fHistMassAntiLambdaCut = new TH1F("h1MassAntiLambdaCut", "#bar{#Lambda}^{0} candidates", 75, 1.05, 1.2);
 127   fHistMassAntiLambdaCut->GetXaxis()->SetTitle("Mass_{#bar{p}#pi^{+}} (GeV/c^{2})");
 128   fHistMassAntiLambdaCut->SetOption("E");
 129   fHistMassAntiLambdaCut->SetMarkerStyle(kFullCircle);
 130   fHistMassAntiLambdaCut->SetMarkerColor(kRed);
 131   fHistMassAntiLambdaCut->SetStats(kFALSE);
 132
 133   // AOD Histograms
 134   fHistPtVsRapK0Short = new TH2F("h2PtVsRapK0Short","K^{0}_{s} phase space",20,-1,1,20,0,10);
 135   fHistPtVsRapK0Short->SetOption("COL2Z");
 136   fHistPtVsRapK0Short->SetXTitle("Rapidity");
 137   fHistPtVsRapK0Short->SetYTitle("p_{t} (GeV/c)");
 138
 139   fHistPtVsRapLambda = new TH2F("h2PtVsRapLambda","#Lambda phase space",20,-1,1,20,0,10);
 140   fHistPtVsRapLambda->SetOption("COL2Z");
 141   fHistPtVsRapLambda->SetXTitle("Rapidity");
 142   fHistPtVsRapLambda->SetYTitle("p_{t} (GeV/c)");
 143
 144   fHistArmenterosPodolanski = new TH2F("h2ArmenterosPodolanski","Armenteros-Podolanski phase space",100,-1.0,1.0,50,0,0.5);
 145   fHistArmenterosPodolanski->SetOption("BOX");
 146   fHistArmenterosPodolanski->SetXTitle("#alpha");
 147   fHistArmenterosPodolanski->SetYTitle("p_{t} arm");
 148 }
 149
 150 void aodV0::Init(TTree *tree)
 151 {
 152   // The Init() function is called when the selector needs to initialize
 153   // a new tree or chain. Typically here the branch addresses of the tree
 154   // will be set. It is normaly not necessary to make changes to the
 155   // generated code, but the routine can be extended by the user if needed.
 156   // Init() will be called many times when running with PROOF.
 157
 158   // Set branch addresses
 159   if (tree == 0) return;
 160   fChain = tree;
 161
 162   fChain->SetBranchAddress("ESD",&fESD);
 163 }
 164
 165 Bool_t aodV0::Notify()
 166 {
 167   // The Notify() function is called when a new file is opened. This
 168   // can be either for a new TTree in a TChain or when when a new TTree
 169   // is started when using PROOF. Typically here the branch pointers
 170   // will be retrieved. It is normaly not necessary to make changes
 171   // to the generated code, but the routine can be extended by the
 172   // user if needed.
 173
 174   return kTRUE;
 175 }
 176
 177
 178 Bool_t aodV0::Process(Long64_t entry)
 179 {
 180   // The Process() function is called for each entry in the tree (or possibly
 181   // keyed object in the case of PROOF) to be processed. The entry argument
 182   // specifies which entry in the currently loaded tree is to be processed.
 183   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
 184   // to read either all or the required parts of the data. When processing
 185   // keyed objects with PROOF, the object is already loaded and is available
 186   // via the fObject pointer.
 187   //
 188   // This function should contain the "body" of the analysis. It can contain
 189   // simple or elaborate selection criteria, run algorithms on the data
 190   // of the event and typically fill histograms.
 191
 192   // WARNING when a selector is used with a TChain, you must use
 193   //  the pointer to the current TTree to call GetEntry(entry).
 194   //  The entry is always the local entry number in the current tree.
 195   //  Assuming that fChain is the pointer to the TChain being processed,
 196   //  use fChain->GetEntry(entry).
 197
 198   fChain->GetTree()->GetEntry(entry);
 199
 200   if (!fESD) return kFALSE;
 201
 202   Int_t mBoDebug = 1;
 203
 204   // Primary vertex and Position
 205   const AliESDVertex *lPrimaryVertex = fESD->GetVertex();
 206   Double_t tPositionPrimaryVertex[3]; lPrimaryVertex->GetXYZ(tPositionPrimaryVertex);
 207   if (mBoDebug) printf("***BoInfo: Primary Vertex position x=%.3f, y=%.3f, z=%.3f \n",tPositionPrimaryVertex[0],tPositionPrimaryVertex[1],tPositionPrimaryVertex[2]);
 208
 209   // Following is commented since ntracks is not used (remove warning)
 210   //  Int_t ntracks = fESD->GetNumberOfTracks();
 211   Int_t nv0s = fESD->GetNumberOfV0s();
 212   if (mBoDebug) printf("***BoInfo: Number of V0s =%d \n",nv0s);
 213   fHistV0PerEvent->Fill(nv0s);
 214
 215   Double_t mass;
 216   Double_t pPrior[5] = {0, 0, 0.85, 0.1, 0.05};
 217   Double_t p[5];
 218   Double_t pSum, pNorm;
 219
 220   AliAODv0 *myAODv0 = new AliAODv0();
 221
 222   // variable for AOD v0
 223   Double_t dcaPosToPrimVertex = 0, dcaNegToPrimVertex = 0;
 224   Double_t dcaV0Daughters     = 0, decayLengthV0      = 0;
 225   Double_t alphaV0            = 0, ptArmV0            = 0;
 226   Double_t rapK0Short         = 0, rapLambda          = 0;
 227   Double_t pt                 = 0;
 228
 229   for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
 230     AliESDv0* v0 = fESD->GetV0(iV0);
 231     if (!v0) continue;
 232     myAODv0->ResetV0();
 233     myAODv0->Fill(v0,fESD);
 234
 235     dcaPosToPrimVertex = myAODv0->DcaPosToPrimVertex();
 236     dcaNegToPrimVertex = myAODv0->DcaNegToPrimVertex();
 237     dcaV0Daughters     = myAODv0->DcaV0Daughters();
 238     decayLengthV0      = myAODv0->DecayLengthV0(tPositionPrimaryVertex);
 239     alphaV0            = myAODv0->AlphaV0();
 240     ptArmV0            = myAODv0->PtArmV0();
 241     rapK0Short         = myAODv0->RapK0Short();
 242     rapLambda          = myAODv0->RapLambda();
 243     pt                 = TMath::Sqrt(myAODv0->Ptot2V0());
 244
 245     fHistPtVsRapK0Short->Fill(rapK0Short,pt);
 246     fHistPtVsRapLambda->Fill(rapLambda,pt);
 247     fHistArmenterosPodolanski->Fill(alphaV0,ptArmV0);
 248
 249     fHistMassK0->Fill(v0->GetEffMass());
 250
 251     v0->ChangeMassHypothesis(kLambda0);
 252     mass = v0->GetEffMass();
 253     fESD->GetTrack(v0->GetPindex())->GetESDpid(p);
 254     pSum = 0;
 255     for (Int_t i = 0; i < AliPID::kSPECIES; i++) pSum += p[i] * pPrior[i];
 256     if (pSum <= 0) pSum = 1.;
 257     pNorm = p[AliPID::kProton] * pPrior[AliPID::kProton] / pSum;
 258     fHistMassLambdaVsProb->Fill(100.*pNorm, mass);
 259     if (pNorm > 0.1) fHistMassLambda->Fill(mass);
 260
 261     v0->ChangeMassHypothesis(kLambda0Bar);
 262     mass = v0->GetEffMass();
 263     fESD->GetTrack(v0->GetNindex())->GetESDpid(p);
 264     pSum = 0;
 265     for (Int_t i = 0; i < AliPID::kSPECIES; i++) pSum += p[i] * pPrior[i];
 266     if (pSum <= 0) pSum = 1.;
 267     pNorm = p[AliPID::kProton] * pPrior[AliPID::kProton] / pSum;
 268     fHistMassLambdaVsProb->Fill(100.*pNorm, mass);
 269     if (pNorm > 0.1) fHistMassAntiLambda->Fill(mass);
 270
 271   }//loop v0s
 272
 273   return kTRUE;
 274 }
 275
 276
 277
 278 void aodV0::SlaveTerminate()
 279 {
 280   // The SlaveTerminate() function is called after all entries or objects
 281   // have been processed. When running with PROOF SlaveTerminate() is called
 282   // on each slave server.
 283
 284   // Add the histograms to the output on each slave server
 285   fOutput->Add(fHistV0PerEvent);
 286   fOutput->Add(fHistMassK0);
 287   fOutput->Add(fHistMassLambda);
 288   fOutput->Add(fHistMassAntiLambda);
 289   fOutput->Add(fHistMassLambdaVsProb);
 290   fOutput->Add(fHistMassLambdaCut);
 291   fOutput->Add(fHistMassAntiLambdaCut);
 292   // Add the AOD histograms to the output on each slave server
 293   fOutput->Add(fHistPtVsRapK0Short);
 294   fOutput->Add(fHistPtVsRapLambda);
 295   fOutput->Add(fHistArmenterosPodolanski);
 296 }
 297
 298 void aodV0::Terminate()
 299 {
 300   // The Terminate() function is the last function to be called during
 301   // a query. It always runs on the client, it can be used to present
 302   // the results graphically or save the results to file.
 303
 304   fHistV0PerEvent = dynamic_cast<TH1F*>(fOutput->FindObject("h1V0PerEvent"));
 305   fHistMassK0 = dynamic_cast<TH1F*>(fOutput->FindObject("h1MassK0"));
 306   fHistMassLambda = dynamic_cast<TH1F*>(fOutput->FindObject("h1MassLambda"));
 307   fHistMassAntiLambda = dynamic_cast<TH1F*>(fOutput->FindObject("h1MassAntiLambda"));
 308   fHistMassLambdaVsProb = dynamic_cast<TH2F*>(fOutput->FindObject("h2MassLambdaVsProb"));
 309   fHistMassLambdaCut = dynamic_cast<TH1F*>(fOutput->FindObject("h1MassLambdaCut"));
 310   fHistMassAntiLambdaCut = dynamic_cast<TH1F*>(fOutput->FindObject("h1MassAntiLambdaCut"));
 311
 312   // AOD Histograms
 313   fHistPtVsRapK0Short = dynamic_cast<TH2F*>(fOutput->FindObject("h2PtVsRapK0Short"));
 314   fHistPtVsRapLambda = dynamic_cast<TH2F*>(fOutput->FindObject("h2PtVsRapLambda"));
 315   fHistArmenterosPodolanski = dynamic_cast<TH2F*>(fOutput->FindObject("h2ArmenterosPodolanski"));
 316
 317   TFile* file = TFile::Open("aodV0.root", "RECREATE");
 318   fHistV0PerEvent->Write();
 319   fHistMassK0->Write();
 320   fHistMassLambda->Write();
 321   fHistMassAntiLambda->Write();
 322   fHistMassLambdaVsProb->Write();
 323   fHistMassLambdaCut->Write();
 324   fHistMassAntiLambdaCut->Write();
 325   // AOD Histograms
 326   fHistPtVsRapK0Short->Write();
 327   fHistPtVsRapLambda->Write();
 328   fHistArmenterosPodolanski->Write();
 329
 330   file->Close();
 331   delete file;
 332
 333   // Define postscript
 334
 335
 336   if (!gROOT->IsBatch())
 337     {
 338       TCanvas *canOnline = new TCanvas("canOnline","V0s",10,10,610,610);
 339       canOnline->SetFillColor(10);
 340       canOnline->SetHighLightColor(10);
 341       canOnline->Divide(2,1);
 342
 343       canOnline->cd(1);
 344       fHistMassK0->DrawCopy("E");
 345       fHistMassK0->Fit("gaus","q","",0.49,0.51);
 346       TVirtualPad* pad = (TVirtualPad*)canOnline->cd(2);
 347       pad->Divide(1,2);
 348       pad->cd(1);
 349       fHistMassLambda->DrawCopy("E");
 350       pad->cd(2);
 351       fHistMassAntiLambda->DrawCopy("E");
 352     }
 353
 354   // Following define a postscript output
 355   int  myAliceRed        = TColor::GetColor(192,0,0);
 356   TDatime *currentDatime = new TDatime();
 357   Char_t info[125];
 358   gStyle->SetPaperSize(20,26);
 359 //   gStyle->SetOptStat(111111);
 360   gStyle->SetOptStat(11);
 361   gStyle->SetPalette(1,0);
 362   TCanvas     *c0 = new TCanvas("c0","AOD V0 Analysis",600,800);
 363   c0->Range(0,0,20,24);
 364
 365   TPostScript *ps = new TPostScript("aodV0.ps");
 366
 367   TPaveLabel  *t0 = new TPaveLabel(4,23.0,16,23.8,"AOD V0 Analysis","br");
 368   t0->SetFillColor(0);
 369   t0->SetBorderSize(2);
 370   t0->SetTextColor(myAliceRed);
 371
 372   TPaveLabel  *d0 = new TPaveLabel(16.5,23.0,19.5,23.6,currentDatime->AsSQLString(),"br");
 373   d0->SetFillColor(0);
 374   d0->SetTextSize(0.4);
 375   d0->SetBorderSize(0);
 376   d0->SetTextColor(1);
 377
 378   TPaveText *ptt1 = new TPaveText(1.5,20.6,8.0,21.4);
 379   ptt1->SetFillColor(0);
 380   ptt1->SetBorderSize(2);
 381   ptt1->SetTextSize(0.03);
 382   ptt1->SetTextColor(myAliceRed);
 383   ptt1->SetTextAlign(12);
 384   ptt1->AddText("Global Information:");
 385
 386   TPaveText *ptx1 = new TPaveText(2,16,18,21);
 387   ptx1->SetFillColor(0);
 388   ptx1->SetBorderSize(1);
 389   ptx1->SetTextSize(0.02);
 390   ptx1->SetTextAlign(13);
 391   ptx1->AddText("");
 392   TString USER("$USER");
 393   TString HOST("$HOST");
 394   TString ROOTSYS("$ROOTSYS");
 395   TString ALICE_LEVEL("$ALICE_LEVEL");
 396   gSystem->ExpandPathName(USER);
 397   gSystem->ExpandPathName(HOST);
 398   gSystem->ExpandPathName(ROOTSYS);
 399   gSystem->ExpandPathName(ALICE_LEVEL);
 400   sprintf(info,"Analysis from %s on %s",USER.Data(),HOST.Data());
 401   ptx1->AddText(info);
 402   sprintf(info,"Root Version: %s",ROOTSYS.Data());
 403   ptx1->AddText(info);
 404   sprintf(info,"Alice Version: %s",ALICE_LEVEL.Data());
 405   ptx1->AddText(info);
 406   sprintf(info,"Number of scanned events : %.0f",fHistV0PerEvent->GetEntries());
 407   ptx1->AddText(info);
 408
 409   ps->NewPage();
 410   t0->Draw();
 411   d0->Draw();
 412
 413   ptx1->Draw();
 414   ptt1->Draw();
 415   c0->Update();
 416
 417   delete ptx1;
 418   delete ptt1;
 419
 420   delete t0;
 421   delete c0;
 422
 423   // Rich Primary Vertex
 424   ps->NewPage();
 425   TCanvas* c1 = new TCanvas("c1","Phase Space Page",500,690);
 426   c1->Update();
 427   c1->cd();
 428   c1->Range(0,0,25,18);
 429
 430   TPaveLabel *ptitle1 = new TPaveLabel(5,17,20,18,"Phase Space Info","br");
 431   ptitle1->SetFillColor(18);
 432   ptitle1->SetTextFont(32);
 433   ptitle1->SetTextSize(0.8);
 434   ptitle1->SetTextColor(myAliceRed);
 435   ptitle1->Draw();
 436
 437   TPad    *c1Pad11 = new TPad("c1Pad11","box",0.01,0.50,0.33,0.90,0);
 438   TPad    *c1Pad12 = new TPad("c1Pad12","box",0.34,0.50,0.66,0.90,0);
 439   TPad    *c1Pad13 = new TPad("c1Pad13","box",0.67,0.50,0.99,0.90,0);
 440   TPad    *c1Pad21 = new TPad("c1Pad21","box",0.01,0.05,0.33,0.45,0);
 441   TPad    *c1Pad22 = new TPad("c1Pad22","box",0.34,0.05,0.66,0.45,0);
 442   TPad    *c1Pad23 = new TPad("c1Pad23","box",0.67,0.05,0.99,0.45,0);
 443   c1Pad11->SetLeftMargin(0.15); c1Pad11->SetBottomMargin(0.15);
 444   c1Pad11->Draw();
 445   c1Pad12->SetLeftMargin(0.15); c1Pad12->SetBottomMargin(0.15);
 446   c1Pad12->Draw();
 447   c1Pad13->SetLeftMargin(0.15); c1Pad13->SetBottomMargin(0.15);
 448   c1Pad13->Draw();
 449   c1Pad21->SetLeftMargin(0.15); c1Pad21->SetBottomMargin(0.15);
 450   c1Pad21->Draw();
 451   c1Pad22->SetLeftMargin(0.15); c1Pad22->SetBottomMargin(0.15);
 452   c1Pad22->Draw();
 453   c1Pad23->SetLeftMargin(0.15); c1Pad23->SetBottomMargin(0.15);
 454   c1Pad23->Draw();
 455
 456   c1Pad11->cd();
 457   fHistV0PerEvent->Draw("H");
 458   c1Pad12->cd();
 459   fHistMassK0->Draw("HE");
 460   c1Pad13->cd();
 461   fHistMassLambda->Draw("HE");
 462   fHistMassAntiLambda->SetLineStyle(2);
 463   fHistMassAntiLambda->Draw("HESAME");
 464   c1Pad13->Update();
 465
 466   gPad->Update();
 467   c1Pad21->cd();
 468   fHistPtVsRapK0Short->Draw("COL2Z");
 469   c1Pad22->cd();
 470   fHistPtVsRapLambda->Draw("COL2Z");
 471   c1Pad23->cd();
 472   fHistArmenterosPodolanski->Draw("BOX");
 473   c1->Update();
 474   ps->NewPage();
 475
 476   delete c1Pad11;
 477   delete c1Pad12;
 478   delete c1Pad13;
 479   delete c1Pad21;
 480   delete c1Pad22;
 481   delete c1Pad23;
 482   delete c1;
 483
 484   ps->Close();
 485 }