[u/mrichter/AliRoot.git] / STEER / AliKFParticleTest.C


//---------------------------------------------------------------------------------
// The example of usage of AliKFParticle & AliKFVertex classes for V0 analysis
// .
// @author  S.Gorbunov, I.Kisel
// @version 1.0
// @since   13.05.07
// 
// The AliKFParticleTest macro contains a toy V0 finder for ESD tracks.
// At the first step, event primary vertex is reconstructed. 
// At the second step, ideal PID hypothesis are assigned to all the particles.
// At the third step, V0 candidates are constructed for each pair 
// of positive and negative particles.
// V0 candidate considered as good when it passes Chi^2 cut and 
// it is placed >= 3 Sigma away from the primary vertex.
// Invariant mass distribution for all good V0 candidates is plotted.
//
//  -= Copyright &copy ALICE HLT Group =-
//_________________________________________________________________________________


#if !defined( __CINT__) || defined(__MAKECINT__)

  //gSystem->AddIncludePath("-I\"$ALICE_ROOT/include\"");

  #include <Riostream.h>
  #include <TTree.h>
  #include <TFile.h>
  #include <TCanvas.h>
  #include <TH1D.h>
  #include <TLatex.h>
  #include <TDatabasePDG.h>
  #include <TParticle.h>

  #include "AliRun.h"
  #include "AliStack.h"
  #include "AliESDEvent.h"
  #include "AliTracker.h"
  #include "AliKFParticle.h"
  #include "AliKFVertex.h"

#endif


Bool_t old_file = 1;
TCanvas *canvas=0;
TH1D  *histoMass[3]={0,0,0};
Int_t histoPDG [4]= {22,310,3122,421};
char* histoName[4]= {"#gamma","K^{0}_{s}","#Lambda","D^{0}"};
TLatex *histoText[4]={0,0,0,0};

void DrawV0() 
{
  //* Draw the invariant mass histogram

  if( !canvas ) return;
  histoMass[0]->Draw();
  for( Int_t i=1; i<3; i++ ) histoMass[i]->Draw("same");
  for( Int_t i=0; i<4; i++ ){
    Double_t mass = TDatabasePDG::Instance()->GetParticle(histoPDG[i])->Mass();    
    Int_t bin = histoMass[2]->FindBin(mass) -5;
    if( bin<0 ) bin =0;
    Double_t max = 0;
    for( Int_t j=bin; j<bin+10; j++ ) 
    if( max<histoMass[2]->GetBinContent(j) ) max = histoMass[2]->GetBinContent(j);
    histoText[i]->SetY(max);
    histoText[i]->SetX(mass+.05);
    if(max>0) histoText[i]->Draw();
  }
  if( canvas ) canvas->Update();
}

void StartV0()
{
  //* (Re)create histograms and fill them from temporary file
  
  TDirectory *curr = gDirectory;
  canvas = new TCanvas();
  histoMass[0] = new TH1D("massAll","AliKFParticle test", 500,0,3);
  histoMass[0]->SetXTitle("V^{0} invariant mass [GeV]");
  histoMass[1] = new TH1D("massMulti","V^{+-} contributions", 500,0,3);
  histoMass[2] = new TH1D("massV0","V^{0} signal", 500,0,3);
  histoMass[0]->SetLineColor(8);
  histoMass[0]->SetFillColor(8);
  histoMass[1]->SetLineColor(kMagenta);
  histoMass[1]->SetFillColor(kMagenta);
  histoMass[2]->SetLineColor(kBlue);
  histoMass[2]->SetFillColor(kBlue);
  for( Int_t i=0; i<4; i++ ){
    histoText[i] = new TLatex(0,0,histoName[i]);
    histoText[i]->SetTextColor(kBlue);
  }
  TFile* file = new TFile("AliKFParticleTest.root", "READ" );
  if( !old_file && file ){
    file->cd();
    for( Int_t i=0; i<3; i++ ){
      TObject *old = file->FindObjectAny(histoMass[i]->GetName());
      if( old && histoMass[i] ) histoMass[i]->Add( (TH1D*) old);
    }
    file->Close();
  }
  curr->cd();
  DrawV0();
}

