[u/mrichter/AliRoot.git] / ITS / UPGRADE / macros / QA / AliITSUComparisonPileup.C

/****************************************************************************
 *  This macro calculates the efficiency of pileup reconstruction.          *
 *  Works only for events generated with the AliGenPileup generator.        *
 *                                                                          *
 *  Before running, load the ITSU libraries:                                *
 *  gSystem->Load("libITSUpgradeBase");gSystem->Load("libITSUpgradeRec");   *
 *                                                                          *
 *  Definifions:                                                            *
 *  1) Reconstructable track: physical primary, charged, pT > pTmin         * 
 *  2) Reconstructable vertex: has at least nMin reconstructable tracks     *
 *  3) Associated vertex: has at least nAssMin correctly associated tracks  *
 *  4) Good associated vertex: the fraction of correctly associated tracks  *
 *        is at least fracMin                                               *
 *  5) Fake associated vertex: not a good associated vertex                 *
 *  6) Efficiency:  the ratio of 4) over 3)                                 *
 *  7) Fake rate:   the ratio of 5) over 3)                                 *
 *                                                                          *
 *           Origin: I.Belikov, IPHC, Iouri.Belikov@iphc.cnrs.fr            *
 ****************************************************************************/

#if !defined(__CINT__) || defined(__MAKECINT__)
  #include <Riostream.h>
  #include <TMath.h>
  #include <TTree.h>
  #include <TParticle.h>
  #include <TParticlePDG.h>
  #include <TCanvas.h>
  #include <TFile.h>
  #include <TLine.h>
  #include <TROOT.h>
  #include <TStyle.h>
  #include <TLegend.h>

  #include "AliStack.h"
  #include "AliHeader.h"
  #include "AliGenCocktailEventHeader.h"
  #include "AliRunLoader.h"
  #include "AliRun.h"
  #include "AliESDEvent.h"
  #include "AliESDtrack.h"
#endif

//**** Parameters used in the definitions
const Float_t pTmin=0.2;   // Minimal pT for a reconstructable track
const Int_t nMin=3;        // Minimal N of reconstructable tracks per vertex
const Int_t nAssMin=2;     // Minimal number of correctly associated tracks
const Float_t fracMin=0.8; // Minimal fraction of correctly associated tracks

extern AliRun *gAlice;
extern TROOT *gROOT;
extern TStyle *gStyle;

