ITS/UPGRADE/macros/QA/AliITSUComparisonPileup.C

   1 /****************************************************************************
   2  *  This macro calculates the efficiency of pileup reconstruction.          *
   3  *  Works only for events generated with the AliGenPileup generator.        *
   4  *                                                                          *
   5  *  Before running, load the ITSU libraries:                                *
   6  *  gSystem->Load("libITSUpgradeBase");gSystem->Load("libITSUpgradeRec");   *
   7  *                                                                          *
   8  *           Origin: I.Belikov, IPHC, Iouri.Belikov@iphc.cnrs.fr            *
   9  ****************************************************************************/
  10
  11 #if !defined(__CINT__) || defined(__MAKECINT__)
  12   #include <Riostream.h>
  13   #include <TMath.h>
  14   #include <TTree.h>
  15   #include <TParticle.h>
  16   #include <TCanvas.h>
  17   #include <TFile.h>
  18   #include <TROOT.h>
  19
  20   #include "AliStack.h"
  21   #include "AliHeader.h"
  22   #include "AliGenCocktailEventHeader.h"
  23   #include "AliRunLoader.h"
  24   #include "AliRun.h"
  25   #include "AliESDEvent.h"
  26   #include "AliESDtrack.h"
  27 #endif
  28
  29 Int_t GoodPileupVertices(const Char_t *dir=".");
  30
  31 extern AliRun *gAlice;
  32 extern TROOT *gROOT;
  33
  34 Int_t AliITSUComparisonPileup(const Char_t *dir=".") {
  35    ::Info("AliITSUComparisonPileup.C","Doing comparison...");
  36
  37    // **** Book histogramms
  38    TH2F *h2=(TH2F*)gROOT->FindObject("h2");
  39    if (!h2) h2=new TH2F("h2",";Number of good vertices;Number of reconstructed vertices",100,0.,100.,10,0.,10.);
  40
  41
  42    // **** Generate a rerefence file with reconstructable vertices
  43    Char_t fname[100];
  44    sprintf(fname,"%s/GoodPileupVertices.root",dir);
  45    TFile *refFile=TFile::Open(fname,"old");
  46    if (!refFile || !refFile->IsOpen()) {
  47       ::Info("AliITSUComparisonPileup.C",
  48       "Marking good pileup vertices (will take a while)...");
  49       if (GoodPileupVertices(dir)) {
  50      ::Error("AliITSUComparisonPileup.C","Can't generate the reference file !");
  51          return 1;
  52       }
  53    }
  54    refFile=TFile::Open(fname,"old");
  55    if (!refFile || !refFile->IsOpen()) {
  56      ::Error("AliITSUComparisonPileup.C","Can't open the reference file !");
  57      return 1;
  58    }
  59    TTree *refTree=(TTree*)refFile->Get("refTree");
  60    if (!refTree) {
  61      ::Error("AliITSUComparisonPileup.C","Can't get the reference tree !");
  62      return 2;
  63    }
  64    TBranch *branch=refTree->GetBranch("Vertices");
  65    if (!branch) {
  66      ::Error("AliITSUComparisonPileup.C","Can't get the vertex branch !");
  67      return 3;
  68    }
  69    TClonesArray dummy("AliESDVertex",100), *refs=&dummy;
  70    branch->SetAddress(&refs);
  71
  72
  73    // **** Open the ESD
  74    sprintf(fname,"%s/AliESDs.root",dir);
  75    TFile *ef=TFile::Open(fname);
  76    if ((!ef)||(!ef->IsOpen())) {
  77       ::Error("AliITSUComparisonPileup.C","Can't open AliESDs.root !");
  78       return 4;
  79    }
  80    AliESDEvent* event = new AliESDEvent();
  81    TTree* esdTree = (TTree*) ef->Get("esdTree");
  82    if (!esdTree) {
  83       ::Error("AliITSUComparisonPileup.C", "no ESD tree found");
  84       return 6;
  85    }
  86    event->ReadFromTree(esdTree);
  87
  88
  89    //******* Loop over events *********
  90    Int_t e=0;
  91    while (esdTree->GetEvent(e)) {
  92      cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
  93      TClonesArray *vertices=event->GetPileupVerticesTracks();
  94      Int_t nrec=vertices->GetEntriesFast();
  95
  96      if (refTree->GetEvent(e++)==0) {
  97         cerr<<"No reconstructable vertices for this event !\n";
  98         continue;
  99      }
 100      Int_t ngood=refs->GetEntriesFast();
 101      cout<<"reconstructed: "<<nrec<<"  good: "<<ngood<<endl;
 102
 103      h2->Fill(ngood,nrec);
 104
 105    } //***** End of the loop over events
 106
 107    delete event;
 108    delete esdTree;
 109    ef->Close();
 110
 111    refFile->Close();
 112
 113    TCanvas *c1=new TCanvas("c1","",0,0,750,500);
 114    h2->Draw("box");
 115
 116    TFile fc("AliITSUComparisonCooked.root","RECREATE");
 117    c1->Write();
 118    fc.Close();
 119
 120    return 0;
 121 }
 122
 123
 124
 125 Int_t GoodPileupVertices(const Char_t *dir) {
 126    if (gAlice) {
 127        delete AliRunLoader::Instance();
 128        delete gAlice;//if everything was OK here it is already NULL
 129        gAlice = 0x0;
 130    }
 131
 132    Char_t fname[100];
 133    sprintf(fname,"%s/galice.root",dir);
 134
 135    AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON");
 136    if (!rl) {
 137       ::Error("GoodPileupVertices","Can't start session !");
 138       return 1;
 139    }
 140
 141    rl->LoadgAlice();
 142    rl->LoadHeader();
 143    rl->LoadKinematics();
 144
 145
 146    Int_t nev=rl->GetNumberOfEvents();
 147    ::Info("GoodPileupVertices","Number of events : %d\n",nev);
 148
 149    sprintf(fname,"%s/GoodPileupVertices.root",dir);
 150    TFile *refFile=TFile::Open(fname,"recreate");
 151    TClonesArray dummy("AliESDVertex",100), *refs=&dummy;
 152    TTree refTree("refTree","Tree with the reconstructable vertices");
 153    refTree.Branch("Vertices",&refs);
 154
 155    //********  Loop over generated events
 156    for (Int_t e=0; e<nev; e++) {
 157      rl->GetEvent(e);  refFile->cd();
 158
 159      AliHeader *ah=rl->GetHeader();
 160      AliGenCocktailEventHeader *cock=
 161             (AliGenCocktailEventHeader*)ah->GenEventHeader();
 162      TList *headers=cock->GetHeaders();
 163      const Int_t nvtx=headers->GetEntries();
 164      cout<<"Event "<<e<<" Number of vertices: "<<nvtx<<endl;
 165
 166      Int_t nv=0;
 167      for (Int_t v=0; v<nvtx; v++) {
 168          AliGenEventHeader *h=(AliGenEventHeader *)headers->At(v);
 169          TArrayF vtx(3); h->PrimaryVertex(vtx);
 170          //if (...) continue; // Check if this vertex is reconstructable
 171          AliESDVertex *vertex=new ((*refs)[nv]) AliESDVertex();
 172          vertex->SetXv(vtx[0]);
 173          vertex->SetYv(vtx[1]);
 174          vertex->SetZv(vtx[2]);
 175          nv++;
 176      }
 177      refTree.Fill();
 178      refs->Clear();
 179    } //*** end of the loop over generated events
 180
 181    refTree.Write();
 182    refFile->Close();
 183
 184    delete rl;
 185    return 0;
 186 }
 187
 188