void EndV0()
{
  //* End of gAlice -> store all histos in temporary file, clean memory

  TDirectory *curr = gDirectory;
  TFile* file = new TFile("AliKFParticleTest.root", "RECREATE" );
  old_file = 0;
  file->cd();
  for( Int_t i=0; i<3; i++ ) histoMass[i]->Write();
  file->Close();
  for( Int_t i=0; i<3; i++ ) delete histoMass[i];
  for( Int_t i=0; i<4; i++ ) delete histoText[i];
  delete canvas;
  curr->cd();
}

class TESDTrackInfo
{
public:
  TESDTrackInfo(){}
  AliKFParticle fParticle; //* assigned KFParticle
  Bool_t fPrimUsedFlag;    //* flag shows that the particle was used for primary vertex fit
  Bool_t fOK;              //* is the track good enough
  Int_t mcPDG;             //* Monte Carlo PDG code
  Int_t mcMotherID;        //* Monte Carlo ID of its mother
};

void RunV0(  AliESDEvent *event )
{
  //* V0 finder

  static Int_t iEvent = 0;
  cout<<"Event "<<iEvent++<<endl;

  if( !gAlice ) return;
  AliRunLoader *rl = AliRunLoader::Instance(); 
 
  AliStack *stack = rl->Stack();
  if( !stack ) return;  

  Int_t nESDTracks=event->GetNumberOfTracks();
  if( nESDTracks>1000 ) nESDTracks=1000;


  TESDTrackInfo ESDTrackInfo[1000]; //* parallel to ESD tracks

  //* Fill ESDTrackInfo array

  for (Int_t iTr=0; iTr<nESDTracks; iTr++)
    {   
      TESDTrackInfo &info = ESDTrackInfo[iTr];
      info.fOK = 0;
      info.fPrimUsedFlag = 0;
      info.mcPDG = -1;
      info.mcMotherID = -1;

      //* track quality check

      AliESDtrack *pTrack = event->GetTrack(iTr);    
      if( !pTrack  ) continue;
      if (pTrack->GetKinkIndex(0)>0) continue;
      if ( !( pTrack->GetStatus()&AliESDtrack::kITSrefit ) ) continue;
      Int_t indi[12];
      if( pTrack->GetITSclusters(indi) <5 ) continue;
      Int_t PDG = ( pTrack->GetSigned1Pt() <0 ) ?321 :211;

      //* take MC PDG  
      { 
	Int_t mcID = TMath::Abs(pTrack->GetLabel());
	TParticle * part = stack->Particle(TMath::Abs(mcID));
	if( part ){
	  info.mcPDG = part->GetPdgCode();
	  PDG = info.mcPDG;
	  if( mcID>=0 ) info.mcMotherID = part->GetFirstMother();
	}    
      }

      //* Construct KFParticle for the track

      info.fParticle = AliKFParticle( *pTrack, PDG );
      info.fOK = 1;   
    }

  //* Find event primary vertex
  
  AliKFVertex primVtx;  
  {
    const AliKFParticle * vSelected[1000]; //* Selected particles for vertex fit
    Int_t vIndex[1000];                    //* Indices of selected particles
    Bool_t vFlag[1000];                    //* Flags returned by the vertex finder

    Int_t nSelected = 0;
    for( Int_t i = 0; i<nESDTracks; i++){ 
      if(ESDTrackInfo[i].fOK ){
	vSelected[nSelected] = &(ESDTrackInfo[i].fParticle);
	vIndex[nSelected] = i;
	nSelected++;
      }
    }
    primVtx.ConstructPrimaryVertex( vSelected, nSelected, vFlag, 3. );
    for( Int_t i = 0; i<nSelected; i++){ 
      if( vFlag[i] ) ESDTrackInfo[vIndex[i]].fPrimUsedFlag = 1;
    }
    if( primVtx.GetNDF() <1 ) return; //* Less then two tracks in primary vertex 
  }

  //* V0 finder

  for( Int_t iTr = 0; iTr<nESDTracks; iTr++ ){ //* first daughter

    TESDTrackInfo &info = ESDTrackInfo[iTr];
    if( !info.fOK ) continue;    

    for( Int_t jTr = iTr+1; jTr<nESDTracks; jTr++ ){  //* second daughter
      TESDTrackInfo &jnfo = ESDTrackInfo[jTr];
      if( !jnfo.fOK ) continue;
      
      //* check for different charge

      if( info.fParticle.GetQ() == jnfo.fParticle.GetQ() ) continue;      

      //* construct V0 mother

      AliKFParticle V0( info.fParticle, jnfo.fParticle );     

      //* check V0 Chi^2

      if( V0.GetNDF()<1 ) continue;
      if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) >3. ) continue;

      //* subtruct daughters from primary vertex 

      AliKFVertex primVtxCopy = primVtx;
       
      if( info.fPrimUsedFlag ) primVtxCopy -= info.fParticle;
      if( jnfo.fPrimUsedFlag ) primVtxCopy -= jnfo.fParticle;

      //* Check V0 Chi^2 deviation from primary vertex 

      if( V0.GetDeviationFromVertex( primVtxCopy ) >3. ) continue;

      //* Add V0 to primary vertex to improve the primary vertex resolution

      primVtxCopy += V0;      

      //* Set production vertex for V0

      V0.SetProductionVertex( primVtxCopy );

      //* Check chi^2 for a case

      if( TMath::Sqrt( TMath::Abs(V0.GetChi2()/V0.GetNDF())) >3. ) continue;

      //* Get V0 decay length with estimated error

      Double_t length, sigmaLength;
      if( V0.GetDecayLength( length, sigmaLength ) ) continue;

      //* Reject V0 if it decays too close to the primary vertex

      if( length  <3.*sigmaLength ) continue;

      //* Get V0 invariant mass and plot it

      Double_t mass, sigmaMass;
      if( V0.GetMass( mass, sigmaMass ) ) continue;
      histoMass[0]->Fill(mass);            

      //* Fill signal histograms using MC information

      if( info.mcMotherID==jnfo.mcMotherID && info.mcMotherID>=0 ){
	TParticle *mother = stack->Particle(info.mcMotherID);
	if( mother && TMath::Abs(mother->GetPdgCode())!=21 ){
	  histoMass[1]->Fill(mass);
	  if( mother->GetNDaughters()==2 ){
	    histoMass[2]->Fill(mass);	
	  }
	  cout<<"PDG V0,pi,pj, ndaughters, mc mass, reco mass = "<<mother->GetPdgCode()<<","<<info.mcPDG<<","<<jnfo.mcPDG<<", "
	      << mother->GetNDaughters()<<", "<<mother->GetMass()<<", "<<mass<<endl;
	}
      }
    }
  }
  if( iEvent %100 == 0 || (iEvent<100 && iEvent %10==0)) DrawV0();     
}