Int_t AliITSUComparisonPileup(const Char_t *dir=".") {
   ::Info("AliITSUComparisonPileup.C","Doing comparison...");
   Int_t GoodPileupVertices(const Char_t *dir=".");

   // **** Book histogramms   
   Int_t nb=35;
   Float_t min=0, max=70.;
   TH2F *h2spd=(TH2F*)gROOT->FindObject("h2spd");
   if (!h2spd) 
     h2spd=new TH2F("h2spd","SPD vertices;Number of good vertices;Number of reconstructed vertices",nb,min,max, nb,min,max);
   h2spd->SetLineColor(2);
   TH2F *h2trk=(TH2F*)gROOT->FindObject("h2trk");
   if (!h2trk) 
     h2trk=new TH2F("h2trk","TRK vertices;Good vertices;Reconstructed vertices",
     nb,min,max, nb,min,max);
   h2trk->SetLineColor(4);

   nb=100;
   min=-0.03; max=0.03;
   TH1F *hzspd=(TH1F*)gROOT->FindObject("hzspd");
   if (!hzspd) 
     hzspd=new TH1F("hzspd","SPD resolution in Z;#DeltaZ (cm);",nb,min,max);
   hzspd->SetLineColor(2);
   TH1F *hztrk=(TH1F*)gROOT->FindObject("hztrk");
   if (!hztrk) 
     hztrk=new TH1F("hztrk","TRK resolution in Z;#DeltaZ (cm);",nb,min,max);
   hztrk->SetLineColor(4);

   
   nb=30;
   min=-10.; max=10.; 
   TH1F *hgood=(TH1F*)gROOT->FindObject("hgood");
   if (!hgood) 
     hgood=new TH1F("hgood",";Z (cm);",nb,min,max);

   TH1F *hfoundspd=(TH1F*)gROOT->FindObject("hfoundspd");
   if (!hfoundspd) 
    hfoundspd=new TH1F("hfoundspd",";Z (cm);",nb,min,max);
   TH1F *heffspd=(TH1F*)gROOT->FindObject("heffspd");
   if (!heffspd) 
      heffspd=new TH1F("heffspd","SPD efficiency + fake rate;Z position of a prim. vertex (cm);Efficiency",nb,min,max);
   heffspd->SetLineColor(2);
   heffspd->Sumw2();

   TH1F *hfoundtrk=(TH1F*)gROOT->FindObject("hfoundtrk");
   if (!hfoundtrk) 
    hfoundtrk=new TH1F("hfoundtrk",";Z (cm);",nb,min,max);
   TH1F *hefftrk=(TH1F*)gROOT->FindObject("hefftrk");
   if (!hefftrk) 
      hefftrk=new TH1F("hefftrk","TRK efficiency;Z (cm);Efficiency",nb,min,max);
   hefftrk->SetLineColor(4);
   hefftrk->Sumw2();

   TH1F *hfaketrk=(TH1F*)gROOT->FindObject("hfaketrk");
   if (!hfaketrk) 
    hfaketrk=new TH1F("hfaketrk",";Z (cm);",nb,min,max);
   TH1F *heffaketrk=(TH1F*)gROOT->FindObject("heffaketrk");
   if (!heffaketrk) 
      heffaketrk=new TH1F("heffaketrk","TRK fake rate;Z (cm);Fake rate",nb,min,max);
   heffaketrk->SetLineColor(4);
   heffaketrk->SetFillColor(590);
   heffaketrk->Sumw2();


   nb=51;
   min=-0.5; max=50.5;
   TH1F *hngood=(TH1F*)gROOT->FindObject("hngood");
   if (!hngood) 
      hngood=new TH1F("hngood",";Z (cm);",nb,min,max);
   TH1F *hnfoundtrk=(TH1F*)gROOT->FindObject("hnfoundtrk");
   if (!hnfoundtrk) 
      hnfoundtrk=new TH1F("hnfoundtrk",";Z (cm);",nb,min,max);
   TH1F *hnefftrk=(TH1F*)gROOT->FindObject("hnefftrk");
   if (!hnefftrk) 
      hnefftrk=new TH1F("hnefftrk","TRK efficiency;Number of tracks;Efficiency",nb,min,max);
   hnefftrk->SetLineColor(4);
   hnefftrk->Sumw2();

   TH1F *hnfaketrk=(TH1F*)gROOT->FindObject("hnfaketrk");
   if (!hnfaketrk) 
      hnfaketrk=new TH1F("hnfaketrk",";Z (cm);",nb,min,max);
   TH1F *hneffaketrk=(TH1F*)gROOT->FindObject("hneffaketrk");
   if (!hneffaketrk) 
      hneffaketrk=new TH1F("hneffaketrk","TRK fake rate;Number of tracks;Efficiency",nb,min,max);
   hneffaketrk->SetLineColor(4);
   hneffaketrk->SetFillColor(590);
   hneffaketrk->Sumw2();


   // **** Generate a rerefence file with reconstructable vertices
   Char_t fname[100];
   sprintf(fname,"%s/GoodPileupVertices.root",dir);
   TFile *refFile=TFile::Open(fname,"old");
   if (!refFile || !refFile->IsOpen()) {
      ::Info("AliITSUComparisonPileup.C",
      "Marking good pileup vertices (will take a while)...");
      if (GoodPileupVertices(dir)) {
     ::Error("AliITSUComparisonPileup.C","Can't generate the reference file !");
         return 1;
      }
   }
   refFile=TFile::Open(fname,"old");
   if (!refFile || !refFile->IsOpen()) {
     ::Error("AliITSUComparisonPileup.C","Can't open the reference file !");
     return 1;
   }   
   TTree *refTree=(TTree*)refFile->Get("refTree");
   if (!refTree) {
     ::Error("AliITSUComparisonPileup.C","Can't get the reference tree !");
     return 2;
   }
   TBranch *branch=refTree->GetBranch("Vertices");
   if (!branch) {
     ::Error("AliITSUComparisonPileup.C","Can't get the vertex branch !");
     return 3;
   }
   TClonesArray dummy("AliESDVertex",100), *refs=&dummy;
   branch->SetAddress(&refs);    


   // **** Open the ESD 
   sprintf(fname,"%s/AliESDs.root",dir);
   TFile *ef=TFile::Open(fname);
   if ((!ef)||(!ef->IsOpen())) {
      ::Error("AliITSUComparisonPileup.C","Can't open AliESDs.root !");
      return 4;
   }
   AliESDEvent* event = new AliESDEvent();
   TTree* esdTree = (TTree*) ef->Get("esdTree");
   if (!esdTree) {
      ::Error("AliITSUComparisonPileup.C", "no ESD tree found");
      return 6;
   }
   event->ReadFromTree(esdTree);


   //******* Loop over reconstructed events *********
   Int_t ntrk=0, ntrkcor=0, ntrkwro=0;
   Int_t e=0;
   while (esdTree->GetEvent(e)) {
     cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
     Int_t nn=event->GetNumberOfTracks();
     ntrk += nn;

     TClonesArray *verticesSPD=event->GetPileupVerticesSPD();
     Int_t nfoundSPD=verticesSPD->GetEntries(); 
     TClonesArray *verticesTRK=event->GetPileupVerticesTracks();
     Int_t nfoundTRK=verticesTRK->GetEntries(); 

     if (refTree->GetEvent(e++)==0) {
        cerr<<"No reconstructable vertices for this event !\n";
        continue;
     }
     Int_t ngood=refs->GetEntriesFast(); 
     cout<<"Found SPD vertices: "<<nfoundSPD<<
           "  Reconstructable vertices: "<<ngood<<endl;

     h2spd->Fill(ngood,nfoundSPD);
     h2trk->Fill(ngood,nfoundTRK);

     Int_t ncor=0, nwro=0;
     for (Int_t g=0; g<ngood; g++) {
         Int_t Associate(const AliESDVertex *g, const AliESDVertex *f, 
            const AliESDEvent *esd); 
         const AliESDVertex *vtxg=(AliESDVertex*)refs->UncheckedAt(g);
         Double_t zg=vtxg->GetZ();
         Double_t ng=vtxg->GetNIndices();
         hgood->Fill(zg);
         hngood->Fill(ng);

         AliESDVertex *vtxf=0;
         Double_t zf=0.;
         Int_t f=0;
         for (; f<nfoundSPD; f++) {
             vtxf=(AliESDVertex*)verticesSPD->UncheckedAt(f);
             if (!vtxf->GetStatus()) continue;
             if (Associate(vtxg,vtxf,event)==0) continue; 
             break;
         }
         if (f>=nfoundSPD) {
	     vtxf=(AliESDVertex *)event->GetPrimaryVertexSPD();
             if (!vtxf->GetStatus()) goto trk;
             if (Associate(vtxg,vtxf,event)==0) goto trk; 
	 }

         zf=vtxf->GetZ();
         hfoundspd->Fill(zg);
         hzspd->Fill(zf-zg);

     trk:
         Int_t n=0;
         for (f=0; f<nfoundTRK; f++) {
             vtxf=(AliESDVertex*)verticesTRK->UncheckedAt(f);
             if (!vtxf->GetStatus()) continue;
             n=Associate(vtxg,vtxf,event);
	     if (n < nAssMin) continue;
             break;
         }
         if (f>=nfoundTRK) {
	     vtxf=(AliESDVertex*)event->GetPrimaryVertexTracks();
             if (!vtxf->GetStatus()) continue;
             n=Associate(vtxg,vtxf,event);
             if (n < nAssMin) continue;
	 }

         ncor+=n;
         nwro+=(vtxf->GetNIndices()-n); 
         zf=vtxf->GetZ();

         if (Float_t(n)/vtxf->GetNIndices() > fracMin) {
	    hfoundtrk->Fill(zg);
	    hnfoundtrk->Fill(ng);
	 } else {
	    hfaketrk->Fill(zg);
	    hnfaketrk->Fill(ng);
	 }
         hztrk->Fill(zf-zg);

         vtxf->SetNContributors(0); // Mark this vertex as already associated

     }
     // Increase the counter of tracks (not)associated with verices
     ntrkcor += ncor;
     ntrkwro += nwro;

   } //***** End of the loop over reconstructed events

   delete event;
   delete esdTree;
   ef->Close();

   refFile->Close();
   
   cout<<"\nTotal number of found tracks: "<<ntrk<<endl;
   cout<<"Number of tracks correctly associated with vertices: "<<ntrkcor<<endl;
   cout<<"Number of tracks wrongly associated with vertices: "<<ntrkwro<<endl;
   if (ntrk != 0) {
     cout<<"Correctly associated/Found:\t"<<Float_t(ntrkcor)/ntrk<<endl;
     cout<<"Wrongly associated/Found:\t"<<Float_t(ntrkwro)/ntrk<<endl;
     cout<<"Not associated/Found:\t\t"<<1.-Float_t(ntrkwro+ntrkcor)/ntrk<<endl;
   }

   gStyle->SetOptStat(0);
   gStyle->SetOptTitle(0);
   TCanvas *c1=new TCanvas("c1","",0,0,700,1000);
   c1->Divide(1,2);
   c1->cd(1);
   gPad->SetGridx(); gPad->SetGridy();
   h2spd->Draw("box");
   h2trk->Draw("boxsame");
   gPad->BuildLegend(0.13,0.65,0.46,0.86)->SetFillColor(0);
   TLine *l=new TLine(0,0,70,70);
   l->Draw("same");

   c1->cd(2);
   gPad->SetGridx(); gPad->SetGridy();
   hztrk->Draw();
   hzspd->Draw("same");
   gPad->BuildLegend(0.13,0.65,0.46,0.86)->SetFillColor(0);


   TCanvas *c2=new TCanvas("c2","",0,0,700,1000);
   c2->Divide(1,2);
   c2->cd(1);
   gPad->SetGridx(); gPad->SetGridy();
   heffspd->Divide(hfoundspd,hgood,1,1,"b");
   heffspd->SetMinimum(0.); heffspd->SetMaximum(1.2);
   heffspd->Draw("ehist");
   hefftrk->Divide(hfoundtrk,hgood,1,1,"b");
   hefftrk->Draw("ehistsame");
   heffaketrk->Divide(hfaketrk,hgood,1,1,"b");
   heffaketrk->Draw("ehistsame");
   gPad->BuildLegend(0.13,0.65,0.46,0.86)->SetFillColor(0);

   c2->cd(2);
   hnefftrk->Divide(hnfoundtrk,hngood,1,1,"b");
   hnefftrk->SetMinimum(0.); hnefftrk->SetMaximum(1.2);
   hneffaketrk->Divide(hnfaketrk,hngood,1,1,"b");
   gPad->SetGridx(); gPad->SetGridy();
   hnefftrk->Draw("ehist");
   hneffaketrk->Draw("ehistsame");
   gPad->BuildLegend(0.13,0.65,0.46,0.86)->SetFillColor(0);


   TFile fc("AliITSUComparisonPileup.root","RECREATE");
   c1->Write();
   c2->Write();
   fc.Close();

   return 0;
}

