[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");   *
 *                                                                          *
 *           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 "AliStack.h"
  #include "AliHeader.h"
  #include "AliGenCocktailEventHeader.h"
  #include "AliRunLoader.h"
  #include "AliRun.h"
  #include "AliESDEvent.h"
  #include "AliESDtrack.h"
#endif

Int_t GoodPileupVertices(const Char_t *dir=".");

extern AliRun *gAlice;
extern TROOT *gROOT;

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

   // **** Book histogramms   
   Int_t nb=20;
   Float_t min=0, max=30.;
   TH2F *h2spd=(TH2F*)gROOT->FindObject("h2spd");
   if (!h2spd) 
      h2spd=new TH2F("h2spd",";Good vertices;Reconstructed vertices",
      nb,min,max, nb,min,max);
   h2spd->SetLineColor(2);
   TH2F *h2trk=(TH2F*)gROOT->FindObject("h2trk");
   if (!h2trk) 
      h2trk=new TH2F("h2trk",";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","Resolution in Z;#DeltaZ (cm);",nb,min,max);
   hzspd->SetLineColor(2);
   TH1F *hztrk=(TH1F*)gROOT->FindObject("hztrk");
   if (!hztrk) 
     hztrk=new TH1F("hztrk","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","Efficiency;Z (cm);Efficiency",nb,min,max);
   heffspd->SetLineColor(2);
   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",";Z (cm);Efficiency",nb,min,max);
   hefftrk->SetLineColor(4);


   // **** 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 e=0;
   while (esdTree->GetEvent(e)) {
     cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
     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: "<<nfoundSPD<<"  good: "<<ngood<<endl;

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

     for (Int_t g=0; g<ngood; g++) {
         const AliESDVertex *vtxg=(AliESDVertex*)refs->UncheckedAt(g);
         Double_t zg=vtxg->GetZv();
         hgood->Fill(zg);

         const AliESDVertex *vtxf=0;
         Double_t zf=0.;
         Int_t f=0;
         for (; f<nfoundSPD; f++) {
             vtxf=(AliESDVertex*)verticesSPD->UncheckedAt(f);
             if (!vtxf->GetStatus()) continue;
             zf=vtxf->GetZv();
             if (TMath::Abs(zf-zg)>2e-2) continue;
             break;
         }
         if (f>=nfoundSPD) {
	     vtxf=event->GetPrimaryVertexSPD();
             if (!vtxf->GetStatus()) goto trk;
             zf=vtxf->GetZv();
             if (TMath::Abs(zf-zg)>2e-2) goto trk;
	 }
         hfoundspd->Fill(zg);
         hzspd->Fill(zf-zg);

     trk:
         for (f=0; f<nfoundTRK; f++) {
             vtxf=(AliESDVertex*)verticesTRK->UncheckedAt(f);
             if (!vtxf->GetStatus()) continue;
             zf=vtxf->GetZv();
             if (TMath::Abs(zf-zg)>2e-2) continue;
             break;
         }
         if (f>=nfoundTRK) {
	     vtxf=event->GetPrimaryVertexTracks();
             if (!vtxf->GetStatus()) continue;
             zf=vtxf->GetZv();
             if (TMath::Abs(zf-zg)>2e-2) continue;
	 }
	 hfoundtrk->Fill(zg);
         hztrk->Fill(zf-zg);

     }

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

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

   refFile->Close();
   
   TCanvas *c1=new TCanvas("c1","",0,0,700,1000);
   c1->Divide(1,3);
   c1->cd(1);
   gPad->SetGridx(); gPad->SetGridy();
   h2spd->Draw("box");
   h2trk->Draw("boxsame");
   TLine *l=new TLine(0,0,30,30);
   l->Draw("same");

   c1->cd(2);
   gPad->SetGridx(); gPad->SetGridy();
   heffspd->Divide(hfoundspd,hgood,1,1,"b");
   heffspd->SetMinimum(0.); heffspd->SetMaximum(1.);
   heffspd->Draw("hist");
   hefftrk->Divide(hfoundtrk,hgood,1,1,"b");
   hefftrk->Draw("histsame");

   c1->cd(3);
   gPad->SetGridx(); gPad->SetGridy();
   hzspd->Draw();
   hztrk->Draw("same");

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

   return 0;
}


Int_t GoodPileupVertices(const Char_t *dir) {
  Int_t FindContributors(Float_t tz, AliStack *stack, Int_t ntrk);
   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();
         Int_t ntrk=FindContributors(t,stack,np);
         if (ntrk<3) continue;
         AliESDVertex *vertex=new ((*refs)[nv]) AliESDVertex();
         vertex->SetXv(vtx[0]);
         vertex->SetYv(vtx[1]);
         vertex->SetZv(vtx[2]);
         vertex->SetNContributors(ntrk);
         nv++;
     }
     refTree.Fill();
     refs->Clear();
   } //*** end of the loop over generated events

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

   delete rl;
   return 0;
}

