PWGLF/FORWARD/analysis2/scripts/TupleSelector.C

   1 /**
   2  * @file   TupleSelector.C
   3  * @author Christian Holm Christensen <cholm@nbi.dk>
   4  * @date   Thu Nov 13 14:12:11 2014
   5  *
   6  * @brief A selector to draw stuff from the a track tuple as made by
   7  * AliFMDMCTrackELoss.
   8  *
   9  * This class assumes that the files with the TTree's created by
  10  * AliFMDMCTrackELoss lives in tuple/forward_tuple_N.root.
  11  *
  12  * The trees are made by AliFMDMCTrackELoss, which in turn is embedded
  13  * in a AliFMDMCTrackInspector object.  This code, which also does
  14  * fits of the MC 'true' energy loss spectra is contained in the task
  15  * AliFMDMCTrackInspectorTask.  This task can be executed via the
  16  * train MakeFMDMCTrackTrain.C
  17  *
  18  * To run this selector do
  19  *
  20  * @code
  21  void
  22  Run(Bool_t proof=true, Long64_t maxEvents=-1)
  23  {
  24    const char* fwd = "${ALICE_ROOT}/PWGLF/FORWARD/analysis2";
  25    gSystem->AddIncludePath("-I${ALICE_ROOT}/include");
  26    gROOT->Macro(Form("%s/scripts/LoadLibs.C"));
  27    gROOT->LoadMacro(Form("%s/TupleSelector.C++g",fwd));
  28
  29    if (proof) TupleSelector::Proof(maxEvents);
  30    else       TupleSelector::Run(maxEvents);
  31  }
  32  * @endcode
  33  *
  34  * Here, $ANA_SRC is assumed to point to the source directory of
  35  * PWGLF/FORWARD/analysis2
  36  */
  37
  38 #ifndef SELECTOR_C
  39 #define SELECTOR_C
  40
  41 #include <TSelector.h>
  42 #include <TQObject.h>
  43 #ifndef __CINT__
  44 # include <TH1.h>
  45 # include <TString.h>
  46 # include <TCanvas.h>
  47 # include <TLegend.h>
  48 # include <TLegendEntry.h>
  49 # include <TError.h>
  50 # include <TMath.h>
  51 # include <THStack.h>
  52 # include <TTree.h>
  53 # include <TClonesArray.h>
  54 # include <TChain.h>
  55 # include <TSystem.h>
  56 # include <TOutputListSelectorDataMap.h>
  57 # include <TLatex.h>
  58 # include <TROOT.h>
  59 # include <TFile.h>
  60 # include <TDirectory.h>
  61 # include <TSystemDirectory.h>
  62 # include <TRegexp.h>
  63 # include <TKey.h>
  64 # include <TFileCollection.h>
  65 # include <THashList.h>
  66 # include <TDSet.h>
  67 # include <TQConnection.h>
  68 # include <iostream>
  69 # include <TProof.h>
  70 # include <TStopwatch.h>
  71 # include "AliFMDMCTrackELoss.h"
  72 #else
  73 // Forward declarations of types in the interface, and to trigger
  74 // autoloading of some libraries
  75 class TTree;
  76 class TChain;
  77 class TCanvas;
  78 class TH1;
  79 class THStack;
  80 class TLegend;
  81 class TCanvas;
  82 class TString;
  83 class TDirectory;
  84 class TSystemDirectory;
  85 class TRegexp;
  86 class TDSet;
  87 class TProof;
  88 class AliFMDMCTrackELoss;
  89 class AliFMDMCTrackELoss::Hit;
  90 #endif
  91
  92 //====================================================================
  93 /**
  94  * Container of primary/secondary histograms
  95  */
  96 struct Spectra : public TObject
  97 {
  98   /** The name */
  99   TString fName;
 100   /** Primaries */
 101   TH1*    fPrimary;
 102   /** Secondaries */
 103   TH1*    fSecondary;
 104   /**
 105    * I/O CTOR
 106    */
 107   Spectra()
 108     : TObject(), fName(""), fPrimary(0), fSecondary(0)
 109   {}
 110   /**
 111    * User CTOR
 112    *
 113    * @param name   Name
 114    * @param title  Title
 115    * @param color  Color
 116    * @param bins   Bin definition
 117    */
 118   Spectra(const char*    name,
 119           const char*    title,
 120           Color_t        color,
 121           const TArrayD& bins)
 122     : TObject(),
 123       fName(name),
 124       fPrimary(0),
 125       fSecondary(0)
 126   {
 127     fPrimary = new TH1D(Form("primary%s", name),
 128                         title, bins.GetSize()-1,
 129                         bins.GetArray());
 130     fPrimary->SetXTitle("#beta#gamma");
 131     fPrimary->SetMarkerColor(color);
 132     fPrimary->SetLineColor(color);
 133     fPrimary->SetMarkerStyle(20);
 134     fPrimary->SetFillStyle(0);
 135     fPrimary->SetDirectory(0);
 136     fPrimary->Sumw2();
 137
 138     fSecondary = static_cast<TH1*>(fPrimary->Clone(Form("secondary%s",name)));
 139     fSecondary->SetMarkerStyle(24);
 140     fSecondary->SetMarkerSize(1.2);
 141     fSecondary->SetDirectory(0);
 142   }
 143   /**
 144    * Copy CTOR
 145    *
 146    * @param o Object to copy from
 147    */
 148   Spectra(const Spectra& o)
 149     : TObject(o),
 150       fName(o.fName),
 151       fPrimary(0),
 152       fSecondary(0)
 153   {
 154     if (o.fPrimary)
 155       fPrimary = static_cast<TH1*>(o.fPrimary->Clone());
 156     if (o.fSecondary)
 157       fSecondary = static_cast<TH1*>(o.fSecondary->Clone());
 158   }
 159   /**
 160    * DTOR
 161    */
 162   virtual ~Spectra()
 163   {
 164     if (fPrimary)   delete fPrimary;
 165     if (fSecondary) delete fSecondary;
 166   }
 167   /**
 168    * Assigment operator
 169    *
 170    * @param o Object to assign from
 171    *
 172    * @return Reference to this
 173    */
 174   Spectra& operator=(const Spectra& o)
 175   {
 176     if (&o == this) return *this;
 177
 178     TObject::operator=(o);
 179     fName = o.fName;
 180     if (fPrimary)   { delete fPrimary;   fPrimary   = 0; }
 181     if (fSecondary) { delete fSecondary; fSecondary = 0; }
 182     if (o.fPrimary)
 183       fPrimary = static_cast<TH1*>(o.fPrimary->Clone());
 184     if (o.fSecondary)
 185       fSecondary = static_cast<TH1*>(o.fSecondary->Clone());
 186
 187     return *this;
 188   }
 189   /**
 190    * Get the name
 191    *
 192    * @return name
 193    */
 194   const char* GetName() const { return fName.Data(); }
 195   /**
 196    * Fill into histogram
 197    *
 198    * @param primary true if for primary
 199    * @param x       What to fill
 200    */
 201   void Fill(Bool_t primary, Double_t x)
 202   {
 203     if (primary) fPrimary->Fill(x);
 204     else         fSecondary->Fill(x);
 205   }
 206   /**
 207    * Get a histogram
 208    *
 209    * @param primary If true, get histogram for primaries, otherwise
 210    * for secondaries.
 211    *
 212    * @return Pointer to histogram (possibly null)
 213    */
 214   TH1* Hist(Bool_t primary) const
 215   {
 216     return (primary ? fPrimary : fSecondary);
 217   }
 218   /**
 219    * Add the histograms to a stack and legend
 220    *
 221    * @param stack Stack to add to
 222    * @param leg   Legend
 223    */
 224   virtual void Stack(THStack* stack, TLegend* leg)
 225   {
 226     TH1* hists[] = { fPrimary, fSecondary, 0 };
 227     for (Int_t i = 0; i < 2; i++) {
 228       if (!hists[i]) continue;
 229
 230       Int_t n = hists[i]->GetEntries();
 231       if (n <= 0) continue;
 232
 233       if (stack) {
 234         hists[i]->Scale(1. / n, "width");
 235         stack->Add(hists[i]);
 236
 237         if (i != 0) continue;
 238         if (!leg) continue;
 239
 240         leg->AddEntry(hists[i], hists[i]->GetTitle(), "p");
 241       }
 242       else if (leg) {
 243         TLegendEntry* e = leg->AddEntry("",
 244                                         (i == 0 ? "Primary" : "Secondary"),
 245                                         "p");
 246         e->SetMarkerStyle(hists[i]->GetMarkerStyle());
 247         e->SetFillStyle(hists[i]->GetFillStyle());
 248         e->SetFillColor(kBlack);
 249       }
 250     }
 251   }
 252   /**
 253    * Merge this object with similar objects in the list.
 254    *
 255    * @param list list of objects
 256    *
 257    * @return Number of merged objects
 258    */
 259   Long64_t Merge(TCollection *list)
 260   {
 261     TIter    nxt(list);
 262     TObject* obj = 0;
 263     Long64_t cnt = 0;
 264     while ((obj = nxt())) {
 265       if (!obj->IsA()->InheritsFrom(this->IsA())) {
 266         Warning("Merge", "Will not merge a Spectra with a %s",
 267                 obj->IsA()->GetName());
 268         continue;
 269       }
 270       Spectra* oth  = static_cast<Spectra*>(obj);
 271       if (Add(oth)) cnt++;
 272     }
 273     // Info("Merge", "%s merged %lld entries", fName.Data(), cnt);
 274     return cnt;
 275   }
 276   Bool_t Add(const Spectra* oth)
 277   {
 278     if (!oth->fName.EqualTo(fName)) {
 279       Warning("Merge", "Will not merge %s with %s",
 280               oth->fName.Data(), fName.Data());
 281       return false;
 282     }
 283     // Int_t nPrmOld = fPrimary->GetEntries();
 284     // Int_t nSecOld = fPrimary->GetEntries();
 285
 286     fPrimary  ->Add(oth->fPrimary);
 287     fSecondary->Add(oth->fSecondary);
 288
 289     // Int_t nPrmNow = fPrimary->GetEntries();
 290     // Int_t nSecNow = fPrimary->GetEntries();
 291
 292     // Info("Add", "%s: %d/%d merged to %d/%d",
 293     //      fName.Data(), nPrmOld, nSecOld, nPrmNow, nSecNow);
 294     return true;
 295   }
 296   /**
 297    * List this object
 298    *
 299    * @param option Not used
 300    */
 301   void ls(Option_t* option="") const
 302   {
 303     TObject::ls(option);
 304     gROOT->IncreaseDirLevel();
 305     TH1* hists[] = { fPrimary,
 306                      fSecondary,
 307                      0 };
 308     for (Int_t i = 0; i < 2; i++) if (hists[i]) hists[i]->ls();
 309     gROOT->DecreaseDirLevel();
 310   }
 311   /**
 312    * Store this on file
 313    *
 314    * @param dir Parent directory
 315    */
 316   void Store(TDirectory* dir)
 317   {
 318     TDirectory* out = dir->mkdir(fName);
 319     out->cd();
 320     fPrimary->Clone("primary")->Write();
 321     fSecondary->Clone("secondary")->Write();
 322     dir->cd();
 323   }
 324   ClassDef(Spectra,1);
 325 };
 326 //====================================================================
 327 /**
 328  * Container to type histograms
 329  */
 330 struct Type : public Spectra
 331 {
 332   /**
 333    * I/O CTOR
 334    */
 335   Type() : Spectra() {}
 336   /**
 337    * User CTOR
 338    *
 339    * @param name   Name
 340    * @param title  Title
 341    * @param color  Color
 342    */
 343   Type(const char*    name,
 344        const char*    title,
 345        Color_t        color)
 346     : Spectra() // Do not use the base class ctor
 347   {
 348     fName = name;
 349     fPrimary = new TH1D(Form("primary%s", name),
 350                         title, 6, .5, 6.5);
 351     fPrimary->SetXTitle("particle type");
 352     fPrimary->SetMarkerColor(color);
 353     fPrimary->SetMarkerStyle(20);
 354     fPrimary->SetFillColor(color);
 355     fPrimary->SetFillStyle(3001);
 356     fPrimary->SetLineColor(color);
 357     fPrimary->SetDirectory(0);
 358     fPrimary->Sumw2();
 359     fPrimary->GetXaxis()->SetBinLabel(1, "e^{#pm}");
 360     fPrimary->GetXaxis()->SetBinLabel(2, "#pi^{#pm}");
 361     fPrimary->GetXaxis()->SetBinLabel(3, "K^{#pm}");
 362     fPrimary->GetXaxis()->SetBinLabel(4, "p/#bar{p}");
 363     fPrimary->GetXaxis()->SetBinLabel(5, "strange");
 364     fPrimary->GetXaxis()->SetBinLabel(6, "other");
 365
 366     fSecondary =
 367       static_cast<TH1*>(fPrimary->Clone(Form("secondary%s",name)));
 368     fSecondary->SetMarkerStyle(24);
 369     fSecondary->SetMarkerSize(1.2);
 370     fSecondary->SetFillStyle(3002);
 371     fSecondary->SetDirectory(0);
 372   }
 373   ClassDef(Type,1);
 374 };
 375
 376 //====================================================================
 377 struct Ring : public TObject
 378 {
 379
 380   /** Detector number */
 381   UShort_t fD;
 382   /** Ring identifier */
 383   Char_t   fR;
 384   /** Cached name */
 385   TString  fName;
 386   /** Collection of spectra */
 387   TObjArray* fSpectra;
 388
 389   /** Enumeration of our spectra */
 390   enum {
 391     kBetaGamma = 0,
 392     kBetaGammaPi,
 393     kBetaGammaElectron,
 394     kTypes,
 395     kNSpectra
 396   };
 397
 398   /**
 399    * I/O constructor
 400    */
 401   Ring()
 402     : TObject(),
 403       fD(0),
 404       fR('\0'),
 405       fName(),
 406       fSpectra(0)
 407   {}
 408   /**
 409    * Constructor
 410    *
 411    * @param d Detector
 412    * @param r Ring
 413    */
 414   Ring(UShort_t d, Char_t r)
 415     : fD(d),
 416       fR(r),
 417       fName(Form("FMD%d%c", d, r)),
 418       fSpectra(0)
 419   {
 420     Color_t color = AliForwardUtil::RingColor(fD, fR);
 421     TArrayD bgArray(601);
 422     AliForwardUtil::MakeLogScale(300, -2, 5, bgArray);
 423
 424     fSpectra   = new TObjArray(kNSpectra);
 425     fSpectra->SetOwner(true);
 426     fSpectra->AddAt(new Spectra("BetaGamma", fName, color, bgArray),
 427                     kBetaGamma);
 428     fSpectra->AddAt(new Spectra("BetaGammaPi", fName, color, bgArray),
 429                     kBetaGammaPi);
 430     fSpectra->AddAt(new Spectra("BetaGammaElectron", fName, color, bgArray),
 431                     kBetaGammaElectron);
 432     fSpectra->AddAt(new Type("Type", fName, color), kTypes);
 433     // fSpectra->ls();
 434   }
 435   /**
 436    * Copy constructor
 437    *
 438    * @param r Other to copy from
 439    */
 440   Ring(const Ring& r)
 441     : TObject(r),
 442       fD(r.fD),
 443       fR(r.fR),
 444       fName(r.fName),
 445       fSpectra(0)
 446   {
 447     if (r.fSpectra) fSpectra = new TObjArray(*r.fSpectra);
 448   }
 449   /**
 450    * Assignment operator
 451    *
 452    * @param r Other to assign from
 453    *
 454    * @return Reference to this obect
 455    */
 456   Ring& operator=(const Ring& r)
 457   {
 458     if (&r == this) return *this;
 459
 460     fD                     = r.fD;
 461     fR                     = r.fR;
 462     fName                  = r.fName;
 463     if (fSpectra) { delete fSpectra; fSpectra = 0; }
 464     if (r.fSpectra) fSpectra = new TObjArray(*r.fSpectra);
 465
 466     return *this;
 467   }
 468   const char* GetName() const { return fName.Data(); }
 469   /**
 470    * Get a spectra
 471    *
 472    * @param which Identifier
 473    *
 474    * @return Pointer to Spectra object or null
 475    */
 476   Spectra* Get(UInt_t which)
 477   {
 478     if (!fSpectra) return 0;
 479     Int_t i = which;
 480     if (i > fSpectra->GetEntriesFast()) return 0;
 481     return static_cast<Spectra*>(fSpectra->At(i));
 482   }
 483   /**
 484    * Get a spectra
 485    *
 486    * @param which Identifier
 487    *
 488    * @return Pointer to Spectra object or null
 489    */
 490   const Spectra* Get(UInt_t which) const
 491   {
 492     if (!fSpectra) return 0;
 493     Int_t i = which;
 494     if (i > fSpectra->GetEntriesFast()) return 0;
 495     return static_cast<Spectra*>(fSpectra->At(i));
 496   }
 497   /**
 498    * Fill histograms
 499    *
 500    * @param hit Hit structure
 501    */
 502   void Fill(AliFMDMCTrackELoss::Hit* hit)
 503   {
 504     Bool_t prim = hit->IsPrimary();
 505     Int_t  apdg = hit->AbsPdg();
 506     Int_t  mfl  = Int_t(apdg/TMath::Power(10,Int_t(TMath::Log10(apdg))));
 507     Int_t  type = (hit->IsElectron() ? 1 :
 508                    hit->IsPion()     ? 2 :
 509                    hit->IsKaon()     ? 3 :
 510                    hit->IsProton()   ? 4 :
 511                    mfl == 3          ? 5 : 6);
 512
 513     Get(kBetaGamma)        ->Fill(prim, hit->BetaGamma());
 514     Get(kTypes)            ->Fill(prim, type);
 515     if (hit->IsPion() || hit->IsProton() || hit->IsKaon() || apdg > 100)
 516       Get(kBetaGammaPi)      ->Fill(prim, hit->BetaGamma());
 517     if (hit->IsElectron())
 518       Get(kBetaGammaElectron)->Fill(prim, hit->BetaGamma());
 519
 520   }
 521   /**
 522    * Get histograms
 523    *
 524    * @param primary If true, for primaries
 525    * @param which   Which histogram to get
 526    *
 527    * @return Histogram or null
 528    */
 529   TH1* Hist(Bool_t primary, UShort_t which) const
 530   {
 531     const Spectra* spe = Get(which);
 532     if (!spe) return 0;
 533     return spe->Hist(primary);
 534   }
 535   /**
 536    * Merge this object with similar objects in the list.
 537    *
 538    * @param list list of objects
 539    *
 540    * @return Number of merged objects
 541    */
 542   Long64_t Merge(TCollection *list)
 543   {
 544     TIter    nxt(list);
 545     TObject* obj = 0;
 546     Long64_t cnt = 0;
 547     while ((obj = nxt())) {
 548       if (!obj->IsA()->InheritsFrom(this->IsA())) {
 549         Warning("Merge", "Will not merge a Ring with a %s",
 550                 obj->IsA()->GetName());
 551         continue;
 552       }
 553       Ring* oth = static_cast<Ring*>(obj);
 554       if (oth->fD != fD || oth->fR != fR) {
 555         Warning("Merge", "Will not merge FMD%d%c with FMD%d%c",
 556                 fD, fR, oth->fD, oth->fR);
 557         continue;
 558       }
 559
 560       if (fSpectra) {
 561         // fSpectra->Merge(oth->fSpectra);
 562         TIter thsNext(fSpectra);
 563         TIter othNext(oth->fSpectra);
 564         Spectra* thsSpec = 0;
 565         Spectra* othSpec = 0;
 566
 567         while ((thsSpec = static_cast<Spectra*>(thsNext())) &&
 568                (othSpec = static_cast<Spectra*>(othNext())))
 569           thsSpec->Add(othSpec);
 570
 571       }
 572       cnt++;
 573     }
 574     // Info("Merge", "FMD%d%c merged %lld entries", fD, fR, cnt);
 575     return cnt;
 576   }
 577   void Draw(Option_t* opt="")
 578   {
 579     Info("Draw", "Should draw rings here");
 580     ls(opt);
 581   }
 582   /**
 583    * List this object
 584    *
 585    * @param option Not used
 586    */
 587   void ls(Option_t* option="") const
 588   {
 589     TObject::ls(option);
 590     gROOT->IncreaseDirLevel();
 591     if (fSpectra) fSpectra->ls(option);
 592     gROOT->DecreaseDirLevel();
 593   }
 594   /**
 595    * Store this on file
 596    *
 597    * @param dir Parent directory
 598    */
 599   void Store(TDirectory* dir)
 600   {
 601     TDirectory* out = dir->mkdir(fName);
 602     out->cd();
 603     TIter next(fSpectra);
 604     Spectra* spec = 0;
 605     while ((spec = static_cast<Spectra*>(next())))
 606       spec->Store(out);
 607     dir->cd();
 608   }
 609
 610   ClassDef(Ring,2);
 611 };
 612
 613 //====================================================================
 614 struct Monitor : public TObject, public TQObject
 615 {
 616   Monitor()
 617   {
 618     if (!gProof) return;
 619
 620     fName = gProof->GetSessionTag();
 621     gDirectory->Add(this);
 622     // _must_ specify signal _exactly_ like below, or we won't be called
 623     Bool_t ret = gProof->Connect("Feedback(TList *objs)", "Monitor", this,
 624                                  "Feedback(TList *objs)");
 625     if (!ret)
 626       Warning("Monitor", "Failed to connect to Proof");
 627   }
 628   virtual ~Monitor()
 629   {
 630     if (!gProof) return;
 631     gProof->Disconnect("Feedback(TList *objs)",this,
 632                        "Feedback(TList* objs)");
 633   }
 634   void SetName(const char* name) { fName = name; }
 635   const char* GetName() const { return fName.Data(); }
 636   void Feedback(TList* objs)
 637   {
 638     Info("Feedback", "Got a list of objects (%p)", objs);
 639     if (!objs) return;
 640     objs->ls();
 641   }
 642   TString fName;
 643   ClassDef(Monitor,1);
 644 };
 645
 646
 647 //====================================================================
 648 struct TupleSelector : public TSelector
 649 {
 650   TString       fTitle; //! Must not be persistent
 651   TTree*        fTree;  //! Must not be persistent
 652   TClonesArray* fHits;  //! Must not be persistent
 653   Int_t         fI;     //! Not persistent
 654   TObjArray*    fRings; //! Not persistent
 655
 656   /**
 657    * Constructor
 658    *
 659    * @param name Optional title
 660    */
 661   TupleSelector(const char* name="") :
 662     fTitle(name), fTree(0), fHits(0), fI(0), fRings(0)
 663   {
 664     if (fTitle.IsNull())
 665       fTitle = gSystem->BaseName(gSystem->WorkingDirectory());
 666     // SetOption("fb=rings");
 667
 668   }
 669   const char* GetName() const { return fTitle.Data(); }
 670   /**
 671    * @{
 672    * @name Setup
 673    */
 674   /**
 675    * Get index for a particular ring
 676    *
 677    * @param d Detector
 678    * @param r Ring
 679    *
 680    * @return Index, or 0xFFFF
 681    */
 682   UShort_t Index(UShort_t d, Char_t r) const
 683   {
 684     Bool_t inner = (r == 'i' || r == 'I');
 685     switch (d) {
 686     case 1: return 0;
 687     case 2: return (inner ? 1 : 2);
 688     case 3: return (inner ? 4 : 3);
 689     }
 690     ::Warning("", "Unknown ring FMD%d%c", d, r);
 691     return 0xFFFF;
 692   }
 693   void Init(TTree* tree)
 694   {
 695     Info("Init", "Got a tree: %p", tree);
 696
 697     if (!tree) {
 698       Warning("Init", "No tree passed");
 699       return;
 700     }
 701
 702     fTree = tree;
 703     fHits = new TClonesArray("AliFMDMCTrackELoss::Hit");
 704     fTree->SetBranchAddress("hits", &fHits);
 705   }
 706   void Begin(TTree* tree)
 707   {
 708     Info("Begin", "Called w/tree @ %p", tree);
 709     // if (tree) SlaveBegin(tree);
 710     new Monitor;
 711   }
 712   /**
 713    * Begin on slave
 714    *
 715    * @param tree Tree to process
 716    */
 717   void SlaveBegin(TTree* tree)
 718   {
 719     Info("SlaveBegin", "Got a tree: %p", tree);
 720     fRings = new TObjArray(5);
 721     fRings->SetName("rings");
 722     fRings->SetOwner(true);
 723     fRings->AddAt(new Ring(1,'I'), Index(1,'I'));
 724     fRings->AddAt(new Ring(2,'I'), Index(2,'I'));
 725     fRings->AddAt(new Ring(2,'O'), Index(2,'O'));
 726     fRings->AddAt(new Ring(3,'O'), Index(3,'O'));
 727     fRings->AddAt(new Ring(3,'I'), Index(3,'I'));
 728
 729     fOutput->Add(fRings);
 730     if (tree) {
 731       // In case of local mode, we get the tree here,
 732       // and init isn't called
 733       Init(tree);
 734     }
 735   }
 736   /* @} */
 737
 738   /**
 739    * @{
 740    * @name Event processing
 741    */
 742   Bool_t Notify()
 743   {
 744     TFile* file = (fTree ? fTree->GetCurrentFile() : 0);
 745     Info("Notify","processing file: %p (%s)",
 746          file, (file ? file->GetName() : "nil"));
 747     return true;
 748   }
 749   /**
 750    * Process an entry
 751    *
 752    * @param entry Entry
 753    *
 754    * @return true on success
 755    */
 756   Bool_t Process(Long64_t entry)
 757   {
 758     if (!fTree) {
 759       Warning("Process", "No tree");
 760       return false;
 761     }
 762     if (!fTree->GetTree()) {
 763       Warning("Process", "No real tree");
 764       return false;
 765     }
 766     fHits->Clear();
 767     fTree->GetTree()->GetEntry(entry);
 768     fI++;
 769     // if (fI % 100 == 0) {
 770     //   printf("Event # %6d (%6d hits)\r", fI, fHits->GetEntries());
 771     //   fflush(stdout);
 772     // }
 773     // Printf("Event # %7d, %6d hits", fI, fHits->GetEntries());
 774
 775     TIter next(fHits);
 776     AliFMDMCTrackELoss::Hit* hit = 0;
 777     while ((hit = static_cast<AliFMDMCTrackELoss::Hit*>(next())))
 778       Fill(hit);
 779    return true;
 780   }
 781   /**
 782    * Get a Ring
 783    *
 784    * @param d Detector
 785    * @param r Ring
 786    *
 787    * @return Pointer to ring object or null
 788    */
 789   Ring* Get(UShort_t d, Char_t r)
 790   {
 791     return static_cast<Ring*>(fRings->At(Index(d,r)));
 792   }
 793   /**
 794    * Get a Ring
 795    *
 796    * @param d Detector
 797    * @param r Ring
 798    *
 799    * @return Pointer to constant ring object or null
 800    */
 801   const Ring* Get(UShort_t d, Char_t r) const
 802   {
 803     return static_cast<Ring*>(fRings->At(Index(d,r)));
 804   }
 805   /**
 806    * Fill in to a ring object
 807    *
 808    * @param hit Hit information
 809    */
 810   void Fill(AliFMDMCTrackELoss::Hit* hit)
 811   {
 812     Ring* r = Get(hit->D(), hit->R());
 813     r->Fill(hit);
 814   }
 815   /* @} */
 816
 817   /**
 818    * @{
 819    * @name Final processing
 820    */
 821   void SlaveTerminate()
 822   {
 823   }
 824   /**
 825    * Terminate the job
 826    *
 827    */
 828   void Terminate()
 829   {
 830     Printf("\nDone (rings=%p)", fRings);
 831
 832     if (!fRings)
 833       fRings = static_cast<TObjArray*>(fOutput->FindObject("rings"));
 834     if (!fRings) {
 835       Error("Terminate", "No rings in output array");
 836       return;
 837     }
 838     // fRings->ls();
 839
 840     TFile* out = TFile::Open("tuple_summary.root", "RECREATE");
 841     Store(out);
 842
 843     PlotOne(Ring::kBetaGamma,        "#beta#gamma - all",     out);
 844     PlotOne(Ring::kBetaGammaPi,      "#beta#gamma - Hadrons", out);
 845     PlotOne(Ring::kBetaGammaElectron,"#beta#gamma - e^{#pm}", out);
 846     PlotOne(Ring::kTypes,            "Relative abundance",    out);
 847   }
 848   /**
 849    * Plot one quantity
 850    *
 851    * @param which
 852    * @param what
 853    */
 854   void PlotOne(UShort_t which, const char* title,
 855                TDirectory* dir=0, Option_t* opt="")
 856   {
 857     static Int_t cnt  = 1;
 858     const char*  name = Form("c%02d", cnt++);
 859
 860     TCanvas* c = new TCanvas(name, title, 1000, 1000);
 861     c->SetFillColor(0);
 862     c->SetFillStyle(0);
 863     c->SetTopMargin(0.01);
 864     c->SetRightMargin(0.03);
 865
 866     TLegend* leg   = new TLegend(0.65, 0.65, 0.975, 0.975);
 867     leg->SetFillStyle(0);
 868     leg->SetBorderSize(0);
 869
 870     TString xTitle = "";
 871     THStack* stack = new THStack(name, title);
 872     for (Int_t i = 0; i < 5; i++) {
 873       Ring*    ring = static_cast<Ring*>(fRings->At(i));
 874       Spectra* spec = ring->Get(which);
 875       if (!spec) {
 876         Warning("PlotOne", "No spectra for %d", which);
 877         continue;
 878       }
 879
 880       // Add to our stack
 881       spec->Stack(stack, leg);
 882
 883       // Get X title
 884       if (xTitle.IsNull() && spec->Hist(true))
 885         xTitle = spec->Hist(true)->GetXaxis()->GetTitle();
 886
 887       // If we're not at the last one continue
 888       if (i < 4) continue;
 889
 890       // Add primary/secondary entry to legend
 891       spec->Stack(0, leg);
 892
 893     }
 894     c->cd();
 895     if (!stack->GetHists() ||
 896         stack->GetHists()->GetEntries() < 0) {
 897       Warning("PlotOne", "No histograms added");
 898       return;
 899     }
 900     stack->Draw(Form("nostack %s", opt));
 901     stack->GetHistogram()->SetXTitle(xTitle);
 902     stack->GetHistogram()->GetListOfFunctions()->Add(leg);
 903
 904     leg->Draw();
 905
 906     if (which != Ring::kTypes) {
 907       c->SetLogx();
 908       c->SetLogy();
 909     }
 910
 911     TLatex* ltx = new TLatex(0.02, 0.02, fTitle);
 912     ltx->SetNDC();
 913     ltx->SetTextColor(kRed+2);
 914     ltx->SetTextAlign(11);
 915     ltx->SetTextSize(0.02);
 916     ltx->Draw();
 917     stack->GetHistogram()->GetListOfFunctions()->Add(ltx);
 918
 919     if (dir) {
 920       dir->cd();
 921       stack->Clone()->Write();
 922     }
 923
 924
 925     c->Modified();
 926     c->Update();
 927     c->cd();
 928
 929     c->Print(Form("%s.pdf", name));
 930   }
 931   /**
 932    * Store this on file
 933    *
 934    * @param dir Parent directory
 935    */
 936   void Store(TDirectory* dir)
 937   {
 938     TDirectory* out = dir->mkdir(fTitle);
 939     out->cd();
 940     TIter next(fRings);
 941     Ring* ring = 0;
 942     while ((ring = static_cast<Ring*>(next())))
 943       ring->Store(out);
 944     dir->cd();
 945   }
 946   /* @} */
 947
 948   /**
 949    * @{
 950    * @name Other interface
 951    */
 952   /**
 953    * Get the version of the selector
 954    *
 955    * @return 1
 956    */
 957   Int_t Version() const { return 1; }
 958   /**
 959    * Get the status
 960    *
 961    * @return Number of processedd events
 962    */
 963   Long64_t GetStatus() const { return fI; }
 964
 965   /* @} */
 966
 967   //------------------------------------------------------------------
 968   /**
 969    * @{
 970    * @name Service functions
 971    */
 972   //------------------------------------------------------------------
 973   /**
 974    * Get a chain or data set from the given file
 975    *
 976    * @param file File to look in
 977    *
 978    * @return Pointer to object or null
 979    */
 980   static TObject* GetChainOrDataSet(const char* file)
 981   {
 982     if (gSystem->AccessPathName(file)) return 0;
 983
 984     TFile* in = TFile::Open(file, "READ");
 985     if (!in) return 0;
 986
 987     TObject* ret = in->Get("tree");
 988     if (!ret) {
 989       in->Close();
 990       return 0;
 991     }
 992
 993     if (ret->IsA()->InheritsFrom(TChain::Class()))
 994       static_cast<TChain*>(ret)->SetDirectory(0);
 995     else if (ret->IsA()->InheritsFrom(TDSet::Class()))
 996       static_cast<TDSet*>(ret)->SetDirectory(0);
 997     else {
 998       ::Warning("GetChainOrDataSet", "Found object is a %s",
 999                 ret->IsA()->GetName());
1000       ret = 0;
1001     }
1002     in->Close();
1003     return ret;
1004   }
1005   //------------------------------------------------------------------
1006   /**
1007    * make our data set
1008    *
1009    * @param src        Source directory
1010    * @param recursive  Whether to scan recursively
1011    * @param verbose    Be verbose
1012    *
1013    * @return Data set or null
1014    */
1015   static TDSet* MakeDataSet(const TString& src=".",
1016                             Bool_t recursive=false,
1017                             Bool_t verbose=false)
1018   {
1019     TString dsFile(Form("%s/dataset.root", src.Data()));
1020     TDSet* dataset = static_cast<TDSet*>(GetChainOrDataSet(dsFile));
1021     if (dataset) {
1022       /// dataset->Print("a");
1023       return dataset;
1024     }
1025
1026     TChain* c = DoMakeChain(src, recursive, verbose);
1027     if (!c) return 0;
1028
1029     dataset = new TDSet(*c, false);
1030     dataset->SetName("tree");
1031     dataset->SetLookedUp();
1032     dataset->Validate();
1033
1034     delete c;
1035     if (dataset) {
1036       TFile* out = TFile::Open(dsFile, "RECREATE");
1037       dataset->Write();
1038       dataset->SetDirectory(0);
1039       out->Close();
1040     }
1041
1042     return dataset;
1043
1044   }
1045   //------------------------------------------------------------------
1046   /**
1047    * Create our chain
1048    *
1049    * @param src        Source directory
1050    * @param recursive  Whether to scan recursively
1051    * @param verbose    Be verbose
1052    *
1053    * @return Chain or null
1054    */
1055   static TChain* MakeChain(const TString& src=".",
1056                            Bool_t recursive=false,
1057                            Bool_t verbose=false)
1058   {
1059     TString chFile(Form("%s/chain.root", src.Data()));
1060     if (!gSystem->AccessPathName(chFile)) {
1061       TFile* in = TFile::Open(chFile, "READ");
1062       if (in) {
1063         TChain* ret = static_cast<TChain*>(in->Get("tree"));
1064         if (ret) {
1065           ret->SetDirectory(0);
1066           // ret->GetListOfFiles()->ls();
1067         }
1068         in->Close();
1069         if (ret) return ret;
1070       }
1071     }
1072
1073     TChain* chain = DoMakeChain(src, recursive, verbose);
1074     if (chain) {
1075       TFile* out = TFile::Open(chFile, "RECREATE");
1076       chain->Write();
1077       chain->SetDirectory(0);
1078       out->Close();
1079     }
1080
1081     return chain;
1082   }
1083   //------------------------------------------------------------------
1084   /**
1085    * Create our chain
1086    *
1087    * @param src        Source directory
1088    * @param recursive  Whether to scan recursively
1089    * @param verbose    Be verbose
1090    *
1091    * @return Chain or null
1092    */
1093   static TChain* DoMakeChain(const TString& src=".",
1094                              Bool_t recursive=false,
1095                              Bool_t verbose=false)
1096   {
1097     TChain* chain = new TChain("tree");
1098     TString savdir(gSystem->WorkingDirectory());
1099     TSystemDirectory d(gSystem->BaseName(src.Data()), src.Data());
1100     if (!ScanDirectory(&d, chain, "forward_tree_*", recursive, verbose)) {
1101       delete chain;
1102       chain = 0;
1103     }
1104     else if (!chain->GetListOfFiles() ||
1105         chain->GetListOfFiles()->GetEntries() <= 0) {
1106       delete chain;
1107       chain = 0;
1108     }
1109     return chain;
1110   }
1111
1112   //------------------------------------------------------------------
1113   /**
1114    * Scan directory @a dir (possibly recursive) for tree files to add
1115    * to the chain.    This does not follow sym-links
1116    *
1117    * @param dir        Directory to scan
1118    * @param chain      Chain to add to
1119    * @param pattern    File name pattern
1120    * @param anchor     Anchor (tree name)
1121    * @param flags      Flags
1122    *
1123    * @return true if any files where added
1124    */
1125   static Bool_t ScanDirectory(TSystemDirectory* dir,
1126                               TChain*           chain,
1127                               const TString&    pattern,
1128                               Bool_t            recursive,
1129                               Bool_t            verbose)
1130   {
1131     if (!dir) {
1132       ::Error("ScanDirectory", "No diretory passed");
1133       return false;
1134     }
1135     if (verbose) ::Info("ScanDirectory", "Scanning %s", dir->GetName());
1136
1137     Bool_t  ret = false;
1138     TRegexp wild(pattern, true);
1139     TString oldDir(gSystem->WorkingDirectory());
1140     TList* files = dir->GetListOfFiles();
1141
1142     if (!gSystem->ChangeDirectory(oldDir)) {
1143       ::Error("ScanDirectory", "Failed to go back to %s",
1144               oldDir.Data());
1145       return false;
1146     }
1147     if (!files) {
1148       ::Warning("ScanDirectory", "No files");
1149       return false;
1150     }
1151
1152     TList toAdd;
1153     toAdd.SetOwner();
1154
1155     // Sort list of files and check if we should add it
1156     files->Sort();
1157     TIter next(files);
1158     TSystemFile* file = 0;
1159     while ((file = static_cast<TSystemFile*>(next()))) {
1160       TString name(file->GetName());
1161       TString title(file->GetTitle());
1162       TString full(gSystem->ConcatFileName(file->GetTitle(), name.Data()));
1163       if (file->IsA()->InheritsFrom(TSystemDirectory::Class())) full = title;
1164
1165
1166       if (name.EqualTo(".")||name.EqualTo(".."))  {
1167         // Ignore special links
1168         if (verbose) ::Info("ScanDirectory", "Ignoring special %s",
1169                             name.Data());
1170         continue;
1171       }
1172
1173       if (verbose) ::Info("ScanDirectory", "Looking at %s", full.Data());
1174
1175       FileStat_t fs;
1176       if (gSystem->GetPathInfo(full.Data(), fs)) {
1177         ::Warning("ScanDirectory", "Cannot stat %s (%s)",
1178                   full.Data(), gSystem->WorkingDirectory());
1179         continue;
1180       }
1181       // Check if this is a directory
1182       if (file->IsDirectory(full)) {
1183         if (verbose) ::Info("ScanDirectory", "Got a directory: %s",
1184                             full.Data());
1185         if (recursive) {
1186           // if (title[0] == '/')
1187           TSystemDirectory* d = new TSystemDirectory(file->GetName(),
1188                                                      full.Data());
1189           if (ScanDirectory(d, chain, pattern, recursive, verbose))
1190             ret = true;
1191           delete d;
1192         }
1193         continue;
1194       }
1195
1196       // If this is not a root file, ignore
1197       if (!name.EndsWith(".root") &&
1198           !name.EndsWith(".zip",TString::kIgnoreCase)) {
1199         if (verbose) ::Info("ScanDirectory", "Ignoring non-ROOT or ZIP %s",
1200                             name.Data());
1201         continue;
1202       }
1203
1204       // If this file does not contain AliESDs, ignore
1205       if (!name.Contains(wild)) {
1206         if (verbose) ::Info("ScanDirectory", "%s does not match %s",
1207                             name.Data(), pattern.Data());
1208         continue;
1209       }
1210       // Add
1211       if (verbose)
1212         ::Info("::ScanDirectory", "Candidate %s", full.Data());
1213       toAdd.Add(new TObjString(full));
1214     }
1215     TIter nextAdd(&toAdd);
1216     TObjString* s = 0;
1217     Int_t added = 0;
1218     while ((s = static_cast<TObjString*>(nextAdd()))) {
1219       // Info("ChainBuilder::ScanDirectory",
1220       //      "Adding %s", s->GetString().Data());
1221       TString fn = s->GetString();
1222       if (!CheckFile(fn, chain, true/*verbose*/)) continue;
1223
1224       added++;
1225     }
1226     if (added > 0) ret = true;
1227     if (verbose) ::Info("ScanDirectory", "Added %d files", added);
1228
1229     gSystem->ChangeDirectory(oldDir);
1230     return ret;
1231   }
1232   //------------------------------------------------------------------
1233   /**
1234    * Check if we can add a file to the chain
1235    *
1236    * @param path   Full path to file
1237    * @param chain  Chain
1238    *
1239    * @return true on success, false otherwise
1240    */
1241   static Bool_t CheckFile(const TString& path,
1242                           TChain* chain,
1243                           Bool_t verbose)
1244   {
1245     // Info("", "Checking %s", path.Data());
1246     TString fn   = path;
1247
1248     gSystem->RedirectOutput("/dev/null", "w");
1249     TFile*  test = TFile::Open(fn, "READ");
1250     gSystem->RedirectOutput(0);
1251     if (!test) {
1252       ::Warning("CheckFile", "Failed to open %s", fn.Data());
1253       return false;
1254     }
1255
1256     Bool_t           ok = false;
1257     TObject*         o  = test->Get(chain->GetName());
1258     TTree*           t  = dynamic_cast<TTree*>(o);
1259     TFileCollection* c  = dynamic_cast<TFileCollection*>(o);
1260     if (t) {
1261       Int_t n = t->GetEntries();
1262       test->Close();
1263       chain->Add(fn, n);
1264       ok = true;
1265       if (verbose) ::Info("CheckFile", "Added file %s (%d)", fn.Data(), n);
1266     } else if (c) {
1267       chain->AddFileInfoList(c->GetList());
1268       ok = true;
1269       if (verbose) ::Info("CheckFile", "Added file collection %s", fn.Data());
1270     } else {
1271       // Let's try to find a TFileCollection
1272       TList* l = test->GetListOfKeys();
1273       TIter next(l);
1274       TKey* k = 0;
1275       while ((k = static_cast<TKey*>(next()))) {
1276         TString cl(k->GetClassName());
1277         if (!cl.EqualTo("TFileCollection")) continue;
1278         c = dynamic_cast<TFileCollection*>(k->ReadObj());
1279         if (!c) {
1280           ::Warning("CheckFile", "Returned collection invalid");
1281           continue;
1282         }
1283         // Info("", "Adding file collection");
1284         chain->AddFileInfoList(c->GetList());
1285         ok = true;
1286         if (verbose) ::Info("CheckFile", "Added file collection %s", fn.Data());
1287       }
1288       test->Close();
1289     }
1290
1291     if (!ok)
1292       ::Warning("CheckFile",
1293                 "The file %s does not contain the tree %s "
1294                 "or a file collection",
1295                 path.Data(), chain->GetName());
1296     return ok;
1297   }
1298
1299   /**
1300    * Run this selector on a chain locally
1301    *
1302    * @param maxEvents Maximum number of events
1303    * @param title     Optional title
1304    * @param maxFiles  Maximum number of files to put in chain
1305    *
1306    * @return true on sucess
1307    */
1308   static Bool_t Run(Long64_t maxEvents=-1,
1309                     const char* title="")
1310   {
1311     TStopwatch timer;
1312     timer.Start();
1313     TChain* chain = MakeChain("tuple");
1314     if (!chain) {
1315       ::Error("Run", "No chain!");
1316       return false;
1317     }
1318
1319     TupleSelector* s = new TupleSelector(title);
1320     Int_t status= chain->Process(s, "", maxEvents);
1321     timer.Print();
1322     return status >= 0;
1323   }
1324   /**
1325    * Run this selector on a chain in Proof
1326    *
1327    * @param maxEvents Maximum number of events
1328    * @param title     Optional title
1329    * @param maxFiles  Maximum number of files to put in chain
1330    *
1331    * @return true on sucess
1332    */
1333   static Bool_t Proof(Long64_t    maxEvents,
1334                       const char* opt="",
1335                       const char* title="")
1336   {
1337     TStopwatch timer;
1338     timer.Start();
1339     TProof::Reset("lite:///?workers=8");
1340     TProof::Open("lite:///?workers=8");
1341     gProof->ClearCache();
1342     TString ali = gSystem->ExpandPathName("$(ALICE_ROOT)");
1343     TString fwd = ali + "/PWGLF/FORWARD/analysis2";
1344     gProof->AddIncludePath(Form("%s/include", ali.Data()));
1345     gProof->Load(Form("%s/scripts/LoadLibs.C",fwd.Data()), true);
1346     gProof->Exec("LoadLibs()");
1347     // gROOT->ProcessLine("gProof->SetLogLevel(5);");
1348     gProof->Load(Form("%s/scripts/TupleSelector.C+%s", fwd.Data(), opt),true);
1349
1350     TDSet* dataset = MakeDataSet("tuple");
1351     if (!dataset) {
1352       ::Error("Proof", "No dataset");
1353       return false;
1354     }
1355
1356     TupleSelector* s = new TupleSelector(title);
1357     gProof->AddFeedback("rings");
1358     gProof->Process(dataset, s, "", maxEvents);
1359     timer.Print();
1360     return true; // status >= 0;
1361   }
1362
1363   /* @} */
1364   ClassDef(TupleSelector,2);
1365
1366 private:
1367   /**
1368    * Copy constructor
1369    */
1370   TupleSelector(const TupleSelector&); // {}
1371   /**
1372    * Assignment operator
1373    *
1374    * @return Reference to this object
1375    */
1376   TupleSelector& operator=(const TupleSelector&); // { return *this; }
1377
1378 };
1379
1380
1381
1382 #endif
1383 // Local Variables:
1384 //   mode: C++
1385 // End:
1386 //
1387 // EOF
1388 //