Int_t 
Associate(const AliESDVertex *g,const AliESDVertex *f,const AliESDEvent *esd) { 
   UShort_t *idxg=g->GetIndices(); Int_t ng=g->GetNIndices(); 
   UShort_t *idxf=f->GetIndices(); Int_t nf=f->GetNIndices();

   if (nf==0) { 
   // SPD vertex
       Double_t zg=g->GetZ();
       Double_t zf=f->GetZ();
       if (TMath::Abs(zf-zg)>2e-2) return 0;
       return 1;
   }
   // TRK vertex
   Int_t nass=0;
   for (Int_t i=0; i<ng; i++) {
       UShort_t labg=idxg[i];
       for (Int_t j=0; j<nf; j++) {
           const AliESDtrack *t=esd->GetTrack(idxf[j]);
           UShort_t labf=TMath::Abs(t->GetLabel());
           if (labg != labf) continue;
           nass++;
           break; 
       }
   } 

   return nass;
}

Int_t GoodPileupVertices(const Char_t *dir) {
  Bool_t FindContributors(Float_t tz, AliStack *stack, UShort_t *idx, Int_t &n);
   if (gAlice) { 
       delete AliRunLoader::Instance();
       delete gAlice;//if everything was OK here it is already NULL
       gAlice = 0x0;
   }

   Char_t fname[100];
   sprintf(fname,"%s/galice.root",dir);

   AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON");
   if (!rl) {
      ::Error("GoodPileupVertices","Can't start session !");
      return 1;
   }

   rl->LoadgAlice();
   rl->LoadHeader();
   rl->LoadKinematics();


   Int_t nev=rl->GetNumberOfEvents();
   ::Info("GoodPileupVertices","Number of events : %d\n",nev);  

   sprintf(fname,"%s/GoodPileupVertices.root",dir);
   TFile *refFile=TFile::Open(fname,"recreate");
   TClonesArray dummy("AliESDVertex",100), *refs=&dummy;
   TTree refTree("refTree","Tree with the reconstructable vertices");
   refTree.Branch("Vertices",&refs);

   //********  Loop over generated events 
   for (Int_t e=0; e<nev; e++) {
     rl->GetEvent(e);  refFile->cd();
     AliStack* stack = rl->Stack();

     AliHeader *ah=rl->GetHeader();
     AliGenCocktailEventHeader *cock=
            (AliGenCocktailEventHeader*)ah->GenEventHeader();
     TList *headers=cock->GetHeaders();
     const Int_t nvtx=headers->GetEntries();
     const Int_t np=ah->GetNtrack();
     cout<<"Event "<<e<<" Number of vertices, particles: "
	 <<nvtx<<' '<<np<<endl;

     Int_t nv=0;
     for (Int_t v=0; v<nvtx; v++) {
         AliGenEventHeader *h=(AliGenEventHeader *)headers->At(v);
         TArrayF vtx(3); h->PrimaryVertex(vtx);
         Float_t t=h->InteractionTime();
         UShort_t *idx=new UShort_t[np];
         Int_t ntrk=0;
         if (!FindContributors(t,stack,idx,ntrk)) {delete[] idx; continue;}
         AliESDVertex *vertex=new ((*refs)[nv]) AliESDVertex();
         vertex->SetXv(vtx[0]);
         vertex->SetYv(vtx[1]);
         vertex->SetZv(vtx[2]);
         vertex->SetNContributors(ntrk);
         vertex->SetIndices(ntrk,idx);
         delete[] idx;
         nv++;
     }
     refTree.Fill();
     refs->Clear();
   } //*** end of the loop over generated events

   refTree.Write();
   refFile->Close();

   delete rl;
   return 0;
}