Int_t AliKFParticleTest(Int_t n1=0,Int_t n2=1000,char *dire="/d/alice10/sma/my_v4-05-Release/pp/"){
  //* Main macro
    
  //  LOOP  OVER  SERIES  OF  DIRECTORIES
  
  for (Int_t ifi=n1; ifi<=n2; ifi++) {
   
    char nstring[5], filename[100], esdfile[100];
    sprintf(nstring,"%3.3d",ifi);
    sprintf(filename,"%s%s/galice.root",dire,nstring);
    sprintf(esdfile,"%s%s/AliESDs.root",dire,nstring);
    
    cout <<"  Opening "<<filename<<"\nand ESD "<<esdfile<<endl;
    
    TFile* ftemp = TFile::Open(filename);
    Printf("ftemp = %p", ftemp);
    if (!ftemp) continue;

    if (gAlice) {
      delete AliRunLoader::Instance();
      delete gAlice;
      gAlice=0;
    }
    
    AliRunLoader *rl = AliRunLoader::Open(filename);
    if ( !rl ) { 
      ::Error("AliKFParticleTest.C","Can not open session !");
      continue;
    }
    if (rl->LoadgAlice()) {
      ::Error("AliKFParticleTest.C","LoadgAlice returned error !");
      delete rl;
      continue;
    }
    if (rl->LoadHeader()) {
      ::Error("AliKFParticleTest.C","LoadHeader returned error !");
      delete rl;
      continue;
    }
    rl->LoadKinematics();

    //---------------------------------------//
    //                                       //
    //               ESD  file               //
    //                                       //
    //---------------------------------------//
    
    // Open file with the ESD
    TFile *ef=TFile::Open(esdfile);

    //Check if the file could be opened
    if (!ef || !ef->IsOpen()) {cerr<<"Error open AliESDs.root !\n"; continue ;}

    //create event object
    AliESDEvent *event = new AliESDEvent;

    //Set pointer to the esd tree in the file
    TTree* tree = (TTree*) ef->Get("esdTree");
    
    //check if the tree exists
    if (!tree) {cerr<<"no ESD tree found\n"; continue;};
    
    //Set pointer to the esd object in the tree
    event->ReadFromTree(tree);
    
    //Number of events
    Int_t nevents=tree->GetEntriesFast();
    cout << "Number of events: " << nevents << endl;

    StartV0();
    for (Int_t iev=0; iev<nevents; iev++){
      tree->GetEvent(iev);
      rl->GetEvent(iev);
      AliKFParticle::SetField(event->GetMagneticField());
      RunV0(event);
    }
    EndV0();
    delete event;
    ef->Close();
  }
  StartV0();
  return 0;
}
Commit	Line	Data
cc7cb084	1
	2	//---------------------------------------------------------------------------------
	3	// The example of usage of AliKFParticle & AliKFVertex classes for V0 analysis
	4	// .
	5	// @author S.Gorbunov, I.Kisel
	6	// @version 1.0
	7	// @since 13.05.07
	8	//
	9	// The AliKFParticleTest macro contains a toy V0 finder for ESD tracks.
	10	// At the first step, event primary vertex is reconstructed.
	11	// At the second step, ideal PID hypothesis are assigned to all the particles.
	12	// At the third step, V0 candidates are constructed for each pair
	13	// of positive and negative particles.
	14	// V0 candidate considered as good when it passes Chi^2 cut and
	15	// it is placed >= 3 Sigma away from the primary vertex.
	16	// Invariant mass distribution for all good V0 candidates is plotted.
	17	//
	18	// -= Copyright &copy ALICE HLT Group =-
	19	//_________________________________________________________________________________
	20
	21
	22
	23	#if !defined( __CINT__) \|\| defined(__MAKECINT__)
	24
	25	//gSystem->AddIncludePath("-I\"$ALICE_ROOT/include\"");
	26
	27	#include <Riostream.h>
	28	#include <TTree.h>
	29	#include <TFile.h>
	30	#include <TCanvas.h>
	31	#include <TH1D.h>
	32	#include <TLatex.h>
	33	#include <TDatabasePDG.h>
	34	#include <TParticle.h>
	35
	36	#include "AliRun.h"
	37	#include "AliStack.h"
