1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 #include "AliRICHReconstructor.h" //class header
17 #include "AliRICH.h" //Reconstruct(...)
18 #include <AliRawReader.h> //Reconstruct(...)
19 #include <AliRun.h> //ConvertDigits uses gAlice
20 #include <AliESD.h> //RichAna()
21 #include <TFile.h> //RichAna()
22 #include <TMinuit.h> //Dig2Clu()
23 ClassImp(AliRICHReconstructor)
25 void AliRICHReconstructor::Reconstruct(AliRunLoader *pAL,AliRawReader* pRR)const
27 //Invoked by AliReconstruction to convert raw digits from DDL files to clusters
28 //Arguments: pAL - ALICE run loader pointer
29 // pRR - ALICE raw reader pointer
31 AliLoader *pRL=pAL->GetDetectorLoader("RICH"); AliRICH *pRich=(AliRICH*)pAL->GetAliRun()->GetDetector("RICH");//get pointers for RICH and RICH loader
33 AliRICHDigit dig; //tmp digit, raw digit will be converted to it
34 TClonesArray *pDigList=new TClonesArray("AliRICHDigit"); Int_t iDigCnt=0; //tmp list of digits for single chamber only
37 while(pRR->NextEvent()){//events loop
38 pAL->GetEvent(iEvtN++);
39 pRL->MakeTree("R"); pRich->MakeBranch("R");
41 for(Int_t iChN=1;iChN<=7;iChN++){//chambers loop
42 pRR->Select(AliRICHDigit::kRichRawId,2*iChN-2,2*iChN-1);//select only DDL files for the current chamber
44 while(pRR->ReadNextInt(w32)){//raw records loop (in selected DDL files)
45 UInt_t ddl=pRR->GetDDLID(); //returns 0,1,2 ... 13
47 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()));
48 new((*pDigList)[iDigCnt++]) AliRICHDigit(dig); //add this digit to the tmp list
50 if(iDigCnt) Dig2Clu(pDigList,pRich->Clusters(iChN));//cluster finder for the current chamber if any digits present
52 pDigList->Delete(); iDigCnt=0;//clean up list of digits for the current chamber
54 pRL->TreeR()->Fill(); //fill tree for current event
55 pRL->WriteRecPoints("OVERWRITE");//write out clusters for current event
56 pRich->ClustersReset();
58 pRL->UnloadRecPoints();
59 }//Reconstruct raw data
60 //__________________________________________________________________________________________________
61 void AliRICHReconstructor::Dig2Clu(TClonesArray *pDigList,TClonesArray *pCluList)const
63 //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.
64 //Puts all found clusters in the given clusters list.
65 //Arguments: pDigList - list of digits provided not empty
67 TMatrixF digMap(1,AliRICHParam::NpadsX(),1,AliRICHParam::NpadsY()); digMap=(Float_t)-1; //digit map for one chamber reseted to -1
68 for(Int_t iDigN=0 ; iDigN < pDigList->GetEntriesFast() ; iDigN++){ //digits loop to fill digits map
69 AliRICHDigit *pDig= (AliRICHDigit*)pDigList->At(iDigN);
70 digMap( pDig->PadX(), pDig->PadY() )=iDigN; //(padx,pady) cell takes digit number
73 AliRICHCluster clu; Int_t iCluCnt=0; //tmp cluster and cluster counter
75 for(Int_t iDigN=0;iDigN<pDigList->GetEntriesFast();iDigN++){//digits loop
76 AliRICHDigit *pDig=(AliRICHDigit*)pDigList->At(iDigN);
77 if(!(pDig=UseDig(pDig->PadX(),pDig->PadY(),pDigList,&digMap))) continue; //this digit is already taken in FormCluster(), go after next digit
78 FormCluster(&clu,pDig,pDigList,&digMap); //form cluster starting from this digit
79 TMinuit *pMinuit=clu.Solve(); //solve this cluster
81 if(pMinuit){//means cluster is solved into local maxima number of clusters, so add all of them in loop
82 Double_t x=-1,y=-1,q=-1;TString str; Double_t b1,b2,b3; Int_t iErrFlg;//tmp to withdraw resulting parameters
83 for(Int_t i=0;i<clu.Nlocmax();i++){//take resulting parameters
84 pMinuit->mnpout(3*i ,str, x, b1, b2, b3, iErrFlg);
85 pMinuit->mnpout(3*i+1,str, y, b1, b2, b3, iErrFlg);
86 pMinuit->mnpout(3*i+2,str, q, b1, b2, b3, iErrFlg);
87 clu.Set(x,y,(Int_t)q);
88 new((*pCluList)[iCluCnt++]) AliRICHCluster(clu);
91 }else//means cluster is solved as simple center of gravity cluster, add it
92 new((*pCluList)[iCluCnt++]) AliRICHCluster(clu);
93 clu.Reset();//make current cluster empty
96 //__________________________________________________________________________________________________
97 void AliRICHReconstructor::FormCluster(AliRICHCluster *pClu,AliRICHDigit *pDig,TClonesArray *pDigList,TMatrixF *pDigMap)const
99 //Forms the initial cluster as a sum of all adjascent digits. Starts from the given digit
100 //then calls itself recursevly for all neighbours.
101 //Arguments: pClu - pointer to cluster being formed
103 pClu->AddDigit(pDig);//take this digit in cluster and mark it as taken
107 Int_t iNnei=AliRICHParam::PadNeighbours(pDig->PadX(),pDig->PadY(),x,y);//returns in x,y all possible neighbours of the given one
108 for (Int_t i=0;i<iNnei;i++)
109 if((pDig=UseDig(x[i],y[i],pDigList,pDigMap))) FormCluster(pClu,pDig,pDigList,pDigMap);
111 //__________________________________________________________________________________________________
112 void AliRICHReconstructor::RichAna(Int_t nNevMax,Bool_t askPatRec)
114 TFile *pFile=TFile::Open("AliESDs.root","read");
115 if(!pFile || !pFile->IsOpen()) {return;} //open AliESDs.root
116 TTree *pTree = (TTree*) pFile->Get("esdTree");
117 if(!pTree){return;} //get ESD tree
119 AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
121 AliMagF * magf = gAlice->Field(); AliTracker::SetFieldMap(magf,kTRUE); // AliTracker::SetFieldMap(magf);
123 TFile *pFileRA = new TFile("$(HOME)/RichAna.root","RECREATE","RICH Pattern Recognition");
124 TNtupleD *hn = new TNtupleD("hn","ntuple","Pmod:Charge:TrackTheta:TrackPhi:MinX:MinY:ThetaCerenkov:NPhotons:ChargeMIP:Chamber:TOF:LengthTOF:"
125 "prob1:prob2:prob3:ErrPar1:ErrPar2:ErrPar3:Th1:Th2:Th3:nPhotBKG");
126 if(nNevMax==0) nNevMax=999999;
127 if(pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents()<nNevMax) nNevMax = pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();
128 AliESD *pESD=new AliESD; pTree->SetBranchAddress("ESD", &pESD);
129 for(Int_t iEvtN=0;iEvtN<nNevMax;iEvtN++) {
130 pTree->GetEvent(iEvtN);
131 pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
132 pRich->GetLoader()->LoadRecPoints();
133 pRich->GetLoader()->TreeR()->GetEntry(0);
134 //Pattern recognition started
135 if(pESD->GetNumberOfTracks()) {
136 Int_t iNtracks=pESD->GetNumberOfTracks();
137 AliInfoClass(Form("Start with %i tracks",iNtracks));
138 for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
139 if(iTrackN%100==0)AliInfoClass(Form("Track %i to be processed",iTrackN));
140 AliRICHTracker *pTrRich = new AliRICHTracker();
141 if(askPatRec==kTRUE) pTrRich->RecWithESD(pESD,pRich,iTrackN);
142 AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
145 pTrack->GetRICHdxdy(dx,dy);
146 hnvec[0]=pTrack->GetP();
147 hnvec[1]=pTrack->GetSign();
149 pTrack->GetRICHthetaPhi(hnvec[2],hnvec[3]);
150 pTrack->GetRICHdxdy(hnvec[4],hnvec[5]);
151 hnvec[6]=pTrack->GetRICHsignal();
152 hnvec[7]=pTrack->GetRICHnclusters();
153 hnvec[9]=pTrack->GetRICHcluster()/1000000;
154 hnvec[8]=pTrack->GetRICHcluster()-hnvec[9]*1000000;
155 hnvec[10]=pTrack->GetTOFsignal();
156 hnvec[11]=pTrack->GetIntegratedLength();
158 pTrack->GetRICHpid(prob);
159 hnvec[12]=prob[0]+prob[1]+prob[2];
162 hnvec[15]=pTrRich->fErrPar[2];
163 hnvec[16]=pTrRich->fErrPar[3];
164 hnvec[17]=pTrRich->fErrPar[4];
165 for(Int_t i=0;i<3;i++) {
166 Double_t mass = AliRICHParam::fgMass[i+2];
167 Double_t refIndex=AliRICHParam::RefIdxC6F14(AliRICHParam::MeanCkovEnergy());
168 Double_t cosThetaTh = TMath::Sqrt(mass*mass+pTrack->GetP()*pTrack->GetP())/(refIndex*pTrack->GetP());
170 if(cosThetaTh>=1) continue;
171 hnvec[18+i]= TMath::ACos(cosThetaTh);
173 hnvec[21]=pTrRich->fnPhotBKG;
177 AliInfoClass(Form("Pattern Recognition done for event %i",iEvtN));
179 pFileRA->Write();pFileRA->Close();// close RichAna.root
180 delete pESD; pFile->Close();//close AliESDs.root
182 //__________________________________________________________________________________________________
183 void AliRICHReconstructor::CheckPR()
185 //Pattern recognition with stack particles
186 TFile *pFile = new TFile("$(HOME)/RPR.root","RECREATE","RICH Pattern Recognition");
187 TNtupleD *hn = new TNtupleD("hn","ntuple","Pmod:Charge:TrackTheta:TrackPhi:TrackX:TrackY:MinX:MinY:ChargeMIP:ThetaCerenkov:NPhotons:MipIndex:Chamber:Particle");
188 AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
189 AliMagF * magf = gAlice->Field();
190 AliTracker::SetFieldMap(magf,kTRUE);
191 for(Int_t iEvtN=0;iEvtN<pRich->GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvtN++) {
192 pRich->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
193 AliRICHTracker *tr = new AliRICHTracker();
194 tr->RecWithStack(hn);
195 AliInfoClass(Form("Pattern Recognition done for event %i \b",iEvtN));
197 pFile->Write();pFile->Close();