Bool_t FindContributors(Float_t tz, AliStack *stack, UShort_t *idx, Int_t &n) {
  Int_t ntrk=0;
  Int_t np=stack->GetNtrack();
  for (Int_t i=0; i<np; i++) {
      if (!stack->IsPhysicalPrimary(i)) continue;
      TParticle *part=stack->Particle(i);
      if (!part) continue;
      TParticlePDG *partPDG = part->GetPDG();
      if (!partPDG) continue;
      if (TMath::Abs(partPDG->Charge())<1e-10) continue;
      Float_t dt=0.5*(tz-part->T())/(tz+part->T());
      if (TMath::Abs(dt)>1e-5) continue;
      idx[n++]=i;
      if (part->Pt() > pTmin) ntrk++;
  }
  return (ntrk<nMin) ? kFALSE : kTRUE;
}
Commit	Line	Data
394f72ea	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	* *
f00f355d	8	* Definifions: *
	9	* 1) Reconstructable track: physical primary, charged, pT > pTmin *
	10	* 2) Reconstructable vertex: has at least nMin reconstructable tracks *
	11	* 3) Associated vertex: has at least nAssMin correctly associated tracks *
	12	* 4) Good associated vertex: the fraction of correctly associated tracks *
	13	* is at least fracMin *
	14	* 5) Fake associated vertex: not a good associated vertex *
	15	* 6) Efficiency: the ratio of 4) over 3) *
	16	* 7) Fake rate: the ratio of 5) over 3) *
	17	* *
394f72ea	18	* Origin: I.Belikov, IPHC, Iouri.Belikov@iphc.cnrs.fr *
	19	****************************************************************************/
	20
	21	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	22	#include <Riostream.h>
	23	#include <TMath.h>
	24	#include <TTree.h>
	25	#include <TParticle.h>
c4b8f157	26	#include <TParticlePDG.h>
394f72ea	27	#include <TCanvas.h>
394f72ea	28	#include <TFile.h>
c4b8f157	29	#include <TLine.h>
394f72ea	30	#include <TROOT.h>
f00f355d	31	#include <TStyle.h>
f00f355d	32	#include <TLegend.h>
394f72ea	33
	34	#include "AliStack.h"
	35	#include "AliHeader.h"
	36	#include "AliGenCocktailEventHeader.h"
	37	#include "AliRunLoader.h"
	38	#include "AliRun.h"
	39	#include "AliESDEvent.h"
	40	#include "AliESDtrack.h"
	41	#endif
	42
f00f355d	43	//**** Parameters used in the definitions
	44	const Float_t pTmin=0.2; // Minimal pT for a reconstructable track
	45	const Int_t nMin=3; // Minimal N of reconstructable tracks per vertex
	46	const Int_t nAssMin=2; // Minimal number of correctly associated tracks
	47	const Float_t fracMin=0.8; // Minimal fraction of correctly associated tracks
