[u/mrichter/AliRoot.git] / RICH / AliRICHReconstructor.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

#include "AliRICHReconstructor.h" //class header
#include "AliRICH.h"              //Reconstruct(...) 
#include <AliRawReader.h>         //Reconstruct(...)
#include <AliRun.h>               //ConvertDigits uses gAlice
#include <AliESD.h>            //RichAna()
#include <AliStack.h>          //RichAna()
#include <TFile.h>             //RichAna()
#include <TMinuit.h>              //Dig2Clu()
#include <TParticle.h>            //TParticle()
ClassImp(AliRICHReconstructor)

void AliRICHReconstructor::Reconstruct(AliRunLoader *pAL,AliRawReader* pRR)const
{
//Invoked  by AliReconstruction to convert raw digits from DDL files to clusters
//Arguments: pAL - ALICE run loader pointer
//           pRR - ALICE raw reader pointer  
//  Returns: none    
  AliLoader *pRL=pAL->GetDetectorLoader("RICH");  AliRICH *pRich=(AliRICH*)pAL->GetAliRun()->GetDetector("RICH");//get pointers for RICH and RICH loader
  
  AliRICHDigit dig; //tmp digit, raw digit will be converted to it
  TClonesArray *pDigList=new TClonesArray("AliRICHDigit"); Int_t iDigCnt=0; //tmp list of digits for single chamber only
  
  Int_t iEvtN=0;
  while(pRR->NextEvent()){//events loop
    pAL->GetEvent(iEvtN++);
    pRL->MakeTree("R");  pRich->MakeBranch("R");
    
    for(Int_t iChN=1;iChN<=7;iChN++){//chambers loop
      pRR->Select(AliRICHDigit::kRichRawId,2*iChN-2,2*iChN-1);//select only DDL files for the current chamber      
      UInt_t w32=0;
      while(pRR->ReadNextInt(w32)){//raw records loop (in selected DDL files)
        UInt_t ddl=pRR->GetDDLID(); //returns 0,1,2 ... 13
        dig.Raw2Dig(ddl,w32);  
        AliDebug(1,Form("Ch=%i DDL=%i raw=0x%x digit=(%3i,%3i,%3i,%3i) Q=%5.2f",iChN,ddl,w32,dig.Chamber(),dig.Sector(),dig.PadX(),dig.PadY(),dig.Qdc()));
        new((*pDigList)[iDigCnt++]) AliRICHDigit(dig); //add this digit to the tmp list
      }//raw records loop
      if(iDigCnt) Dig2Clu(pDigList,pRich->Clusters(iChN));//cluster finder for the current chamber if any digits present
      pRR->Reset();        
      pDigList->Delete();  iDigCnt=0;//clean up list of digits for the current chamber
    }//chambers loop
    pRL->TreeR()->Fill();            //fill tree for current event
    pRL->WriteRecPoints("OVERWRITE");//write out clusters for current event
    pRich->ClustersReset();
  }//events loop  
  pRL->UnloadRecPoints();  
}//Reconstruct raw data
//__________________________________________________________________________________________________
void AliRICHReconstructor::Dig2Clu(TClonesArray *pDigList,TClonesArray *pCluList)const
{
//Finds all clusters for a given digits list provided not empty. Currently digits list provided is a list of all digits for a single chamber.
//Puts all found clusters in the given clusters list. 
//Arguments: pDigList - list of digits provided not empty
//  Returns: none    
  TMatrixF digMap(1,AliRICHParam::NpadsX(),1,AliRICHParam::NpadsY());  digMap=(Float_t)-1; //digit map for one chamber reseted to -1
  for(Int_t iDigN=0 ; iDigN < pDigList->GetEntriesFast() ; iDigN++){ //digits loop to fill digits map
    AliRICHDigit *pDig= (AliRICHDigit*)pDigList->At(iDigN);
    digMap( pDig->PadX(), pDig->PadY() )=iDigN;                     //(padx,pady) cell takes digit number
  }
  
  AliRICHCluster clu; Int_t iCluCnt=0; //tmp cluster and cluster counter
  
  for(Int_t iDigN=0;iDigN<pDigList->GetEntriesFast();iDigN++){//digits loop    
    AliRICHDigit *pDig=(AliRICHDigit*)pDigList->At(iDigN);
    if(!(pDig=UseDig(pDig->PadX(),pDig->PadY(),pDigList,&digMap))) continue;  //this digit is already taken in FormCluster(), go after next digit
    FormCluster(&clu,pDig,pDigList,&digMap);                            //form cluster starting from this digit
    TMinuit *pMinuit=clu.Solve();                              //solve this cluster
    
    if(pMinuit){//means cluster is solved into local maxima number of clusters, so add all of them in loop
      Double_t x=-1,y=-1,q=-1;TString str; Double_t b1,b2,b3; Int_t iErrFlg;//tmp to withdraw resulting parameters
      for(Int_t i=0;i<clu.Nlocmax();i++){//take resulting parameters 
        pMinuit->mnpout(3*i  ,str, x, b1, b2, b3, iErrFlg);
        pMinuit->mnpout(3*i+1,str, y, b1, b2, b3, iErrFlg);
        pMinuit->mnpout(3*i+2,str, q, b1, b2, b3, iErrFlg);
        clu.Set(x,y,(Int_t)q);
        new((*pCluList)[iCluCnt++]) AliRICHCluster(clu);
      }
      delete pMinuit;
    }else//means cluster is solved as simple center of gravity cluster, add it
      new((*pCluList)[iCluCnt++]) AliRICHCluster(clu); 
    clu.Reset();//make current cluster empty
  }//digits loop
}//Dig2Clu()
//__________________________________________________________________________________________________
void  AliRICHReconstructor::FormCluster(AliRICHCluster *pClu,AliRICHDigit *pDig,TClonesArray *pDigList,TMatrixF *pDigMap)const
{
//Forms the initial cluster as a sum of all adjascent digits. Starts from the given digit
//then calls itself recursevly  for all neighbours.
//Arguments: pClu - pointer to cluster being formed
//  Returns: none
  pClu->AddDigit(pDig);//take this digit in cluster and mark it as taken

  Int_t x[4],y[4];
  
  Int_t iNnei=AliRICHParam::PadNeighbours(pDig->PadX(),pDig->PadY(),x,y);//returns in x,y all possible neighbours of the given one
  for (Int_t i=0;i<iNnei;i++)
    if((pDig=UseDig(x[i],y[i],pDigList,pDigMap))) FormCluster(pClu,pDig,pDigList,pDigMap);    
}//FormCluster()
//__________________________________________________________________________________________________
void AliRICHReconstructor::CheckPR()
{
//Pattern recognition with stack particles
  TFile *pFile = new TFile("$(HOME)/RPR.root","RECREATE","RICH Pattern Recognition");
  TNtupleD *hn = new TNtupleD("hn","ntuple","Pmod:Charge:TrackTheta:TrackPhi:TrackX:TrackY:MinX:MinY:ChargeMIP:ThetaCerenkov:NPhotons:MipIndex:Chamber:Particle");
//  printf("\n\n");
//  printf("Pattern Recognition done for event %5i",0);
  AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
  AliMagF * magf = gAlice->Field();
  AliTracker::SetFieldMap(magf,kTRUE);
  for(Int_t iEvtN=0;iEvtN<pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvtN++) {
    pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
    AliRICHTracker *tr = new AliRICHTracker();
    tr->RecWithStack(hn);
    AliInfoClass(Form("Pattern Recognition done for event %i \b",iEvtN));
//    printf("\b\b\b\b\b%5i",iEvtN+1);
  }
  printf("\n\n");
  pFile->Write();pFile->Close();
}
//__________________________________________________________________________________________________
void AliRICHReconstructor::RichAna(Int_t iNevMin,Int_t iNevMax,Bool_t askPatRec)
{
  TFile *pFile=TFile::Open("AliESDs.root","read");
  if(!pFile || !pFile->IsOpen()) {AliInfoClass("ESD file not open.");return;}      //open AliESDs.root                                                                    
  TTree *pTree = (TTree*) pFile->Get("esdTree");
  if(!pTree){AliInfoClass("ESD not found.");return;}                               //get ESD tree  
  AliInfoClass("ESD found. Go ahead!");
  
  AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));

  AliMagF * magf = gAlice->Field();
  AliTracker::SetFieldMap(magf,kTRUE);
  pRich->GetLoader()->GetRunLoader()->LoadHeader();
  pRich->GetLoader()->GetRunLoader()->LoadKinematics();
  TString var1 = "Pmod:Charge:TrackTheta:TrackPhi:MinX:MinY:ThetaCerenkov:NPhotons:";
  TString var2 = "ChargeMIP:Chamber:TOF:LengthTOF:prob1:prob2:prob3:";
  TString var3 = "ErrPar1:ErrPar2:ErrPar3:Th1:Th2:Th3:nPhotBKG:pdgCode";
  TString varList = var1+var2+var3;
  
  Double_t dx,dy;
  Double_t hnvec[30];