Int_t FindContributors(Float_t tz, AliStack *stack, Int_t np) {
  Int_t ntrk=0;
  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;
      if (part->Pt()<1) continue;
      Float_t dt=0.5*(tz-part->T())/(tz+part->T());
      if (TMath::Abs(dt)>1e-5) continue;
      ntrk++;
  }
  return ntrk;
}
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	* *
	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>
c4b8f157	16	#include <TParticlePDG.h>
394f72ea	17	#include <TCanvas.h>
394f72ea	18	#include <TFile.h>
c4b8f157	19	#include <TLine.h>
394f72ea	20	#include <TROOT.h>
	21
	22	#include "AliStack.h"
	23	#include "AliHeader.h"
	24	#include "AliGenCocktailEventHeader.h"
	25	#include "AliRunLoader.h"
	26	#include "AliRun.h"
	27	#include "AliESDEvent.h"
	28	#include "AliESDtrack.h"
	29	#endif
	30
	31	Int_t GoodPileupVertices(const Char_t *dir=".");
	32
	33	extern AliRun *gAlice;
	34	extern TROOT *gROOT;
	35
	36	Int_t AliITSUComparisonPileup(const Char_t *dir=".") {
	37	::Info("AliITSUComparisonPileup.C","Doing comparison...");
	38
	39	// **** Book histogramms
c4b8f157	40	Int_t nb=20;
	41	Float_t min=0, max=30.;
	42	TH2F h2spd=(TH2F)gROOT->FindObject("h2spd");
	43	if (!h2spd)
	44	h2spd=new TH2F("h2spd",";Good vertices;Reconstructed vertices",
	45	nb,min,max, nb,min,max);
	46	h2spd->SetLineColor(2);
	47	TH2F h2trk=(TH2F)gROOT->FindObject("h2trk");
	48	if (!h2trk)
	49	h2trk=new TH2F("h2trk",";Good vertices;Reconstructed vertices",
	50	nb,min,max, nb,min,max);
	51	h2trk->SetLineColor(4);
	52
	53	nb=100;
	54	min=-0.03; max=0.03;
	55	TH1F hzspd=(TH1F)gROOT->FindObject("hzspd");
	56	if (!hzspd)
	57	hzspd=new TH1F("hzspd","Resolution in Z;#DeltaZ (cm);",nb,min,max);
	58	hzspd->SetLineColor(2);
	59	TH1F hztrk=(TH1F)gROOT->FindObject("hztrk");
	60	if (!hztrk)
	61	hztrk=new TH1F("hztrk","Resolution in Z;#DeltaZ (cm);",nb,min,max);
	62	hztrk->SetLineColor(4);
	63
	64
	65	nb=30;
	66	min=-10.; max=10.;
	67	TH1F hgood=(TH1F)gROOT->FindObject("hgood");
	68	if (!hgood)
	69	hgood=new TH1F("hgood",";Z (cm);",nb,min,max);
	70	TH1F hfoundspd=(TH1F)gROOT->FindObject("hfoundspd");
	71	if (!hfoundspd)
	72	hfoundspd=new TH1F("hfoundspd",";Z (cm);",nb,min,max);
	73	TH1F heffspd=(TH1F)gROOT->FindObject("heffspd");
	74	if (!heffspd)
	75	heffspd=new TH1F("heffspd","Efficiency;Z (cm);Efficiency",nb,min,max);
	76	heffspd->SetLineColor(2);
	77	TH1F hfoundtrk=(TH1F)gROOT->FindObject("hfoundtrk");
	78	if (!hfoundtrk)
	79	hfoundtrk=new TH1F("hfoundtrk",";Z (cm);",nb,min,max);
	80	TH1F hefftrk=(TH1F)gROOT->FindObject("hefftrk");
	81	if (!hefftrk)
	82	hefftrk=new TH1F("hefftrk",";Z (cm);Efficiency",nb,min,max);
	83	hefftrk->SetLineColor(4);
394f72ea	84
	85
	86	// **** Generate a rerefence file with reconstructable vertices
	87	Char_t fname[100];
	88	sprintf(fname,"%s/GoodPileupVertices.root",dir);
	89	TFile *refFile=TFile::Open(fname,"old");
	90	if (!refFile \|\| !refFile->IsOpen()) {
	91	::Info("AliITSUComparisonPileup.C",
	92	"Marking good pileup vertices (will take a while)...");
	93	if (GoodPileupVertices(dir)) {
	94	::Error("AliITSUComparisonPileup.C","Can't generate the reference file !");
	95	return 1;
	96	}
	97	}
	98	refFile=TFile::Open(fname,"old");
	99	if (!refFile \|\| !refFile->IsOpen()) {
	100	::Error("AliITSUComparisonPileup.C","Can't open the reference file !");
	101	return 1;
	102	}
	103	TTree refTree=(TTree)refFile->Get("refTree");
	104	if (!refTree) {
	105	::Error("AliITSUComparisonPileup.C","Can't get the reference tree !");
	106	return 2;
	107	}
	108	TBranch *branch=refTree->GetBranch("Vertices");
	109	if (!branch) {
	110	::Error("AliITSUComparisonPileup.C","Can't get the vertex branch !");
	111	return 3;
	112	}
	113	TClonesArray dummy("AliESDVertex",100), *refs=&dummy;
	114	branch->SetAddress(&refs);
	115
	116
	117	// **** Open the ESD
	118	sprintf(fname,"%s/AliESDs.root",dir);
	119	TFile *ef=TFile::Open(fname);
	120	if ((!ef)\|\|(!ef->IsOpen())) {
	121	::Error("AliITSUComparisonPileup.C","Can't open AliESDs.root !");
	122	return 4;
	123	}
	124	AliESDEvent* event = new AliESDEvent();
	125	TTree* esdTree = (TTree*) ef->Get("esdTree");
	126	if (!esdTree) {
	127	::Error("AliITSUComparisonPileup.C", "no ESD tree found");
	128	return 6;
	129	}
	130	event->ReadFromTree(esdTree);
	131
	132
c4b8f157	133	//***** Loop over reconstructed events *******
394f72ea	134	Int_t e=0;
	135	while (esdTree->GetEvent(e)) {
	136	cout<<endl<<endl<<"******* Processing event number: "<<e<<"*****\n";
c4b8f157	137	TClonesArray *verticesSPD=event->GetPileupVerticesSPD();
	138	Int_t nfoundSPD=verticesSPD->GetEntries();
	139	TClonesArray *verticesTRK=event->GetPileupVerticesTracks();
	140	Int_t nfoundTRK=verticesTRK->GetEntries();
394f72ea	141
	142	if (refTree->GetEvent(e++)==0) {
	143	cerr<<"No reconstructable vertices for this event !\n";
	144	continue;
	145	}
	146	Int_t ngood=refs->GetEntriesFast();
c4b8f157	147	cout<<"Found SPD: "<<nfoundSPD<<" good: "<<ngood<<endl;
	148
	149	h2spd->Fill(ngood,nfoundSPD);
	150	h2trk->Fill(ngood,nfoundTRK);
	151
	152	for (Int_t g=0; g<ngood; g++) {
	153	const AliESDVertex vtxg=(AliESDVertex)refs->UncheckedAt(g);
	154	Double_t zg=vtxg->GetZv();
	155	hgood->Fill(zg);
	156
	157	const AliESDVertex *vtxf=0;
	158	Double_t zf=0.;
	159	Int_t f=0;
	160	for (; f<nfoundSPD; f++) {
	161	vtxf=(AliESDVertex*)verticesSPD->UncheckedAt(f);
	162	if (!vtxf->GetStatus()) continue;
	163	zf=vtxf->GetZv();
	164	if (TMath::Abs(zf-zg)>2e-2) continue;
	165	break;
	166	}
	167	if (f>=nfoundSPD) {
	168	vtxf=event->GetPrimaryVertexSPD();
	169	if (!vtxf->GetStatus()) goto trk;
	170	zf=vtxf->GetZv();
	171	if (TMath::Abs(zf-zg)>2e-2) goto trk;
	172	}
	173	hfoundspd->Fill(zg);
	174	hzspd->Fill(zf-zg);
	175
	176	trk:
	177	for (f=0; f<nfoundTRK; f++) {
	178	vtxf=(AliESDVertex*)verticesTRK->UncheckedAt(f);
	179	if (!vtxf->GetStatus()) continue;
	180	zf=vtxf->GetZv();
	181	if (TMath::Abs(zf-zg)>2e-2) continue;
	182	break;
	183	}
	184	if (f>=nfoundTRK) {
	185	vtxf=event->GetPrimaryVertexTracks();
	186	if (!vtxf->GetStatus()) continue;
	187	zf=vtxf->GetZv();
	188	if (TMath::Abs(zf-zg)>2e-2) continue;
	189	}
	190	hfoundtrk->Fill(zg);
	191	hztrk->Fill(zf-zg);
394f72ea	192
c4b8f157	193	}
394f72ea	194
c4b8f157	195	} //***** End of the loop over reconstructed events
394f72ea	196
	197	delete event;
	198	delete esdTree;
	199	ef->Close();
	200
	201	refFile->Close();
	202
c4b8f157	203	TCanvas *c1=new TCanvas("c1","",0,0,700,1000);
	204	c1->Divide(1,3);
	205	c1->cd(1);
	206	gPad->SetGridx(); gPad->SetGridy();
	207	h2spd->Draw("box");
	208	h2trk->Draw("boxsame");
	209	TLine *l=new TLine(0,0,30,30);
	210	l->Draw("same");
	211
	212	c1->cd(2);
	213	gPad->SetGridx(); gPad->SetGridy();
	214	heffspd->Divide(hfoundspd,hgood,1,1,"b");
	215	heffspd->SetMinimum(0.); heffspd->SetMaximum(1.);
	216	heffspd->Draw("hist");
	217	hefftrk->Divide(hfoundtrk,hgood,1,1,"b");
	218	hefftrk->Draw("histsame");
	219
	220	c1->cd(3);
	221	gPad->SetGridx(); gPad->SetGridy();
	222	hzspd->Draw();
	223	hztrk->Draw("same");
	224
	225	TFile fc("AliITSUComparisonPileup.root","RECREATE");
394f72ea	226	c1->Write();
	227	fc.Close();
	228
	229	return 0;
	230	}
	231
	232
	233
	234	Int_t GoodPileupVertices(const Char_t *dir) {
c4b8f157	235	Int_t FindContributors(Float_t tz, AliStack *stack, Int_t ntrk);
394f72ea	236	if (gAlice) {
	237	delete AliRunLoader::Instance();
	238	delete gAlice;//if everything was OK here it is already NULL
	239	gAlice = 0x0;
	240	}
	241
	242	Char_t fname[100];
	243	sprintf(fname,"%s/galice.root",dir);
	244
	245	AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON");
	246	if (!rl) {
	247	::Error("GoodPileupVertices","Can't start session !");
	248	return 1;
	249	}
	250
	251	rl->LoadgAlice();
	252	rl->LoadHeader();
	253	rl->LoadKinematics();
	254
	255
	256	Int_t nev=rl->GetNumberOfEvents();
	257	::Info("GoodPileupVertices","Number of events : %d\n",nev);
	258
	259	sprintf(fname,"%s/GoodPileupVertices.root",dir);
	260	TFile *refFile=TFile::Open(fname,"recreate");
	261	TClonesArray dummy("AliESDVertex",100), *refs=&dummy;
	262	TTree refTree("refTree","Tree with the reconstructable vertices");
	263	refTree.Branch("Vertices",&refs);
	264
	265	//******** Loop over generated events
	266	for (Int_t e=0; e<nev; e++) {
	267	rl->GetEvent(e); refFile->cd();
c4b8f157	268	AliStack* stack = rl->Stack();
394f72ea	269
	270	AliHeader *ah=rl->GetHeader();
	271	AliGenCocktailEventHeader *cock=
	272	(AliGenCocktailEventHeader*)ah->GenEventHeader();
	273	TList *headers=cock->GetHeaders();
	274	const Int_t nvtx=headers->GetEntries();
c4b8f157	275	const Int_t np=ah->GetNtrack();
	276	cout<<"Event "<<e<<" Number of vertices, particles: "
	277	<<nvtx<<' '<<np<<endl;
394f72ea	278
	279	Int_t nv=0;
	280	for (Int_t v=0; v<nvtx; v++) {
	281	AliGenEventHeader h=(AliGenEventHeader )headers->At(v);
	282	TArrayF vtx(3); h->PrimaryVertex(vtx);
c4b8f157	283	Float_t t=h->InteractionTime();
	284	Int_t ntrk=FindContributors(t,stack,np);
	285	if (ntrk<3) continue;
394f72ea	286	AliESDVertex vertex=new ((refs)[nv]) AliESDVertex();
	287	vertex->SetXv(vtx[0]);
	288	vertex->SetYv(vtx[1]);
	289	vertex->SetZv(vtx[2]);
c4b8f157	290	vertex->SetNContributors(ntrk);
394f72ea	291	nv++;
	292	}
	293	refTree.Fill();
	294	refs->Clear();
	295	} //*** end of the loop over generated events
	296
	297	refTree.Write();
	298	refFile->Close();
	299
	300	delete rl;
	301	return 0;
	302	}
	303
c4b8f157	304	Int_t FindContributors(Float_t tz, AliStack *stack, Int_t np) {
	305	Int_t ntrk=0;
	306	for (Int_t i=0; i<np; i++) {
	307	if (!stack->IsPhysicalPrimary(i)) continue;
	308	TParticle *part=stack->Particle(i);
	309	if (!part) continue;
	310	TParticlePDG *partPDG = part->GetPDG();
	311	if (!partPDG) continue;
	312	if (TMath::Abs(partPDG->Charge())<1e-10) continue;
	313	if (part->Pt()<1) continue;
	314	Float_t dt=0.5*(tz-part->T())/(tz+part->T());
	315	if (TMath::Abs(dt)>1e-5) continue;
	316	ntrk++;
	317	}
	318	return ntrk;
	319	}