394f72ea	48
	49	extern AliRun *gAlice;
	50	extern TROOT *gROOT;
f00f355d	51	extern TStyle *gStyle;
394f72ea	52
	53	Int_t AliITSUComparisonPileup(const Char_t *dir=".") {
	54	::Info("AliITSUComparisonPileup.C","Doing comparison...");
f00f355d	55	Int_t GoodPileupVertices(const Char_t *dir=".");
394f72ea	56
394f72ea	57	// **** Book histogramms
e4607b46	58	Int_t nb=35;
e4607b46	59	Float_t min=0, max=70.;
c4b8f157	60	TH2F h2spd=(TH2F)gROOT->FindObject("h2spd");
c4b8f157	61	if (!h2spd)
f00f355d	62	h2spd=new TH2F("h2spd","SPD vertices;Number of good vertices;Number of reconstructed vertices",nb,min,max, nb,min,max);
c4b8f157	63	h2spd->SetLineColor(2);
	64	TH2F h2trk=(TH2F)gROOT->FindObject("h2trk");
	65	if (!h2trk)
f00f355d	66	h2trk=new TH2F("h2trk","TRK vertices;Good vertices;Reconstructed vertices",
f00f355d	67	nb,min,max, nb,min,max);
c4b8f157	68	h2trk->SetLineColor(4);
	69
	70	nb=100;
	71	min=-0.03; max=0.03;
	72	TH1F hzspd=(TH1F)gROOT->FindObject("hzspd");
	73	if (!hzspd)
f00f355d	74	hzspd=new TH1F("hzspd","SPD resolution in Z;#DeltaZ (cm);",nb,min,max);
c4b8f157	75	hzspd->SetLineColor(2);
	76	TH1F hztrk=(TH1F)gROOT->FindObject("hztrk");
	77	if (!hztrk)
f00f355d	78	hztrk=new TH1F("hztrk","TRK resolution in Z;#DeltaZ (cm);",nb,min,max);
c4b8f157	79	hztrk->SetLineColor(4);
	80
	81
	82	nb=30;
	83	min=-10.; max=10.;
	84	TH1F hgood=(TH1F)gROOT->FindObject("hgood");
	85	if (!hgood)
	86	hgood=new TH1F("hgood",";Z (cm);",nb,min,max);
f00f355d	87
c4b8f157	88	TH1F hfoundspd=(TH1F)gROOT->FindObject("hfoundspd");
	89	if (!hfoundspd)
	90	hfoundspd=new TH1F("hfoundspd",";Z (cm);",nb,min,max);
	91	TH1F heffspd=(TH1F)gROOT->FindObject("heffspd");
	92	if (!heffspd)
f00f355d	93	heffspd=new TH1F("heffspd","SPD efficiency + fake rate;Z position of a prim. vertex (cm);Efficiency",nb,min,max);
c4b8f157	94	heffspd->SetLineColor(2);
35b6fd78	95	heffspd->Sumw2();
f00f355d	96
c4b8f157	97	TH1F hfoundtrk=(TH1F)gROOT->FindObject("hfoundtrk");
	98	if (!hfoundtrk)
	99	hfoundtrk=new TH1F("hfoundtrk",";Z (cm);",nb,min,max);
	100	TH1F hefftrk=(TH1F)gROOT->FindObject("hefftrk");
	101	if (!hefftrk)
f00f355d	102	hefftrk=new TH1F("hefftrk","TRK efficiency;Z (cm);Efficiency",nb,min,max);
c4b8f157	103	hefftrk->SetLineColor(4);
35b6fd78	104	hefftrk->Sumw2();
394f72ea	105
f00f355d	106	TH1F hfaketrk=(TH1F)gROOT->FindObject("hfaketrk");
	107	if (!hfaketrk)
	108	hfaketrk=new TH1F("hfaketrk",";Z (cm);",nb,min,max);
	109	TH1F heffaketrk=(TH1F)gROOT->FindObject("heffaketrk");
	110	if (!heffaketrk)
	111	heffaketrk=new TH1F("heffaketrk","TRK fake rate;Z (cm);Fake rate",nb,min,max);
	112	heffaketrk->SetLineColor(4);
	113	heffaketrk->SetFillColor(590);
35b6fd78	114	heffaketrk->Sumw2();
f00f355d	115
f00f355d	116
394f72ea	117
40473af2	118	nb=51;
	119	min=-0.5; max=50.5;
	120	TH1F hngood=(TH1F)gROOT->FindObject("hngood");
	121	if (!hngood)
	122	hngood=new TH1F("hngood",";Z (cm);",nb,min,max);
	123	TH1F hnfoundtrk=(TH1F)gROOT->FindObject("hnfoundtrk");
	124	if (!hnfoundtrk)
	125	hnfoundtrk=new TH1F("hnfoundtrk",";Z (cm);",nb,min,max);
	126	TH1F hnefftrk=(TH1F)gROOT->FindObject("hnefftrk");
	127	if (!hnefftrk)
	128	hnefftrk=new TH1F("hnefftrk","TRK efficiency;Number of tracks;Efficiency",nb,min,max);
	129	hnefftrk->SetLineColor(4);
35b6fd78	130	hnefftrk->Sumw2();
40473af2	131
	132	TH1F hnfaketrk=(TH1F)gROOT->FindObject("hnfaketrk");
	133	if (!hnfaketrk)
	134	hnfaketrk=new TH1F("hnfaketrk",";Z (cm);",nb,min,max);
	135	TH1F hneffaketrk=(TH1F)gROOT->FindObject("hneffaketrk");
	136	if (!hneffaketrk)
	137	hneffaketrk=new TH1F("hneffaketrk","TRK fake rate;Number of tracks;Efficiency",nb,min,max);
	138	hneffaketrk->SetLineColor(4);
	139	hneffaketrk->SetFillColor(590);
35b6fd78	140	hneffaketrk->Sumw2();
40473af2	141
40473af2	142
394f72ea	143	// **** Generate a rerefence file with reconstructable vertices
	144	Char_t fname[100];
	145	sprintf(fname,"%s/GoodPileupVertices.root",dir);
	146	TFile *refFile=TFile::Open(fname,"old");
	147	if (!refFile \|\| !refFile->IsOpen()) {
	148	::Info("AliITSUComparisonPileup.C",
	149	"Marking good pileup vertices (will take a while)...");
	150	if (GoodPileupVertices(dir)) {
	151	::Error("AliITSUComparisonPileup.C","Can't generate the reference file !");
	152	return 1;
	153	}
	154	}
	155	refFile=TFile::Open(fname,"old");
	156	if (!refFile \|\| !refFile->IsOpen()) {
	157	::Error("AliITSUComparisonPileup.C","Can't open the reference file !");
	158	return 1;
	159	}
	160	TTree refTree=(TTree)refFile->Get("refTree");
	161	if (!refTree) {
	162	::Error("AliITSUComparisonPileup.C","Can't get the reference tree !");
	163	return 2;
	164	}
	165	TBranch *branch=refTree->GetBranch("Vertices");
	166	if (!branch) {
	167	::Error("AliITSUComparisonPileup.C","Can't get the vertex branch !");
	168	return 3;
	169	}
	170	TClonesArray dummy("AliESDVertex",100), *refs=&dummy;
	171	branch->SetAddress(&refs);
	172
	173
	174	// **** Open the ESD
	175	sprintf(fname,"%s/AliESDs.root",dir);
	176	TFile *ef=TFile::Open(fname);
	177	if ((!ef)\|\|(!ef->IsOpen())) {
	178	::Error("AliITSUComparisonPileup.C","Can't open AliESDs.root !");
	179	return 4;
	180	}
	181	AliESDEvent* event = new AliESDEvent();
	182	TTree* esdTree = (TTree*) ef->Get("esdTree");
	183	if (!esdTree) {
	184	::Error("AliITSUComparisonPileup.C", "no ESD tree found");
	185	return 6;
	186	}
	187	event->ReadFromTree(esdTree);
	188
	189
c4b8f157	190	//***** Loop over reconstructed events *******
f00f355d	191	Int_t ntrk=0, ntrkcor=0, ntrkwro=0;
394f72ea	192	Int_t e=0;
	193	while (esdTree->GetEvent(e)) {
	194	cout<<endl<<endl<<"******* Processing event number: "<<e<<"*****\n";
f00f355d	195	Int_t nn=event->GetNumberOfTracks();
	196	ntrk += nn;
	197
c4b8f157	198	TClonesArray *verticesSPD=event->GetPileupVerticesSPD();
	199	Int_t nfoundSPD=verticesSPD->GetEntries();
	200	TClonesArray *verticesTRK=event->GetPileupVerticesTracks();
	201	Int_t nfoundTRK=verticesTRK->GetEntries();
394f72ea	202
	203	if (refTree->GetEvent(e++)==0) {
	204	cerr<<"No reconstructable vertices for this event !\n";
	205	continue;
	206	}
	207	Int_t ngood=refs->GetEntriesFast();
f00f355d	208	cout<<"Found SPD vertices: "<<nfoundSPD<<
9d6e1a08	209	" Reconstructable vertices: "<<ngood<<endl;
c4b8f157	210
	211	h2spd->Fill(ngood,nfoundSPD);
	212	h2trk->Fill(ngood,nfoundTRK);
	213
f00f355d	214	Int_t ncor=0, nwro=0;
c4b8f157	215	for (Int_t g=0; g<ngood; g++) {
f00f355d	216	Int_t Associate(const AliESDVertex g, const AliESDVertex f,
e4607b46	217	const AliESDEvent *esd);
c4b8f157	218	const AliESDVertex vtxg=(AliESDVertex)refs->UncheckedAt(g);
15eb37d8	219	Double_t zg=vtxg->GetZ();
40473af2	220	Double_t ng=vtxg->GetNIndices();
c4b8f157	221	hgood->Fill(zg);
40473af2	222	hngood->Fill(ng);
c4b8f157	223
f00f355d	224	AliESDVertex *vtxf=0;
c4b8f157	225	Double_t zf=0.;
	226	Int_t f=0;
	227	for (; f<nfoundSPD; f++) {
	228	vtxf=(AliESDVertex*)verticesSPD->UncheckedAt(f);
	229	if (!vtxf->GetStatus()) continue;
f00f355d	230	if (Associate(vtxg,vtxf,event)==0) continue;
c4b8f157	231	break;
	232	}
	233	if (f>=nfoundSPD) {
f00f355d	234	vtxf=(AliESDVertex *)event->GetPrimaryVertexSPD();
c4b8f157	235	if (!vtxf->GetStatus()) goto trk;
f00f355d	236	if (Associate(vtxg,vtxf,event)==0) goto trk;
c4b8f157	237	}
f00f355d	238
15eb37d8	239	zf=vtxf->GetZ();
c4b8f157	240	hfoundspd->Fill(zg);
	241	hzspd->Fill(zf-zg);
	242
	243	trk:
f00f355d	244	Int_t n=0;
c4b8f157	245	for (f=0; f<nfoundTRK; f++) {
	246	vtxf=(AliESDVertex*)verticesTRK->UncheckedAt(f);
	247	if (!vtxf->GetStatus()) continue;
f00f355d	248	n=Associate(vtxg,vtxf,event);
f00f355d	249	if (n < nAssMin) continue;
c4b8f157	250	break;
	251	}
	252	if (f>=nfoundTRK) {
f00f355d	253	vtxf=(AliESDVertex*)event->GetPrimaryVertexTracks();
c4b8f157	254	if (!vtxf->GetStatus()) continue;
f00f355d	255	n=Associate(vtxg,vtxf,event);
	256	if (n < nAssMin) continue;
	257	}
	258
	259	ncor+=n;
	260	nwro+=(vtxf->GetNIndices()-n);
15eb37d8	261	zf=vtxf->GetZ();
f00f355d	262
	263	if (Float_t(n)/vtxf->GetNIndices() > fracMin) {
	264	hfoundtrk->Fill(zg);
40473af2	265	hnfoundtrk->Fill(ng);
f00f355d	266	} else {
f00f355d	267	hfaketrk->Fill(zg);
40473af2	268	hnfaketrk->Fill(ng);
c4b8f157	269	}
c4b8f157	270	hztrk->Fill(zf-zg);
394f72ea	271
f00f355d	272	vtxf->SetNContributors(0); // Mark this vertex as already associated
f00f355d	273
c4b8f157	274	}
f00f355d	275	// Increase the counter of tracks (not)associated with verices
	276	ntrkcor += ncor;
	277	ntrkwro += nwro;
394f72ea	278
c4b8f157	279	} //***** End of the loop over reconstructed events
394f72ea	280
	281	delete event;
	282	delete esdTree;
	283	ef->Close();
	284
	285	refFile->Close();
	286
f00f355d	287	cout<<"\nTotal number of found tracks: "<<ntrk<<endl;
	288	cout<<"Number of tracks correctly associated with vertices: "<<ntrkcor<<endl;
	289	cout<<"Number of tracks wrongly associated with vertices: "<<ntrkwro<<endl;
	290	if (ntrk != 0) {
	291	cout<<"Correctly associated/Found:\t"<<Float_t(ntrkcor)/ntrk<<endl;
	292	cout<<"Wrongly associated/Found:\t"<<Float_t(ntrkwro)/ntrk<<endl;
	293	cout<<"Not associated/Found:\t\t"<<1.-Float_t(ntrkwro+ntrkcor)/ntrk<<endl;
	294	}
	295
	296	gStyle->SetOptStat(0);
	297	gStyle->SetOptTitle(0);
c4b8f157	298	TCanvas *c1=new TCanvas("c1","",0,0,700,1000);
40473af2	299	c1->Divide(1,2);
c4b8f157	300	c1->cd(1);
	301	gPad->SetGridx(); gPad->SetGridy();
	302	h2spd->Draw("box");
	303	h2trk->Draw("boxsame");
f00f355d	304	gPad->BuildLegend(0.13,0.65,0.46,0.86)->SetFillColor(0);
f00f355d	305	TLine *l=new TLine(0,0,70,70);
c4b8f157	306	l->Draw("same");
	307
	308	c1->cd(2);
	309	gPad->SetGridx(); gPad->SetGridy();
40473af2	310	hztrk->Draw();
	311	hzspd->Draw("same");
	312	gPad->BuildLegend(0.13,0.65,0.46,0.86)->SetFillColor(0);
	313
	314
	315	TCanvas *c2=new TCanvas("c2","",0,0,700,1000);
	316	c2->Divide(1,2);
	317	c2->cd(1);
	318	gPad->SetGridx(); gPad->SetGridy();
c4b8f157	319	heffspd->Divide(hfoundspd,hgood,1,1,"b");
e4607b46	320	heffspd->SetMinimum(0.); heffspd->SetMaximum(1.2);
35b6fd78	321	heffspd->Draw("ehist");
c4b8f157	322	hefftrk->Divide(hfoundtrk,hgood,1,1,"b");
35b6fd78	323	hefftrk->Draw("ehistsame");
f00f355d	324	heffaketrk->Divide(hfaketrk,hgood,1,1,"b");
35b6fd78	325	heffaketrk->Draw("ehistsame");
f00f355d	326	gPad->BuildLegend(0.13,0.65,0.46,0.86)->SetFillColor(0);
c4b8f157	327
40473af2	328	c2->cd(2);
	329	hnefftrk->Divide(hnfoundtrk,hngood,1,1,"b");
	330	hnefftrk->SetMinimum(0.); hnefftrk->SetMaximum(1.2);
	331	hneffaketrk->Divide(hnfaketrk,hngood,1,1,"b");
c4b8f157	332	gPad->SetGridx(); gPad->SetGridy();
35b6fd78	333	hnefftrk->Draw("ehist");
35b6fd78	334	hneffaketrk->Draw("ehistsame");
f00f355d	335	gPad->BuildLegend(0.13,0.65,0.46,0.86)->SetFillColor(0);
c4b8f157	336
40473af2	337
c4b8f157	338	TFile fc("AliITSUComparisonPileup.root","RECREATE");
394f72ea	339	c1->Write();
40473af2	340	c2->Write();
394f72ea	341	fc.Close();
	342
	343	return 0;
	344	}
	345
