PWG2/FORWARD/analysis2/SPDComparison.C

   1 /**
   2  * @mainpage Forward/Backward Correlations
   3  *
   4  * A script containing a class ForwardBackwardTask and a
   5  * function to run the analysis.
   6  *
   7  * The class ForwardBackwardTask is an AliAnalysisTaskSE.  That means
   8  * that it have facilities for analysing ESD, AOD, and MC input.  The
   9  * process of running the code is handled by an AliAnalysisManager
  10  * (created in the function).  It uses a TSelector to loop over a
  11  * TChain of data.
  12  *
  13  * The flow of the code is
  14  * @verbatim
  15  *    +-----------------------+
  16  *    | Create Analysis Train |-> ForwardBackwardTask constructor
  17  *    +-----------------------+
  18  *                |
  19  *                V
  20  *    +-----------------------+
  21  *    | Intialise all tasks   |-> ForwardBackwardTask::Init (not implemented)
  22  *    +-----------------------+
  23  *                |
  24  *                V
  25  *    +-----------------------+
  26  *    | Split job on workers  |
  27  *    +-----------------------+
  28  *                |
  29  *                V
  30  *    +-----------------------+
  31  *    | Create output objects |-> ForwardBackwardTask::CreateOutputObjects
  32  *    | on each worker        |
  33  *    +-----------------------+
  34  *                |
  35  *                V
  36  *    +-----------------------+
  37  *    | More events on this   |<-----+
  38  *    | worker node?          |--+   |
  39  *    +-----------------------+  |   |
  40  *                | no           |   |
  41  *                |              V   |
  42  *                |   +-------------------+
  43  *                |   | Process one event |->ForwardBackwardTask::UserExec
  44  *                |   +-------------------+
  45  *                |
  46  *                V
  47  *    +-----------------------+
  48  *    | Merge output of each  |
  49  *    | worker node           |
  50  *    +-----------------------+
  51  *                |
  52  *                V
  53  *    +-----------------------+
  54  *    | End of job processing |-> ForwardBackwardTask::Terminate
  55  *    +-----------------------+
  56  * @endverbatim
  57  *
  58  * Since the class ForwardBackwardTask derives from a compiled class
  59  * (AliAnalysisTaskSE) we need to compile that code.  The script will,
  60  * when executed in the AliROOT prompt load it self again and byte
  61  * compile it with the preprocessor flag BUILD defined.  \
  62  *
  63  * The preprocessor define BUILD is not defined at first, when the
  64  * script is loaded using
  65  *
  66  * @verbatim
  67  *   Root> .x ForwardBackward.C
  68  * @endverbatim
  69  *
  70  * which means that CINT will only see the function ForwardBackward.
  71  * In that function, we define the BUILD preprocessor symbol
  72  *
  73  * @code
  74  *   gSystem->AddIncludePath("-DBUILD=1 ...");
  75  * @endcode
  76  *
  77  * and then ACLic compile ourselves
  78  *
  79  * @code
  80  *   gROOT->LoadMacro("./SPDComparison.C++");
  81  * @endcode
  82  *
  83  * But since BUILD is now defined, it means that ACLic will only see
  84  * the class and not the function SPDComparison.
  85  *
  86  * This trick hinges on that when you initially load the script and
  87  * when it is done inside the script it is done using two distinct
  88  * paths - otherwise ROOT will try to unload the script first, and
  89  * that fails.
  90  *
  91  */
  92
  93 #ifdef BUILD
  94 #include <AliAnalysisTaskSE.h>
  95 #include <AliESDEvent.h>
  96 #include <AliLog.h>
  97 #include <AliMultiplicity.h>
  98 #include "AliFMDEventInspector.h"
  99 #include "AliAODForwardMult.h"
 100 #include <TH1D.h>
 101 #include <TH2D.h>
 102 #include <TGraphErrors.h>
 103 #include <TList.h>
 104 #include <TObjArray.h>
 105 #include <TMath.h>
 106
 107
 108 /**
 109  * @class SPDComparisonTask
 110  *
 111  * Task to do forward/backward correlations using FMD and SPD
 112  * (cluster) data.
 113  *
 114  * The class contains the sub-structure Bin that represents an
 115  * @f$\eta@f$ range.  One can add (possibly overlapping) @f$\eta@f$
 116  * ranges by calling the member function AddBin
 117  *
 118  */
 119 class SPDComparisonTask : public AliAnalysisTaskSE
 120 {
 121 public:
 122   /**
 123    * A single vertex bin - contains summed histogram of clusters and
 124    * tracklets
 125    *
 126    */
 127   struct VtxBin : public TNamed
 128   {
 129     /**
 130      * Static function to get a bin name
 131      *
 132      * @param low   Low limit
 133      * @param high  High limit
 134      *
 135      * @return Pointer to static array
 136      */
 137     static const char* BinName(Double_t low, Double_t high)
 138     {
 139       static TString n;
 140       n = (Form("vtx%+05.1f_%+05.1f", low, high));
 141       n.ReplaceAll("+", "p");
 142       n.ReplaceAll("-", "m");
 143       n.ReplaceAll(".", "d");
 144       return n.Data();
 145     }
 146     /**
 147      * Constructor
 148      *
 149      * @param low   Low limit
 150      * @param high  High limit
 151      */
 152     VtxBin(Double_t low, Double_t high)
 153       : TNamed(BinName(low,high), ""),
 154         fClusters(0),
 155         fTracklets(0)
 156     {
 157     }
 158     /**
 159      * Constructor
 160      */
 161     VtxBin() : TNamed(), fClusters(0), fTracklets(0) {}
 162     /**
 163      * Copy constructor
 164      *
 165      * @param o Object to copy from
 166      */
 167     VtxBin(const VtxBin& o) :
 168       TNamed(o), fClusters(o.fClusters), fTracklets(o.fTracklets)
 169     {}
 170     /**
 171      * Assignment operator
 172      *
 173      * @param o Object to assign from
 174      *
 175      * @return Reference to this.
 176      */
 177     VtxBin& operator=(const VtxBin& o)
 178     {
 179       TNamed::operator=(o);
 180       fClusters = o.fClusters;
 181       fTracklets = o.fTracklets;
 182       return *this;
 183     }
 184     /**
 185      * Define outputs
 186      *
 187      * @param l    Output list
 188      * @param eta  Eta axis to use
 189      */
 190     void DefineOutput(TList* l, const TAxis& eta)
 191     {
 192       TList* ll = new TList;
 193       ll->SetName(GetName());
 194       ll->SetOwner();
 195       l->Add(ll);
 196
 197       fClusters = new TH1D("clusters", "Clusters",
 198                            eta.GetNbins(),
 199                            eta.GetXmin(),
 200                            eta.GetXmax());
 201       fClusters->SetXTitle("#eta");
 202       fClusters->SetYTitle("dN/d#eta");
 203       fClusters->SetDirectory(0);
 204       fClusters->Sumw2();
 205       ll->Add(fClusters);
 206
 207       fTracklets = static_cast<TH1D*>(fClusters->Clone("tracklets"));
 208       fTracklets->SetTitle("Tracklets");
 209       fTracklets->SetDirectory(0);
 210       ll->Add(fTracklets);
 211     }
 212     /**
 213      * Process the input
 214      *
 215      * @param spdmult  Input
 216      */
 217     void Process(const AliMultiplicity* spdmult)
 218     {
 219       //Filling clusters in layer 1 used for tracklets...
 220       for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) {
 221         Double_t eta = spdmult->GetEta(j);
 222         fClusters->Fill(eta);
 223         fTracklets->Fill(eta);
 224       }
 225
 226       //...and then the unused ones in layer 1
 227       for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
 228         Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
 229         fClusters->Fill(eta);
 230       }
 231     }
 232     /**
 233      * End of job processing
 234      *
 235      * @param l        Output list
 236      * @param nEvents  Number of events to normalise to
 237      */
 238     void Finish(TList* l, Int_t nEvents)
 239     {
 240       TList* ll  = static_cast<TList*>(l->FindObject(GetName()));
 241       fClusters  = static_cast<TH1D*>(ll->FindObject("clusters"));
 242       fTracklets = static_cast<TH1D*>(ll->FindObject("tracklets"));
 243       fClusters->Scale(1. / nEvents, "width");
 244       fTracklets->Scale(1. / nEvents, "width");
 245       TH1* ratio = static_cast<TH1D*>(fClusters->Clone("ratio"));
 246       ratio->SetDirectory(0);
 247       ratio->SetTitle("Clusters/Tracklets");
 248       ratio->Divide(fTracklets);
 249       ll->Add(ratio);
 250     }
 251     TH1D* fClusters;  // Cluster histogram
 252     TH1D* fTracklets; // Tracklet histogram
 253   };
 254
 255   //__________________________________________________________________
 256   /**
 257    * Default constructor - do not use
 258    */
 259   SPDComparisonTask()
 260     : AliAnalysisTaskSE(),
 261       fInspector(),
 262       fBins(0),
 263       fEvents(0),
 264       fList(0),
 265       fFirstEvent(true),
 266       fVertexAxis(20, -20, 20)
 267   {}
 268   /**
 269    * Constructor with a single argument.
 270    */
 271   SPDComparisonTask(const char*)
 272     : AliAnalysisTaskSE("spdComparision"),
 273       fInspector("eventInspector"),
 274       fBins(0),
 275       fEvents(0),
 276       fList(0),
 277       fFirstEvent(true),
 278       fVertexAxis(20, -20, 20)
 279   {
 280     // Declare our output container
 281     DefineOutput(1,TList::Class());
 282   }
 283   /**
 284    * Copy constructor
 285    *
 286    * @param o Object to copy from
 287    */
 288   SPDComparisonTask(const SPDComparisonTask& o)
 289     : AliAnalysisTaskSE(o),
 290       fInspector(o.fInspector),
 291       fBins(o.fBins),
 292       fEvents(o.fEvents),
 293       fList(o.fList),
 294       fFirstEvent(o.fFirstEvent),
 295       fVertexAxis(20, -20, 20)
 296   {
 297     SetVertexAxis(o.fVertexAxis);
 298   }
 299   /**
 300    * Assigment operator
 301    *
 302    * @param o Object to assign from
 303    *
 304    * @return Reference to this
 305    */
 306   SPDComparisonTask& operator=(const SPDComparisonTask& o)
 307   {
 308     AliAnalysisTaskSE::operator=(o);
 309     fInspector = o.fInspector;
 310     fEvents    = o.fEvents;
 311     fList      = o.fList;
 312     fFirstEvent = o.fFirstEvent;
 313     SetVertexAxis(o.fVertexAxis);
 314     return *this;
 315   }
 316   /**
 317    * Destructor
 318    */
 319   ~SPDComparisonTask() {}
 320   /**
 321    * Set the vertex axis to use
 322    *
 323    * @param a Axis to set from
 324    */
 325   void SetVertexAxis(const TAxis& a)
 326   {
 327     SetVertexAxis(a.GetNbins(),
 328                   a.GetXmin(),
 329                   a.GetXmax());
 330   }
 331   /**
 332    * Set the vertex axis to use
 333    *
 334    * @param n     Number of bins
 335    * @param xmin  Least @f$v_z@f$
 336    * @param xmax  Most @f$v_z@f$
 337    */
 338   void SetVertexAxis(Int_t n, Double_t xmin, Double_t xmax)
 339   {
 340     fVertexAxis.Set(n, xmin, xmax);
 341   }
 342   /**
 343    * Create output objects
 344    *
 345    */
 346   void UserCreateOutputObjects()
 347   {
 348     fList = new TList;
 349     fList->SetName(GetName());
 350     fList->SetOwner();
 351
 352
 353     fEvents = new TH1D("events", "Events",
 354                        fVertexAxis.GetNbins(),
 355                        fVertexAxis.GetXmin(),
 356                        fVertexAxis.GetXmax());
 357     fEvents->SetDirectory(0);
 358     fEvents->SetXTitle("v_{z} [cm]");
 359     fList->Add(fEvents);
 360
 361     TAxis& vtxAxis = fVertexAxis;
 362     fBins = new TObjArray(vtxAxis.GetNbins(),1);
 363
 364     TAxis etaAxis(120, -3, 3);
 365     for (Int_t i = 1; i <= vtxAxis.GetNbins(); i++) {
 366       VtxBin* bin = new VtxBin(vtxAxis.GetBinLowEdge(i),
 367                                vtxAxis.GetBinUpEdge(i));
 368       bin->DefineOutput(fList, etaAxis);
 369       fBins->AddAt(bin,i);
 370     }
 371
 372
 373     fInspector.DefineOutput(fList);
 374     fInspector.Init(*(fEvents->GetXaxis()));
 375
 376     PostData(1, fList);
 377   }
 378   /**
 379    * Process one event
 380    *
 381    * @param option Not used
 382    */
 383   void UserExec(Option_t* option="")
 384   {
 385     // Get the input data - ESD event
 386     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 387     if (!esd) {
 388       AliWarning("No ESD event found for input event");
 389       return;
 390     }
 391
 392     //--- Read run information -----------------------------------------
 393     if (fFirstEvent && esd->GetESDRun()) {
 394       fInspector.ReadRunDetails(esd);
 395
 396       AliInfo(Form("Initializing with parameters from the ESD:\n"
 397                    "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 398                    "         AliESDEvent::GetBeamType()     ->%s\n"
 399                    "         AliESDEvent::GetCurrentL3()    ->%f\n"
 400                    "         AliESDEvent::GetMagneticField()->%f\n"
 401                    "         AliESDEvent::GetRunNumber()    ->%d\n",
 402                    esd->GetBeamEnergy(),
 403                    esd->GetBeamType(),
 404                    esd->GetCurrentL3(),
 405                    esd->GetMagneticField(),
 406                    esd->GetRunNumber()));
 407
 408       Print();
 409       fFirstEvent = false;
 410     }
 411
 412     // Some variables
 413     UInt_t   triggers; // Trigger bits
 414     Bool_t   lowFlux;  // Low flux flag
 415     UShort_t iVz;      // Vertex bin from ESD
 416     Double_t vZ;       // Z coordinate from ESD
 417     Double_t cent;     // Centrality
 418     UShort_t nClusters;// Number of SPD clusters
 419     // Process the data
 420     UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ,
 421                                        cent, nClusters);
 422     Bool_t isInel   = triggers & AliAODForwardMult::kB;
 423     Bool_t hasVtx   = retESD == AliFMDEventInspector::kOk;
 424
 425     if (!isInel || !hasVtx) return;
 426     fEvents->Fill(vZ);
 427
 428     const AliMultiplicity* spdmult = esd->GetMultiplicity();
 429
 430     VtxBin* bin = static_cast<VtxBin*>(fBins->At(iVz));
 431     if (!bin) {
 432       AliError(Form("No bin @ %d (%fcm)", iVz, vZ));
 433       return;
 434     }
 435     bin->Process(spdmult);
 436
 437     // Post data to output container
 438     PostData(1,fList);
 439   }
 440   /**
 441    * Called at the end of the processing
 442    *
 443    * @param option
 444    */
 445   void Terminate(Option_t* option="")
 446   {
 447     fList = dynamic_cast<TList*>(GetOutputData(1));
 448     if (!fList) {
 449       AliError("No output list defined");
 450       return;
 451     }
 452
 453     fEvents    = static_cast<TH1D*>(fList->FindObject("events"));
 454
 455     Int_t nEvents = fEvents->GetEntries();
 456     AliInfo(Form("Got a total of %d events", nEvents));
 457
 458     TH1* clusters = 0;
 459     TH1* tracklets = 0;
 460     TIter next(fBins);
 461     VtxBin* bin = 0;
 462     Int_t i = 1;
 463     while ((bin = static_cast<VtxBin*>(next()))) {
 464       bin->Finish(fList, fEvents->GetBinContent(i++));
 465       if (!clusters) clusters = static_cast<TH1D*>(bin->fClusters->Clone());
 466       else clusters->Add(bin->fClusters);
 467       if (!tracklets) tracklets = static_cast<TH1D*>(bin->fTracklets->Clone());
 468       else tracklets->Add(bin->fTracklets);
 469     }
 470     clusters->SetDirectory(0);
 471     tracklets->SetDirectory(0);
 472     clusters->Scale(1. / i);
 473     tracklets->Scale(1. / i);
 474
 475
 476     TH1D* ratio = static_cast<TH1D*>(clusters->Clone("ratio"));
 477     ratio->SetDirectory(0);
 478     ratio->SetTitle("Clusters/Tracklets");
 479     ratio->Divide(tracklets);
 480
 481     fList->Add(clusters);
 482     fList->Add(tracklets);
 483     fList->Add(ratio);
 484   }
 485 protected:
 486   AliFMDEventInspector fInspector; // Inspector
 487   TObjArray*           fBins;      // Vertex bins
 488   TH1D*                fEvents;    // Events
 489   TList*               fList;      // Output list
 490   Bool_t               fFirstEvent;// First event flag
 491   TAxis                fVertexAxis;
 492
 493   ClassDef(SPDComparisonTask,1); // BF analysis
 494 };
 495 #else
 496
 497 //====================================================================
 498 /**
 499  * Run the analysis
 500  *
 501  * @param aoddir  Input directory
 502  * @param nEvents Number of events, negative means all
 503  */
 504 void
 505 SPDComparison(const char* esddir, Int_t nEvents=-1)
 506 {
 507   // --- Libraries to load -------------------------------------------
 508   gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
 509
 510   // --- Our data chain ----------------------------------------------
 511   gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeChain.C");
 512   TChain* chain = MakeChain("ESD", esddir, true);
 513   // If 0 or less events is select, choose all
 514   if (nEvents <= 0) nEvents = chain->GetEntries();
 515
 516   // --- Manager -----------------------------------------------------
 517   AliAnalysisManager* mgr = new AliAnalysisManager("FB", "FB train");
 518   AliAnalysisManager::SetCommonFileName("spd_comps.root");
 519
 520   // --- AOD input handler -------------------------------------------
 521   AliESDInputHandler *inputHandler = new AliESDInputHandler();
 522   mgr->SetInputEventHandler(inputHandler);
 523
 524   // --- compile our code --------------------------------------------
 525   gSystem->AddIncludePath("-I${ALICE_ROOT}/PWG2/FORWARD/analysis2 "
 526                           "-I${ALICE_ROOT}/ANALYSIS "
 527                           "-I${ALICE_ROOT}/include -DBUILD=1");
 528   gROOT->LoadMacro("./SPDComparison.C++g");
 529
 530   // --- Make our object ---------------------------------------------
 531   SPDComparisonTask* task = new SPDComparisonTask("SPD_COMP");
 532   task->SetVertexAxis(10, -10, 10);
 533   mgr->AddTask(task);
 534
 535   // --- create containers for input/output --------------------------
 536   AliAnalysisDataContainer *sums =
 537     mgr->CreateContainer("spdComp", TList::Class(),
 538                          AliAnalysisManager::kOutputContainer,
 539                          AliAnalysisManager::GetCommonFileName());
 540   // --- connect input/output ----------------------------------------
 541   mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
 542   mgr->ConnectOutput(task, 1, sums);
 543
 544   // --- Run the analysis --------------------------------------------
 545   TStopwatch t;
 546   if (!mgr->InitAnalysis()) {
 547     Error("SPDComparison", "Failed to initialize analysis train!");
 548     return;
 549   }
 550   // Some informative output
 551   mgr->PrintStatus();
 552   // mgr->SetDebugLevel(3);
 553   if (mgr->GetDebugLevel() < 1) mgr->SetUseProgressBar(kTRUE,100);
 554
 555   // Run the train
 556   t.Start();
 557   Printf("=== RUNNING ANALYSIS on %9d events ==================", nEvents);
 558   mgr->StartAnalysis("local", chain, nEvents);
 559   t.Stop();
 560   t.Print();
 561
 562   DrawSPDComparison();
 563 }
 564
 565 //====================================================================
 566 /**
 567  * Draw results
 568  *
 569  * @param filename
 570  */
 571 void
 572 DrawSPDComparison(const char* filename="spd_comps.root")
 573 {
 574   // --- Open the file -----------------------------------------------
 575   TFile* file = TFile::Open(filename, "READ");
 576   if (!file) {
 577     Error("DrawSPDComparison", "Failed to open file %s", filename);
 578     return;
 579   }
 580
 581   // --- Get our list ------------------------------------------------
 582   TList* spd = static_cast<TList*>(file->Get("spdComp"));
 583   if (!spd) {
 584     Error("DrawSPDComparison", "Failed to get list SPD_COMP from file %s",
 585           filename);
 586     return;
 587   }
 588
 589   // --- Loop over list and get correlation plots --------------------
 590   TH1* clusters  = static_cast<TH1*>(spd->FindObject("clusters"));
 591   TH1* tracklets = static_cast<TH1*>(spd->FindObject("tracklets"));
 592   TH1* ratio     = static_cast<TH1*>(spd->FindObject("ratio"));
 593   TH1* events    = static_cast<TH1*>(spd->FindObject("events"));
 594
 595   // --- Set style parameters ----------------------------------------
 596   gStyle->SetPalette(1);
 597   gStyle->SetTitleFillColor(0);
 598   gStyle->SetTitleStyle(0);
 599   gStyle->SetTitleBorderSize(0);
 600   gStyle->SetOptStat(0);
 601   gStyle->SetTitleX(0.69);
 602   gStyle->SetTitleY(0.99);
 603   gStyle->SetTitleW(0.30);
 604   gStyle->SetTitleH(0.10);
 605
 606   // --- Make canvas for correlation plots and plot them -------------
 607   TCanvas* c1 = new TCanvas("c1", "c1", 900, 600);
 608   c1->SetTopMargin(0.05);
 609   c1->SetRightMargin(0.05);
 610   c1->SetFillColor(0);
 611   c1->SetFillStyle(0);
 612   c1->SetBorderSize(0);
 613   c1->SetBorderMode(0);
 614
 615   c1->cd();
 616   TPad* p1 = new TPad("p1","p1", 0.0, 0.3, 1.0, 1.0, 0, 0, 0);
 617   p1->SetFillColor(0);
 618   p1->SetFillStyle(0);
 619   p1->SetBorderSize(0);
 620   p1->SetBorderMode(0);
 621   p1->SetTopMargin(0.05);
 622   p1->SetBottomMargin(0.001);
 623   p1->SetRightMargin(0.05);
 624   p1->Draw();
 625   p1->cd();
 626   clusters->SetMarkerStyle(20);
 627   clusters->SetMarkerColor(kRed+1);
 628   clusters->Draw();
 629
 630   tracklets->SetMarkerStyle(20);
 631   tracklets->SetMarkerColor(kBlue+1);
 632   tracklets->Draw("same");
 633
 634   TString v(Form("%+5.1f<|v_{z}|<%+5.1f",
 635                  events->GetXaxis()->GetXmin(),
 636                  events->GetXaxis()->GetXmax()));
 637   TLatex* ltx = new TLatex(.2,.80, v.Data());
 638   ltx->Draw();
 639
 640   TLegend* l = p1->BuildLegend(.6,.6,.94,.94);
 641   l->SetFillColor(0);
 642   l->SetFillStyle(0);
 643   l->SetBorderSize(0);
 644
 645   c1->cd();
 646   TPad* p2 = new TPad("p2","p2", 0.0, 0.0, 1.0, 0.3, 0, 0, 0);
 647   p2->SetFillColor(0);
 648   p2->SetFillStyle(0);
 649   p2->SetBorderSize(0);
 650   p2->SetBorderMode(0);
 651   p2->SetTopMargin(0.001);
 652   p2->SetBottomMargin(0.2);
 653   p2->SetRightMargin(0.05);
 654   p2->Draw();
 655   p2->cd();
 656   ratio->SetMarkerStyle(20);
 657   ratio->SetMarkerColor(kGray+1);
 658   ratio->Draw();
 659   ratio->GetYaxis()->SetRangeUser(0,2);
 660
 661   c1->SaveAs("SPDComparison.png");
 662   c1->SaveAs("SPDComparison.root");
 663 }
 664
 665
 666
 667
 668
 669 #endif
 670 //
 671 // EOF
 672 //
 673