86ad5fcb	38	#include "AliESDEvent.h"
cc7cb084	39	#include "AliTracker.h"
	40	#include "AliKFParticle.h"
	41	#include "AliKFVertex.h"
	42
	43	#endif
	44
	45
	46	Bool_t old_file = 1;
	47	TCanvas *canvas=0;
	48	TH1D *histoMass[3]={0,0,0};
	49	Int_t histoPDG [4]= {22,310,3122,421};
	50	char* histoName[4]= {"#gamma","K^{0}_{s}","#Lambda","D^{0}"};
	51	TLatex *histoText[4]={0,0,0,0};
	52
	53	void DrawV0()
	54	{
	55	//* Draw the invariant mass histogram
	56
	57	if( !canvas ) return;
	58	histoMass[0]->Draw();
	59	for( Int_t i=1; i<3; i++ ) histoMass[i]->Draw("same");
	60	for( Int_t i=0; i<4; i++ ){
	61	Double_t mass = TDatabasePDG::Instance()->GetParticle(histoPDG[i])->Mass();
	62	Int_t bin = histoMass[2]->FindBin(mass) -5;
	63	if( bin<0 ) bin =0;
	64	Double_t max = 0;
	65	for( Int_t j=bin; j<bin+10; j++ )
	66	if( max<histoMass[2]->GetBinContent(j) ) max = histoMass[2]->GetBinContent(j);
	67	histoText[i]->SetY(max);
	68	histoText[i]->SetX(mass+.05);
	69	if(max>0) histoText[i]->Draw();
	70	}
	71	if( canvas ) canvas->Update();
	72	}
	73
	74	void StartV0()
	75	{
	76	//* (Re)create histograms and fill them from temporary file
	77
	78	TDirectory *curr = gDirectory;
	79	canvas = new TCanvas();
	80	histoMass[0] = new TH1D("massAll","AliKFParticle test", 500,0,3);
	81	histoMass[0]->SetXTitle("V^{0} invariant mass [GeV]");
	82	histoMass[1] = new TH1D("massMulti","V^{+-} contributions", 500,0,3);
	83	histoMass[2] = new TH1D("massV0","V^{0} signal", 500,0,3);
	84	histoMass[0]->SetLineColor(8);
	85	histoMass[0]->SetFillColor(8);
	86	histoMass[1]->SetLineColor(kMagenta);
	87	histoMass[1]->SetFillColor(kMagenta);
	88	histoMass[2]->SetLineColor(kBlue);
	89	histoMass[2]->SetFillColor(kBlue);
	90	for( Int_t i=0; i<4; i++ ){
	91	histoText[i] = new TLatex(0,0,histoName[i]);
	92	histoText[i]->SetTextColor(kBlue);
	93	}
	94	TFile* file = new TFile("AliKFParticleTest.root", "READ" );
	95	if( !old_file && file ){
	96	file->cd();
	97	for( Int_t i=0; i<3; i++ ){
	98	TObject *old = file->FindObjectAny(histoMass[i]->GetName());
	99	if( old && histoMass[i] ) histoMass[i]->Add( (TH1D*) old);
	100	}
	101	file->Close();
	102	}
103	curr->cd();
104	DrawV0();
105	}
106
107	void EndV0()
108	{
109	//* End of gAlice -> store all histos in temporary file, clean memory
110
111	TDirectory *curr = gDirectory;
112	TFile* file = new TFile("AliKFParticleTest.root", "RECREATE" );
113	old_file = 0;
114	file->cd();
115	for( Int_t i=0; i<3; i++ ) histoMass[i]->Write();
116	file->Close();
117	for( Int_t i=0; i<3; i++ ) delete histoMass[i];
118	for( Int_t i=0; i<4; i++ ) delete histoText[i];
119	delete canvas;
120	curr->cd();
121	}
122
123	class TESDTrackInfo
124	{
125	public:
126	TESDTrackInfo(){}
127	AliKFParticle fParticle; //* assigned KFParticle
128	Bool_t fPrimUsedFlag; //* flag shows that the particle was used for primary vertex fit
129	Bool_t fOK; //* is the track good enough
130	Int_t mcPDG; //* Monte Carlo PDG code
131	Int_t mcMotherID; //* Monte Carlo ID of its mother
132	};
133
86ad5fcb	134	void RunV0( AliESDEvent *event )
cc7cb084	135	{
	136	//* V0 finder
	137
	138	static Int_t iEvent = 0;
	139	cout<<"Event "<<iEvent++<<endl;
	140
	141	if( !gAlice ) return;
33c3c91a	142	AliRunLoader *rl = AliRunLoader::Instance();
cc7cb084	143
	144	AliStack *stack = rl->Stack();
	145	if( !stack ) return;
	146
	147	Int_t nESDTracks=event->GetNumberOfTracks();
	148	if( nESDTracks>1000 ) nESDTracks=1000;
	149
	150
	151	TESDTrackInfo ESDTrackInfo[1000]; //* parallel to ESD tracks
	152
	153	//* Fill ESDTrackInfo array
	154
	155	for (Int_t iTr=0; iTr<nESDTracks; iTr++)
	156	{
	157	TESDTrackInfo &info = ESDTrackInfo[iTr];
	158	info.fOK = 0;
	159	info.fPrimUsedFlag = 0;
	160	info.mcPDG = -1;
	161	info.mcMotherID = -1;
	162
	163	//* track quality check
	164
	165	AliESDtrack *pTrack = event->GetTrack(iTr);
	166	if( !pTrack ) continue;
	167	if (pTrack->GetKinkIndex(0)>0) continue;
	168	if ( !( pTrack->GetStatus()&AliESDtrack::kITSrefit ) ) continue;
	169	Int_t indi[12];
	170	if( pTrack->GetITSclusters(indi) <5 ) continue;
6c23ffed	171	Int_t PDG = ( pTrack->GetSigned1Pt() <0 ) ?321 :211;
cc7cb084	172
	173	//* take MC PDG
	174	{
	175	Int_t mcID = TMath::Abs(pTrack->GetLabel());
	176	TParticle * part = stack->Particle(TMath::Abs(mcID));
	177	if( part ){
	178	info.mcPDG = part->GetPdgCode();
	179	PDG = info.mcPDG;
	180	if( mcID>=0 ) info.mcMotherID = part->GetFirstMother();
	181	}
	182	}
	183
	184	//* Construct KFParticle for the track
	185
	186	info.fParticle = AliKFParticle( *pTrack, PDG );
	187	info.fOK = 1;
	188	}
	189
	190	//* Find event primary vertex
	191
	192	AliKFVertex primVtx;
	193	{
	194	const AliKFParticle * vSelected[1000]; //* Selected particles for vertex fit
	195	Int_t vIndex[1000]; //* Indices of selected particles
	196	Bool_t vFlag[1000]; //* Flags returned by the vertex finder
	197
	198	Int_t nSelected = 0;
	199	for( Int_t i = 0; i<nESDTracks; i++){
	200	if(ESDTrackInfo[i].fOK ){
	201	vSelected[nSelected] = &(ESDTrackInfo[i].fParticle);
	202	vIndex[nSelected] = i;
	203	nSelected++;
	204	}
	205	}
	206	primVtx.ConstructPrimaryVertex( vSelected, nSelected, vFlag, 3. );
	207	for( Int_t i = 0; i<nSelected; i++){
	208	if( vFlag[i] ) ESDTrackInfo[vIndex[i]].fPrimUsedFlag = 1;
	209	}
	210	if( primVtx.GetNDF() <1 ) return; //* Less then two tracks in primary vertex
	211	}
	212
	213	//* V0 finder
	214
	215	for( Int_t iTr = 0; iTr<nESDTracks; iTr++ ){ //* first daughter
	216
	217	TESDTrackInfo &info = ESDTrackInfo[iTr];
	218	if( !info.fOK ) continue;
	219
	220	for( Int_t jTr = iTr+1; jTr<nESDTracks; jTr++ ){ //* second daughter
	221	TESDTrackInfo &jnfo = ESDTrackInfo[jTr];
	222	if( !jnfo.fOK ) continue;
	223
	224	//* check for different charge
	225
	226	if( info.fParticle.GetQ() == jnfo.fParticle.GetQ() ) continue;
	227
	228	//* construct V0 mother
	229
616ffc76	230	AliKFParticle V0( info.fParticle, jnfo.fParticle );
cc7cb084	231
	232	//* check V0 Chi^2
	233
	234	if( V0.GetNDF()<1 ) continue;
	235	if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) >3. ) continue;
	236
	237	//* subtruct daughters from primary vertex
	238
	239	AliKFVertex primVtxCopy = primVtx;
	240
	241	if( info.fPrimUsedFlag ) primVtxCopy -= info.fParticle;
	242	if( jnfo.fPrimUsedFlag ) primVtxCopy -= jnfo.fParticle;
	243
	244	//* Check V0 Chi^2 deviation from primary vertex
	245
	246	if( V0.GetDeviationFromVertex( primVtxCopy ) >3. ) continue;
	247
	248	//* Add V0 to primary vertex to improve the primary vertex resolution
	249
	250	primVtxCopy += V0;
	251
	252	//* Set production vertex for V0
	253
	254	V0.SetProductionVertex( primVtxCopy );
	255
	256	//* Check chi^2 for a case
	257
0a7eb1ed	258	if( TMath::Sqrt( TMath::Abs(V0.GetChi2()/V0.GetNDF())) >3. ) continue;
cc7cb084	259
	260	//* Get V0 decay length with estimated error
	261
	262	Double_t length, sigmaLength;
	263	if( V0.GetDecayLength( length, sigmaLength ) ) continue;
	264
	265	//* Reject V0 if it decays too close to the primary vertex
	266
	267	if( length <3.*sigmaLength ) continue;
	268
	269	//* Get V0 invariant mass and plot it
	270
	271	Double_t mass, sigmaMass;
	272	if( V0.GetMass( mass, sigmaMass ) ) continue;
	273	histoMass[0]->Fill(mass);
	274
	275	//* Fill signal histograms using MC information
	276
	277	if( info.mcMotherID==jnfo.mcMotherID && info.mcMotherID>=0 ){
	278	TParticle *mother = stack->Particle(info.mcMotherID);
	279	if( mother && TMath::Abs(mother->GetPdgCode())!=21 ){
	280	histoMass[1]->Fill(mass);
	281	if( mother->GetNDaughters()==2 ){
	282	histoMass[2]->Fill(mass);
	283	}
	284	cout<<"PDG V0,pi,pj, ndaughters, mc mass, reco mass = "<<mother->GetPdgCode()<<","<<info.mcPDG<<","<<jnfo.mcPDG<<", "
	285	<< mother->GetNDaughters()<<", "<<mother->GetMass()<<", "<<mass<<endl;
	286	}
	287	}
	288	}
	289	}
	290	if( iEvent %100 == 0 \|\| (iEvent<100 && iEvent %10==0)) DrawV0();
	291	}
	292
	293
	294