//  TFile *pFileRA = new TFile("$(HOME)/RichAna.root","RECREATE","RICH Pattern Recognition");
  TFile *pFileRA = new TFile("./RichAna.root","RECREATE","RICH Pattern Recognition");
  TNtupleD *hn = new TNtupleD("hn","ntuple",varList);
  if(iNevMin<0) iNevMin=0;
  if(iNevMin>iNevMax) {iNevMin=0;iNevMax=0;}  
  if(iNevMax==0) iNevMax=999999;
  if(pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents()<iNevMax) iNevMax = pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();
  AliESD *pESD=new AliESD;  pTree->SetBranchAddress("ESD", &pESD);
  for(Int_t iEvtN=iNevMin;iEvtN<iNevMax;iEvtN++) {
    pTree->GetEvent(iEvtN);
    pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
    pRich->GetLoader()->LoadRecPoints();
    pRich->GetLoader()->TreeR()->GetEntry(0);
    AliStack *pStack = pRich->GetLoader()->GetRunLoader()->Stack();
//Pattern recognition started
    if(pESD->GetNumberOfTracks()) {
      Int_t iNtracks=pESD->GetNumberOfTracks();
      AliInfoClass(Form("Start with %i tracks",iNtracks));
      for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
        if(iTrackN%100==0)AliInfoClass(Form("Track %i to be processed",iTrackN));
        AliRICHTracker *pTrRich = new AliRICHTracker();
        if(askPatRec==kTRUE) pTrRich->RecWithESD(pESD,pRich,iTrackN);
        AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
        
        Int_t lab=TMath::Abs(pTrack->GetLabel());
        TParticle *pPart=pStack->Particle(lab);
        Int_t code=pPart->GetPdgCode();

        pTrack->GetRICHdxdy(dx,dy);
        hnvec[0]=pTrack->GetP();
        hnvec[1]=pTrack->GetSign();
//  cout << " Track momentum " << pTrack->GetP() << " charge " << pTrack->GetSign() << endl;
  
        pTrack->GetRICHthetaPhi(hnvec[2],hnvec[3]);
        pTrack->GetRICHdxdy(hnvec[4],hnvec[5]);
        hnvec[6]=pTrack->GetRICHsignal();
        hnvec[7]=pTrack->GetRICHnclusters();
        hnvec[9]=pTrack->GetRICHcluster()/1000000;
        hnvec[8]=pTrack->GetRICHcluster()-hnvec[9]*1000000;
        hnvec[10]=pTrack->GetTOFsignal();
        hnvec[11]=pTrack->GetIntegratedLength();
        Double_t prob[5];
        pTrack->GetRICHpid(prob);
        hnvec[12]=prob[0]+prob[1]+prob[2];
        hnvec[13]=prob[3];
        hnvec[14]=prob[4];
        hnvec[15]=pTrRich->fErrPar[2];
        hnvec[16]=pTrRich->fErrPar[3];
        hnvec[17]=pTrRich->fErrPar[4];
        for(Int_t i=0;i<3;i++) {
          Double_t mass = AliRICHParam::fgMass[i+2];
          Double_t refIndex=AliRICHParam::RefIdxC6F14(AliRICHParam::MeanCkovEnergy());
          Double_t cosThetaTh = TMath::Sqrt(mass*mass+pTrack->GetP()*pTrack->GetP())/(refIndex*pTrack->GetP());
          hnvec[18+i]=0;
          if(cosThetaTh>=1) continue;
          hnvec[18+i]= TMath::ACos(cosThetaTh);
        }
        if(askPatRec==kTRUE) hnvec[21]=pTrRich->fnPhotBKG; else hnvec[21]=0;
        hnvec[22]=code;
        hn->Fill(hnvec);
      }
    }
    AliInfoClass(Form("Pattern Recognition done for event %i",iEvtN));
  }
  pFileRA->Write();pFileRA->Close();// close RichAna.root
  delete pESD;  pFile->Close();//close AliESDs.root
}
//__________________________________________________________________________________________________
Commit	Line	Data
121a60bd	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
0422a446	16	#include "AliRICHReconstructor.h" //class header
122efc54	17	#include "AliRICH.h" //Reconstruct(...)
0422a446	18	#include <AliRawReader.h> //Reconstruct(...)
0422a446	19	#include <AliRun.h> //ConvertDigits uses gAlice
122efc54	20	#include <AliESD.h> //RichAna()
01c81d9d	21	#include <AliStack.h> //RichAna()
122efc54	22	#include <TFile.h> //RichAna()
0422a446	23	#include <TMinuit.h> //Dig2Clu()
01c81d9d	24	#include <TParticle.h> //TParticle()
121a60bd	25	ClassImp(AliRICHReconstructor)
121a60bd	26
0422a446	27	void AliRICHReconstructor::Reconstruct(AliRunLoader pAL,AliRawReader pRR)const
0422a446	28	{
122efc54	29	//Invoked by AliReconstruction to convert raw digits from DDL files to clusters
	30	//Arguments: pAL - ALICE run loader pointer
	31	// pRR - ALICE raw reader pointer
	32	// Returns: none
	33	AliLoader pRL=pAL->GetDetectorLoader("RICH"); AliRICH pRich=(AliRICH*)pAL->GetAliRun()->GetDetector("RICH");//get pointers for RICH and RICH loader
0422a446	34
	35	AliRICHDigit dig; //tmp digit, raw digit will be converted to it
	36	TClonesArray *pDigList=new TClonesArray("AliRICHDigit"); Int_t iDigCnt=0; //tmp list of digits for single chamber only
	37
	38	Int_t iEvtN=0;
	39	while(pRR->NextEvent()){//events loop
	40	pAL->GetEvent(iEvtN++);
	41	pRL->MakeTree("R"); pRich->MakeBranch("R");
	42
	43	for(Int_t iChN=1;iChN<=7;iChN++){//chambers loop
	44	pRR->Select(AliRICHDigit::kRichRawId,2iChN-2,2iChN-1);//select only DDL files for the current chamber
	45	UInt_t w32=0;
	46	while(pRR->ReadNextInt(w32)){//raw records loop (in selected DDL files)
	47	UInt_t ddl=pRR->GetDDLID(); //returns 0,1,2 ... 13
	48	dig.Raw2Dig(ddl,w32);
	49	AliDebug(1,Form("Ch=%i DDL=%i raw=0x%x digit=(%3i,%3i,%3i,%3i) Q=%5.2f",iChN,ddl,w32,dig.Chamber(),dig.Sector(),dig.PadX(),dig.PadY(),dig.Qdc()));
	50	new((*pDigList)[iDigCnt++]) AliRICHDigit(dig); //add this digit to the tmp list
	51	}//raw records loop
	52	if(iDigCnt) Dig2Clu(pDigList,pRich->Clusters(iChN));//cluster finder for the current chamber if any digits present
	53	pRR->Reset();
	54	pDigList->Delete(); iDigCnt=0;//clean up list of digits for the current chamber
	55	}//chambers loop
	56	pRL->TreeR()->Fill(); //fill tree for current event
	57	pRL->WriteRecPoints("OVERWRITE");//write out clusters for current event
	58	pRich->ClustersReset();
	59	}//events loop
	60	pRL->UnloadRecPoints();
	61	}//Reconstruct raw data
	62	//__________________________________________________________________________________________________
	63	void AliRICHReconstructor::Dig2Clu(TClonesArray pDigList,TClonesArray pCluList)const
	64	{
	65	//Finds all clusters for a given digits list provided not empty. Currently digits list provided is a list of all digits for a single chamber.
	66	//Puts all found clusters in the given clusters list.
	67	//Arguments: pDigList - list of digits provided not empty
	68	// Returns: none
	69	TMatrixF digMap(1,AliRICHParam::NpadsX(),1,AliRICHParam::NpadsY()); digMap=(Float_t)-1; //digit map for one chamber reseted to -1
	70	for(Int_t iDigN=0 ; iDigN < pDigList->GetEntriesFast() ; iDigN++){ //digits loop to fill digits map
	71	AliRICHDigit pDig= (AliRICHDigit)pDigList->At(iDigN);
	72	digMap( pDig->PadX(), pDig->PadY() )=iDigN; //(padx,pady) cell takes digit number
	73	}
	74
	75	AliRICHCluster clu; Int_t iCluCnt=0; //tmp cluster and cluster counter
	76
	77	for(Int_t iDigN=0;iDigN<pDigList->GetEntriesFast();iDigN++){//digits loop
	78	AliRICHDigit pDig=(AliRICHDigit)pDigList->At(iDigN);
	79	if(!(pDig=UseDig(pDig->PadX(),pDig->PadY(),pDigList,&digMap))) continue; //this digit is already taken in FormCluster(), go after next digit
	80	FormCluster(&clu,pDig,pDigList,&digMap); //form cluster starting from this digit
	81	TMinuit *pMinuit=clu.Solve(); //solve this cluster
	82
	83	if(pMinuit){//means cluster is solved into local maxima number of clusters, so add all of them in loop
	84	Double_t x=-1,y=-1,q=-1;TString str; Double_t b1,b2,b3; Int_t iErrFlg;//tmp to withdraw resulting parameters
	85	for(Int_t i=0;i<clu.Nlocmax();i++){//take resulting parameters
	86	pMinuit->mnpout(3*i ,str, x, b1, b2, b3, iErrFlg);
	87	pMinuit->mnpout(3*i+1,str, y, b1, b2, b3, iErrFlg);
	88	pMinuit->mnpout(3*i+2,str, q, b1, b2, b3, iErrFlg);
	89	clu.Set(x,y,(Int_t)q);
	90	new((*pCluList)[iCluCnt++]) AliRICHCluster(clu);
	91	}
	92	delete pMinuit;
	93	}else//means cluster is solved as simple center of gravity cluster, add it
	94	new((*pCluList)[iCluCnt++]) AliRICHCluster(clu);
	95	clu.Reset();//make current cluster empty
	96	}//digits loop
	97	}//Dig2Clu()
