PWG2/FORWARD/analysis2/SPDComparison.C

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