9c12e1cf	295	Int_t AliKFParticleTest(Int_t n1=0,Int_t n2=1000,char *dire="/d/alice10/sma/my_v4-05-Release/pp/"){
cc7cb084	296	//* Main macro
	297
	298	// LOOP OVER SERIES OF DIRECTORIES
	299
	300	for (Int_t ifi=n1; ifi<=n2; ifi++) {
	301
	302	char nstring[5], filename[100], esdfile[100];
	303	sprintf(nstring,"%3.3d",ifi);
	304	sprintf(filename,"%s%s/galice.root",dire,nstring);
	305	sprintf(esdfile,"%s%s/AliESDs.root",dire,nstring);
	306
	307	cout <<" Opening "<<filename<<"\nand ESD "<<esdfile<<endl;
	308
0a7eb1ed	309	TFile* ftemp = TFile::Open(filename);
	310	Printf("ftemp = %p", ftemp);
	311	if (!ftemp) continue;
	312
cc7cb084	313	if (gAlice) {
33c3c91a	314	delete AliRunLoader::Instance();
cc7cb084	315	delete gAlice;
	316	gAlice=0;
	317	}
	318
	319	AliRunLoader *rl = AliRunLoader::Open(filename);
	320	if ( !rl ) {
	321	::Error("AliKFParticleTest.C","Can not open session !");
	322	continue;
	323	}
	324	if (rl->LoadgAlice()) {
	325	::Error("AliKFParticleTest.C","LoadgAlice returned error !");
	326	delete rl;
	327	continue;
	328	}
	329	if (rl->LoadHeader()) {
	330	::Error("AliKFParticleTest.C","LoadHeader returned error !");
	331	delete rl;
	332	continue;
	333	}
	334	rl->LoadKinematics();
cc7cb084	335
	336	//---------------------------------------//
	337	// //
	338	// ESD file //
	339	// //
	340	//---------------------------------------//
	341
	342	// Open file with the ESD
	343	TFile *ef=TFile::Open(esdfile);
	344
	345	//Check if the file could be opened
	346	if (!ef \|\| !ef->IsOpen()) {cerr<<"Error open AliESDs.root !\n"; continue ;}
	347
	348	//create event object
0a7eb1ed	349	AliESDEvent *event = new AliESDEvent;
cc7cb084	350
	351	//Set pointer to the esd tree in the file
	352	TTree* tree = (TTree*) ef->Get("esdTree");
	353
	354	//check if the tree exists
	355	if (!tree) {cerr<<"no ESD tree found\n"; continue;};
	356
	357	//Set pointer to the esd object in the tree
0a7eb1ed	358	event->ReadFromTree(tree);
cc7cb084	359
	360	//Number of events
	361	Int_t nevents=tree->GetEntriesFast();
	362	cout << "Number of events: " << nevents << endl;
	363
	364	StartV0();
	365	for (Int_t iev=0; iev<nevents; iev++){
	366	tree->GetEvent(iev);
	367	rl->GetEvent(iev);
0a7eb1ed	368	AliKFParticle::SetField(event->GetMagneticField());
cc7cb084	369	RunV0(event);
	370	}
	371	EndV0();
	372	delete event;
	373	ef->Close();
	374	}
	375	StartV0();
	376	return 0;
	377	}