98	//__________________________________________________________________________________________________
99	void AliRICHReconstructor::FormCluster(AliRICHCluster pClu,AliRICHDigit pDig,TClonesArray pDigList,TMatrixF pDigMap)const
100	{
101	//Forms the initial cluster as a sum of all adjascent digits. Starts from the given digit
102	//then calls itself recursevly for all neighbours.
103	//Arguments: pClu - pointer to cluster being formed
104	// Returns: none
105	pClu->AddDigit(pDig);//take this digit in cluster and mark it as taken
106
107	Int_t x[4],y[4];
108
109	Int_t iNnei=AliRICHParam::PadNeighbours(pDig->PadX(),pDig->PadY(),x,y);//returns in x,y all possible neighbours of the given one
110	for (Int_t i=0;i<iNnei;i++)
111	if((pDig=UseDig(x[i],y[i],pDigList,pDigMap))) FormCluster(pClu,pDig,pDigList,pDigMap);
112	}//FormCluster()
122efc54	113	//__________________________________________________________________________________________________
01c81d9d	114	void AliRICHReconstructor::CheckPR()
	115	{
	116	//Pattern recognition with stack particles
	117	TFile *pFile = new TFile("$(HOME)/RPR.root","RECREATE","RICH Pattern Recognition");
	118	TNtupleD *hn = new TNtupleD("hn","ntuple","Pmod:Charge:TrackTheta:TrackPhi:TrackX:TrackY:MinX:MinY:ChargeMIP:ThetaCerenkov:NPhotons:MipIndex:Chamber:Particle");
	119	// printf("\n\n");
	120	// printf("Pattern Recognition done for event %5i",0);
	121	AliRICH pRich=((AliRICH)gAlice->GetDetector("RICH"));
	122	AliMagF * magf = gAlice->Field();
	123	AliTracker::SetFieldMap(magf,kTRUE);
	124	for(Int_t iEvtN=0;iEvtN<pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvtN++) {
	125	pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
	126	AliRICHTracker *tr = new AliRICHTracker();
	127	tr->RecWithStack(hn);
	128	AliInfoClass(Form("Pattern Recognition done for event %i \b",iEvtN));
	129	// printf("\b\b\b\b\b%5i",iEvtN+1);
	130	}
	131	printf("\n\n");
	132	pFile->Write();pFile->Close();
	133	}
	134	//__________________________________________________________________________________________________
	135	void AliRICHReconstructor::RichAna(Int_t iNevMin,Int_t iNevMax,Bool_t askPatRec)