f00f355d	346	Int_t
e4607b46	347	Associate(const AliESDVertex g,const AliESDVertex f,const AliESDEvent *esd) {
	348	UShort_t *idxg=g->GetIndices(); Int_t ng=g->GetNIndices();
	349	UShort_t *idxf=f->GetIndices(); Int_t nf=f->GetNIndices();
	350
	351	if (nf==0) {
	352	// SPD vertex
15eb37d8	353	Double_t zg=g->GetZ();
15eb37d8	354	Double_t zf=f->GetZ();
f00f355d	355	if (TMath::Abs(zf-zg)>2e-2) return 0;
f00f355d	356	return 1;
e4607b46	357	}
	358	// TRK vertex
	359	Int_t nass=0;
	360	for (Int_t i=0; i<ng; i++) {
	361	UShort_t labg=idxg[i];
	362	for (Int_t j=0; j<nf; j++) {
	363	const AliESDtrack *t=esd->GetTrack(idxf[j]);
	364	UShort_t labf=TMath::Abs(t->GetLabel());
	365	if (labg != labf) continue;
	366	nass++;
	367	break;
	368	}
	369	}
	370
f00f355d	371	return nass;
e4607b46	372	}
394f72ea	373
394f72ea	374	Int_t GoodPileupVertices(const Char_t *dir) {
475afe2c	375	Bool_t FindContributors(Float_t tz, AliStack stack, UShort_t idx, Int_t &n);
394f72ea	376	if (gAlice) {
	377	delete AliRunLoader::Instance();
	378	delete gAlice;//if everything was OK here it is already NULL
	379	gAlice = 0x0;
	380	}
	381
	382	Char_t fname[100];
	383	sprintf(fname,"%s/galice.root",dir);
	384
	385	AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON");
	386	if (!rl) {
	387	::Error("GoodPileupVertices","Can't start session !");
	388	return 1;
	389	}
	390
	391	rl->LoadgAlice();
	392	rl->LoadHeader();
	393	rl->LoadKinematics();
	394
	395
	396	Int_t nev=rl->GetNumberOfEvents();
	397	::Info("GoodPileupVertices","Number of events : %d\n",nev);
	398
	399	sprintf(fname,"%s/GoodPileupVertices.root",dir);
	400	TFile *refFile=TFile::Open(fname,"recreate");
	401	TClonesArray dummy("AliESDVertex",100), *refs=&dummy;
	402	TTree refTree("refTree","Tree with the reconstructable vertices");
	403	refTree.Branch("Vertices",&refs);
	404
	405	//******** Loop over generated events
	406	for (Int_t e=0; e<nev; e++) {
	407	rl->GetEvent(e); refFile->cd();
c4b8f157	408	AliStack* stack = rl->Stack();
394f72ea	409
	410	AliHeader *ah=rl->GetHeader();
	411	AliGenCocktailEventHeader *cock=
	412	(AliGenCocktailEventHeader*)ah->GenEventHeader();
	413	TList *headers=cock->GetHeaders();
	414	const Int_t nvtx=headers->GetEntries();
c4b8f157	415	const Int_t np=ah->GetNtrack();
	416	cout<<"Event "<<e<<" Number of vertices, particles: "
	417	<<nvtx<<' '<<np<<endl;
394f72ea	418
	419	Int_t nv=0;
	420	for (Int_t v=0; v<nvtx; v++) {
	421	AliGenEventHeader h=(AliGenEventHeader )headers->At(v);
	422	TArrayF vtx(3); h->PrimaryVertex(vtx);
c4b8f157	423	Float_t t=h->InteractionTime();
e4607b46	424	UShort_t *idx=new UShort_t[np];
475afe2c	425	Int_t ntrk=0;
475afe2c	426	if (!FindContributors(t,stack,idx,ntrk)) {delete[] idx; continue;}
394f72ea	427	AliESDVertex vertex=new ((refs)[nv]) AliESDVertex();
	428	vertex->SetXv(vtx[0]);
	429	vertex->SetYv(vtx[1]);
	430	vertex->SetZv(vtx[2]);
c4b8f157	431	vertex->SetNContributors(ntrk);
e4607b46	432	vertex->SetIndices(ntrk,idx);
9d6e1a08	433	delete[] idx;
394f72ea	434	nv++;
	435	}
	436	refTree.Fill();
	437	refs->Clear();
	438	} //*** end of the loop over generated events
	439
	440	refTree.Write();
	441	refFile->Close();
	442
	443	delete rl;
	444	return 0;
	445	}
	446
475afe2c	447	Bool_t FindContributors(Float_t tz, AliStack stack, UShort_t idx, Int_t &n) {
c4b8f157	448	Int_t ntrk=0;
e4607b46	449	Int_t np=stack->GetNtrack();
c4b8f157	450	for (Int_t i=0; i<np; i++) {
	451	if (!stack->IsPhysicalPrimary(i)) continue;
	452	TParticle *part=stack->Particle(i);
	453	if (!part) continue;
	454	TParticlePDG *partPDG = part->GetPDG();
	455	if (!partPDG) continue;
	456	if (TMath::Abs(partPDG->Charge())<1e-10) continue;
c4b8f157	457	Float_t dt=0.5*(tz-part->T())/(tz+part->T());
c4b8f157	458	if (TMath::Abs(dt)>1e-5) continue;
475afe2c	459	idx[n++]=i;
475afe2c	460	if (part->Pt() > pTmin) ntrk++;
c4b8f157	461	}
475afe2c	462	return (ntrk<nMin) ? kFALSE : kTRUE;
c4b8f157	463	}