New optical properties in. New AliRICHreco.Small mods in AliRICHDisplFast
[u/mrichter/AliRoot.git] / RICH / AliRICHClusterFinder.cxx
CommitLineData
8265fa96 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
8265fa96 16
17#include "AliRICHClusterFinder.h"
543d5224 18#include "AliRICHMap.h"
9d6f9427 19#include <TMinuit.h>
ed3ceb24 20#include <TParticle.h>
9d6f9427 21#include <TVector3.h>
cbaf35fb 22#include <AliLoader.h>
ed3ceb24 23#include <AliStack.h>
cbaf35fb 24#include <AliRun.h>
8265fa96 25
9d6f9427 26void RICHMinMathieson(Int_t &npar, Double_t *gin, Double_t &chi2, Double_t *par, Int_t iflag);
8265fa96 27
28ClassImp(AliRICHClusterFinder)
543d5224 29//__________________________________________________________________________________________________
30AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH)
31{//main ctor
32 Info("main ctor","Start.");
33
cbaf35fb 34 fRICH = pRICH;
543d5224 35
cbaf35fb 36 fHitMap = 0;
543d5224 37
543d5224 38}//main ctor
39//__________________________________________________________________________________________________
9d6f9427 40void AliRICHClusterFinder::FindLocalMaxima()
543d5224 41{// Split the cluster according to the number of maxima inside
9d6f9427 42 Info("FindLocalMaxima","Start.");
cbaf35fb 43 Int_t Nlocal = 0;
44 Int_t localX[100],localY[100];
9d6f9427 45 for(Int_t iDig1=0;iDig1<fCurrentCluster.Size();iDig1++) {
cbaf35fb 46 Int_t iNotMax = 0;
9d6f9427 47 AliRICHdigit *pDig1 = (AliRICHdigit *)fCurrentCluster.Digits()->At(iDig1);
cbaf35fb 48 Int_t padX1 = pDig1->X();
49 Int_t padY1 = pDig1->Y();
50 Double_t padQ1 = pDig1->Q();
9d6f9427 51 for(Int_t iDig2=0;iDig2<fCurrentCluster.Size();iDig2++) {
52 AliRICHdigit *pDig2 = (AliRICHdigit *)fCurrentCluster.Digits()->At(iDig2);
cbaf35fb 53 Int_t padX2 = pDig2->X();
54 Int_t padY2 = pDig2->Y();
55 Double_t padQ2 = pDig2->Q();
56 if(iDig1==iDig2) continue;
57 Int_t diffx = TMath::Sign(padX1-padX2,1);
58 Int_t diffy = TMath::Sign(padY1-padY2,1);
59 if((diffx+diffy)<=1) {
60 if(padQ2>padQ1) iNotMax++;
61 }
8265fa96 62 }
cbaf35fb 63 if(iNotMax==0) {
64 localX[Nlocal] = padX1;
65 localY[Nlocal] = padY1;
66 Nlocal++;
8265fa96 67 }
cbaf35fb 68 }
69}//FindLocalMaxima()
543d5224 70//__________________________________________________________________________________________________
71void AliRICHClusterFinder::Exec()
237c933d 72{
543d5224 73 Info("Exec","Start.");
74
cbaf35fb 75
543d5224 76 Rich()->GetLoader()->LoadDigits();
77
ed3ceb24 78 Rich()->GetLoader()->GetRunLoader()->LoadHeader();
79 Rich()->GetLoader()->GetRunLoader()->LoadKinematics();
80
543d5224 81 for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop
82 gAlice->GetRunLoader()->GetEvent(iEventN);
237c933d 83
543d5224 84 Rich()->GetLoader()->MakeTree("R"); Rich()->MakeBranch("R");
cbaf35fb 85 Rich()->ResetDigits(); Rich()->ResetClusters();
543d5224 86
87 Rich()->GetLoader()->TreeD()->GetEntry(0);
cbaf35fb 88 for(Int_t iChamber=1;iChamber<=kNCH;iChamber++){//chambers loop
ed3ceb24 89 FindClusters(iChamber);
543d5224 90 }//chambers loop
543d5224 91 Rich()->GetLoader()->TreeR()->Fill();
92 Rich()->GetLoader()->WriteRecPoints("OVERWRITE");
93 }//events loop
94 Rich()->GetLoader()->UnloadDigits(); Rich()->GetLoader()->UnloadRecPoints();
cbaf35fb 95 Rich()->ResetDigits(); Rich()->ResetClusters();
ed3ceb24 96
97 Rich()->GetLoader()->GetRunLoader()->UnloadHeader();
98 Rich()->GetLoader()->GetRunLoader()->UnloadKinematics();
99
543d5224 100 Info("Exec","Stop.");
101}//Exec()
102//__________________________________________________________________________________________________
ed3ceb24 103void AliRICHClusterFinder::FindClusters(Int_t iChamber)
104{
105 //finds neighbours and fill the tree with raw clusters
cbaf35fb 106 Int_t nDigits=Rich()->Digits(iChamber)->GetEntriesFast();
ed3ceb24 107 Info("FindClusters","Start for Chamber %i with %i digits.",iChamber,nDigits);
cbaf35fb 108 if(nDigits==0)return;
109
eaf390d9 110 fHitMap=new AliRICHMap(Rich()->Digits(iChamber));//create digit map for the given chamber
cbaf35fb 111
cbaf35fb 112 for(Int_t iDig=0;iDig<nDigits;iDig++){
e33758d8 113 AliRICHdigit *dig=(AliRICHdigit*)Rich()->Digits(iChamber)->At(iDig);
cbaf35fb 114 Int_t i=dig->X(); Int_t j=dig->Y();
115 if(fHitMap->TestHit(i,j)==kUsed) continue;
116
9d6f9427 117 FormRawCluster(i,j);
cbaf35fb 118
119 if(AliRICHParam::IsResolveClusters()) {
9d6f9427 120 ResolveCluster(); // ResolveCluster serialization will happen inside
cbaf35fb 121 } else {
9d6f9427 122 WriteRawCluster(); // simply output of the RawCluster found without deconvolution
ed3ceb24 123 }
9d6f9427 124 fCurrentCluster.Reset();
cbaf35fb 125 }//digits loop
126
127 delete fHitMap;
ed3ceb24 128 Info("FindClusters","Stop.");
cbaf35fb 129
ed3ceb24 130}//FindClusters()
131//__________________________________________________________________________________________________
132void AliRICHClusterFinder::FindClusterContribs()
133{
134 //finds CombiPid for a given cluster
135 Info("FindClusterContribs","Start");
136
137 TObjArray *pDigits = fCurrentCluster.Digits();
138 Int_t iNmips=0,iNckovs=0,iNfeeds=0;
139 TArrayI contribs(3*fCurrentCluster.Size());
140 Int_t *pindex = new Int_t[3*fCurrentCluster.Size()];
141 for(Int_t iDigN=0;iDigN<fCurrentCluster.Size();iDigN++) {//loop on digits of a given cluster
142 contribs[3*iDigN] =((AliRICHdigit*)pDigits->At(iDigN))->Tid(0);
143 contribs[3*iDigN+1]=((AliRICHdigit*)pDigits->At(iDigN))->Tid(1);
144 contribs[3*iDigN+2]=((AliRICHdigit*)pDigits->At(iDigN))->Tid(2);
145 cout << "TID1 " << contribs[3*iDigN] << " TID2 " << contribs[3*iDigN+1] << " TID3 " << contribs[3*iDigN+2] << endl;
146 }//loop on digits of a given cluster
147 TMath::Sort(contribs.GetSize(),contribs.GetArray(),pindex);
148 for(Int_t iDigN=0;iDigN<3*fCurrentCluster.Size()-1;iDigN++) {//loop on digits to sort Tid
149 cout << " contrib " << contribs[pindex[iDigN]] << " digit " << iDigN <<
150 " contrib " << contribs[pindex[iDigN+1]] << " digit " << iDigN+1 << endl;
151 if(contribs[pindex[iDigN]]!=contribs[pindex[iDigN+1]]) {
152 Int_t code = Rich()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[iDigN]])->GetPdgCode();
153 Double_t charge = Rich()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[iDigN]])->GetPDG()->Charge();
154
155 if(code==50000050) iNckovs++;
156 else if(code==50000051) iNfeeds++;
157 else if(charge!=0) iNmips++;
158 if (contribs[pindex[iDigN+1]]==kBad) break;
159 }
160 }//loop on digits to sort Tid
161 fCurrentCluster.SetCombiPid(iNckovs,iNfeeds,iNmips);
162 fCurrentCluster.Print();
163 delete [] pindex;
164}// FindClusterContribs()
cbaf35fb 165//__________________________________________________________________________________________________
9d6f9427 166void AliRICHClusterFinder::FormRawCluster(Int_t i, Int_t j)
ed3ceb24 167{
168 // Builder of the final Raw Cluster (before deconvolution)
cbaf35fb 169 Info("FormRawCluster","Start with digit(%i,%i)",i,j);
170
9d6f9427 171 fCurrentCluster.AddDigit((AliRICHdigit*) fHitMap->GetHit(i,j));
cbaf35fb 172 fHitMap->FlagHit(i,j);// Flag hit as taken
173
174 Int_t listX[4], listY[4]; // Now look recursively for all neighbours
175 for (Int_t iNeighbour=0;iNeighbour<Rich()->Param()->PadNeighbours(i,j,listX,listY);iNeighbour++)
176 if(fHitMap->TestHit(listX[iNeighbour],listY[iNeighbour])==kUnused)
9d6f9427 177 FormRawCluster(listX[iNeighbour],listY[iNeighbour]);
cbaf35fb 178}//AddDigit2Cluster()
179//__________________________________________________________________________________________________
9d6f9427 180void AliRICHClusterFinder::ResolveCluster()
cbaf35fb 181{// Decluster algorithm
182 Info("ResolveCluster","Start.");
9d6f9427 183
184 fCurrentCluster.CoG(); // first initial approxmation of the CoG...to start minimization.
185 fCurrentCluster.Print();
186 switch (fCurrentCluster.Size()) {
cbaf35fb 187
9d6f9427 188 case 1:
189 WriteRawCluster(); break;
190 case 2:
191 FitCoG();
192 WriteRawCluster(); break;
193
194 default:
195 WriteRawCluster(); break;
196 }
cbaf35fb 197}//ResolveCluster()
198//__________________________________________________________________________________________________
9d6f9427 199void AliRICHClusterFinder::WriteRawCluster()
cbaf35fb 200{// out the current RawCluster
eaf390d9 201 Info("WriteRawCluster","Start.");
cbaf35fb 202
ed3ceb24 203 FindClusterContribs();
9d6f9427 204 Rich()->AddCluster(fCurrentCluster);
cbaf35fb 205
206}//WriteRawCluster()
207//__________________________________________________________________________________________________
9d6f9427 208void AliRICHClusterFinder::FitCoG()
209{// Fit cluster size 2 by single Mathieson
210 Info("FitCoG","Start.");
211
212 TMinuit *pMinuit = new TMinuit(2);
213
214 Double_t arglist;
215 Int_t ierflag = 0;
216
217 static Double_t vstart[2];
218 static Double_t lower[2], upper[2];
219 static Double_t step[2]={0.001,0.001};
220
221 TString chname;
222 Int_t ierflg;
223
224 pMinuit->SetObjectFit((TObject*)this);
225 pMinuit->SetFCN(RICHMinMathieson);
226 pMinuit->mninit(5,10,7);
227
228 vstart[0] = fCurrentCluster.X();
229 vstart[1] = fCurrentCluster.Y();
230
231 lower[0] = vstart[0] - 2*AliRICHParam::PadSizeX();
232 upper[0] = vstart[0] + 2*AliRICHParam::PadSizeX();
233 lower[1] = vstart[1] - 2*AliRICHParam::PadSizeY();
234 upper[1] = vstart[1] + 2*AliRICHParam::PadSizeY();
235
236
237 pMinuit->mnparm(0," x position ",vstart[0],step[0],lower[0],upper[0],ierflag);
238 pMinuit->mnparm(1," y position ",vstart[1],step[1],lower[1],upper[1],ierflag);
239
240 arglist = -1;
241
242 pMinuit->SetPrintLevel(-1);
243 pMinuit->mnexcm("SET NOGR",&arglist, 1, ierflag);
244 pMinuit->mnexcm("SET NOW",&arglist, 1, ierflag);
245 arglist = 1;
246 pMinuit->mnexcm("SET ERR", &arglist, 1,ierflg);
247 arglist = -1;
248 pMinuit->mnexcm("SIMPLEX",&arglist, 0, ierflag);
249 pMinuit->mnexcm("MIGRAD",&arglist, 0, ierflag);
250 pMinuit->mnexcm("EXIT" ,&arglist, 0, ierflag);
251 Double_t xCoG,yCoG;
252 Double_t eps, b1, b2;
253 pMinuit->mnpout(0,chname, xCoG, eps , b1, b2, ierflg);
254 pMinuit->mnpout(1,chname, yCoG, eps , b1, b2, ierflg);
255 delete pMinuit;
256}
257//__________________________________________________________________________________________________
258void RICHMinMathieson(Int_t &, Double_t *, Double_t &chi2, Double_t *par, Int_t iflag)
259{// Minimization function of Mathieson
260
261 AliRICHcluster *pRawCluster = ((AliRICHClusterFinder*)gMinuit->GetObjectFit())->GetCurrentCluster();
262
263 TVector3 centroid(par[0],par[1],0);
264
265 chi2 = 0;
266 Int_t qtot = pRawCluster->Q();
267 for(Int_t i=0;i<pRawCluster->Size();i++) {
268 Int_t padX = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->X();
269 Int_t padY = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->Y();
270 Double_t padQ = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->Q();
271 chi2 += TMath::Power((qtot*AliRICHParam::Loc2PadFrac(centroid,padX,padY)-padQ),2)/padQ;
272 }
273 if(iflag == 3)
274 {
275 cout << " --- end convergence...summary --- " << endl;
276 cout << " x position " << par[0] << endl;
277 cout << " y position " << par[1] << endl;
278 cout << " chi2 " << chi2 << endl;
279 }
280}//RICHMinMathieson()