122efc54	136	{
122efc54	137	TFile *pFile=TFile::Open("AliESDs.root","read");
01c81d9d	138	if(!pFile \|\| !pFile->IsOpen()) {AliInfoClass("ESD file not open.");return;} //open AliESDs.root
122efc54	139	TTree pTree = (TTree) pFile->Get("esdTree");
01c81d9d	140	if(!pTree){AliInfoClass("ESD not found.");return;} //get ESD tree
	141	AliInfoClass("ESD found. Go ahead!");
	142
122efc54	143	AliRICH pRich=((AliRICH)gAlice->GetDetector("RICH"));
01c81d9d	144
	145	AliMagF * magf = gAlice->Field();
	146	AliTracker::SetFieldMap(magf,kTRUE);
	147	pRich->GetLoader()->GetRunLoader()->LoadHeader();
	148	pRich->GetLoader()->GetRunLoader()->LoadKinematics();
	149	TString var1 = "Pmod:Charge:TrackTheta:TrackPhi:MinX:MinY:ThetaCerenkov:NPhotons:";
	150	TString var2 = "ChargeMIP:Chamber:TOF:LengthTOF:prob1:prob2:prob3:";
	151	TString var3 = "ErrPar1:ErrPar2:ErrPar3:Th1:Th2:Th3:nPhotBKG:pdgCode";
	152	TString varList = var1+var2+var3;
122efc54	153
01c81d9d	154	Double_t dx,dy;
01c81d9d	155	Double_t hnvec[30];
122efc54	156
01c81d9d	157	// TFile *pFileRA = new TFile("$(HOME)/RichAna.root","RECREATE","RICH Pattern Recognition");
	158	TFile *pFileRA = new TFile("./RichAna.root","RECREATE","RICH Pattern Recognition");
	159	TNtupleD *hn = new TNtupleD("hn","ntuple",varList);
	160	if(iNevMin<0) iNevMin=0;
	161	if(iNevMin>iNevMax) {iNevMin=0;iNevMax=0;}
	162	if(iNevMax==0) iNevMax=999999;
	163	if(pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents()<iNevMax) iNevMax = pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();
122efc54	164	AliESD *pESD=new AliESD; pTree->SetBranchAddress("ESD", &pESD);
01c81d9d	165	for(Int_t iEvtN=iNevMin;iEvtN<iNevMax;iEvtN++) {
122efc54	166	pTree->GetEvent(iEvtN);
	167	pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
	168	pRich->GetLoader()->LoadRecPoints();
	169	pRich->GetLoader()->TreeR()->GetEntry(0);
01c81d9d	170	AliStack *pStack = pRich->GetLoader()->GetRunLoader()->Stack();
122efc54	171	//Pattern recognition started
	172	if(pESD->GetNumberOfTracks()) {
	173	Int_t iNtracks=pESD->GetNumberOfTracks();
	174	AliInfoClass(Form("Start with %i tracks",iNtracks));
	175	for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
	176	if(iTrackN%100==0)AliInfoClass(Form("Track %i to be processed",iTrackN));
	177	AliRICHTracker *pTrRich = new AliRICHTracker();
	178	if(askPatRec==kTRUE) pTrRich->RecWithESD(pESD,pRich,iTrackN);
	179	AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
01c81d9d	180
	181	Int_t lab=TMath::Abs(pTrack->GetLabel());
	182	TParticle *pPart=pStack->Particle(lab);
	183	Int_t code=pPart->GetPdgCode();
	184
122efc54	185	pTrack->GetRICHdxdy(dx,dy);
	186	hnvec[0]=pTrack->GetP();
	187	hnvec[1]=pTrack->GetSign();
01c81d9d	188	// cout << " Track momentum " << pTrack->GetP() << " charge " << pTrack->GetSign() << endl;
122efc54	189
	190	pTrack->GetRICHthetaPhi(hnvec[2],hnvec[3]);
	191	pTrack->GetRICHdxdy(hnvec[4],hnvec[5]);
	192	hnvec[6]=pTrack->GetRICHsignal();
	193	hnvec[7]=pTrack->GetRICHnclusters();
	194	hnvec[9]=pTrack->GetRICHcluster()/1000000;
	195	hnvec[8]=pTrack->GetRICHcluster()-hnvec[9]*1000000;
	196	hnvec[10]=pTrack->GetTOFsignal();
	197	hnvec[11]=pTrack->GetIntegratedLength();
	198	Double_t prob[5];
	199	pTrack->GetRICHpid(prob);
	200	hnvec[12]=prob[0]+prob[1]+prob[2];
	201	hnvec[13]=prob[3];
	202	hnvec[14]=prob[4];
	203	hnvec[15]=pTrRich->fErrPar[2];
	204	hnvec[16]=pTrRich->fErrPar[3];
	205	hnvec[17]=pTrRich->fErrPar[4];
	206	for(Int_t i=0;i<3;i++) {
	207	Double_t mass = AliRICHParam::fgMass[i+2];
	208	Double_t refIndex=AliRICHParam::RefIdxC6F14(AliRICHParam::MeanCkovEnergy());
	209	Double_t cosThetaTh = TMath::Sqrt(massmass+pTrack->GetP()pTrack->GetP())/(refIndex*pTrack->GetP());
	210	hnvec[18+i]=0;
	211	if(cosThetaTh>=1) continue;
	212	hnvec[18+i]= TMath::ACos(cosThetaTh);
	213	}
01c81d9d	214	if(askPatRec==kTRUE) hnvec[21]=pTrRich->fnPhotBKG; else hnvec[21]=0;
01c81d9d	215	hnvec[22]=code;
122efc54	216	hn->Fill(hnvec);
	217	}
	218	}
	219	AliInfoClass(Form("Pattern Recognition done for event %i",iEvtN));
	220	}
	221	pFileRA->Write();pFileRA->Close();// close RichAna.root
	222	delete pESD; pFile->Close();//close AliESDs.root
122efc54	223	}
01c81d9d	224	//__________________________________________________________________________________________________