ClassImp(AliRICHClusterFinder)
//__________________________________________________________________________________________________
-AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH)
+AliRICHClusterFinder::AliRICHClusterFinder(AliRunLoader *pRunLoader)
{//main ctor
- fRICH = pRICH;
+ fRICH = (AliRICH*) pRunLoader->GetAliRun()->GetDetector("RICH");
+
AliDebug(1,"main ctor Start.");
fDigitMap = 0;
AliDebug(1,"main ctor Stop.");
}//main ctor
//__________________________________________________________________________________________________
-void AliRICHClusterFinder::Exec()
+void AliRICHClusterFinder::Exec(const Option_t *)
{
//Main method of cluster finder. Loops on events and chambers, everything else is done in FindClusters()
AliDebug(1,"Exec Start.");
R()->GetLoader() ->LoadDigits();
// R()->GetLoader()->GetRunLoader()->LoadHeader();
- R()->GetLoader()->GetRunLoader()->LoadKinematics(); //header is already loaded
+ if(!R()->GetLoader()->GetRunLoader()->TreeK()) R()->GetLoader()->GetRunLoader()->LoadKinematics();
for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop
AliDebug(1,Form("Processing event %i...",iEventN));
R()->GetLoader()->GetRunLoader()->GetEvent(iEventN);
R()->GetLoader()->MakeTree("R"); R()->MakeBranch("R");
- R()->ResetDigits(); R()->ResetClusters();
+ R()->DigitsReset(); R()->ClustersReset();
R()->GetLoader()->TreeD()->GetEntry(0);
for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++){//chambers loop
R()->GetLoader()->TreeR()->Fill(); R()->GetLoader()->WriteRecPoints("OVERWRITE");//write out clusters for current event
}//events loop
- R()->ResetDigits();//reset and unload everything
- R()->ResetClusters();
+ R()->DigitsReset();//reset and unload everything
+ R()->ClustersReset();
R()->GetLoader() ->UnloadDigits();
R()->GetLoader() ->UnloadRecPoints();
// R()->GetLoader()->GetRunLoader()->UnloadHeader();
fDigitMap=new AliRICHMap(R()->Digits(iChamber));//create digit map for the given chamber
for(Int_t iDigN=0;iDigN<iNdigits;iDigN++){//digits loop for a given chamber
- AliRICHdigit *dig=(AliRICHdigit*)R()->Digits(iChamber)->At(iDigN);
+ AliRICHDigit *dig=(AliRICHDigit*)R()->Digits(iChamber)->At(iDigN);
Int_t i=dig->X(); Int_t j=dig->Y();
if(fDigitMap->TestHit(i,j)==kUsed) continue;//this digit is already taken, go after next digit
AliDebug(1,"Stop.");
}//FindClusters()
//__________________________________________________________________________________________________
-void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster)
+void AliRICHClusterFinder::FindClusterContribs(AliRICHCluster *pCluster)
{
//Finds cerenkov-feedback-mip mixture for a given cluster
AliDebug(1,"Start.");ToAliDebug(1,pCluster->Print());
TArrayI contribs(3*pCluster->Size());
Int_t *pindex = new Int_t[3*pCluster->Size()];
for(Int_t iDigN=0;iDigN<pCluster->Size();iDigN++) {//loop on digits of a given cluster
- contribs[3*iDigN] =((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(0);
+ contribs[3*iDigN] =((AliRICHDigit*)pDigits->At(iDigN))->GetTrack(0);
if (contribs[3*iDigN] >= 10000000) contribs[3*iDigN] = 0;
- contribs[3*iDigN+1]=((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(1);
+ contribs[3*iDigN+1]=((AliRICHDigit*)pDigits->At(iDigN))->GetTrack(1);
if (contribs[3*iDigN+1] >= 10000000) contribs[3*iDigN+1] = 0;
- contribs[3*iDigN+2]=((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(2);
+ contribs[3*iDigN+2]=((AliRICHDigit*)pDigits->At(iDigN))->GetTrack(2);
if (contribs[3*iDigN+2] >= 10000000) contribs[3*iDigN+2] = 0;
}//loop on digits of a given cluster
TMath::Sort(contribs.GetSize(),contribs.GetArray(),pindex);
for(Int_t iDigN=0;iDigN<3*pCluster->Size()-1;iDigN++) {//loop on digits to sort tids
AliDebug(1,Form("%4i for digit n. %4i",contribs[pindex[iDigN]],iDigN));
if(contribs[pindex[iDigN]]!=contribs[pindex[iDigN+1]]) {
- TParticle* particle = pStack->Particle(contribs[pindex[iDigN]]);
+ Int_t thecontrib = contribs[pindex[iDigN]];
+ if (thecontrib>=pStack->GetNtrack()) continue;//PH this should not happen
+ TParticle* particle = pStack->Particle(thecontrib);
+ if (!particle) continue;//PH this should not happen
Int_t code = particle->GetPdgCode();
Double_t charge = 0;
if(particle->GetPDG()) charge=particle->GetPDG()->Charge();
}
}//loop on digits to sort Tid
- if (contribs[pindex[3*pCluster->Size()-1]]!=kBad) {
+ if (contribs[pindex[3*pCluster->Size()-1]]!=-1) {
+ Int_t thecontrib = contribs[pindex[3*pCluster->Size()-1]];
+ if (thecontrib<pStack->GetNtrack()){
+ //PH the opposite should not happen
- TParticle* particle = pStack->Particle(contribs[pindex[3*pCluster->Size()-1]]);
- Int_t code = particle->GetPdgCode();
- Double_t charge = 0;
- if(particle->GetPDG()) charge=particle->GetPDG()->Charge();
- AliDebug(1,Form(" charge of particle %f",charge));
- if(code==50000050) iNckovs++;
- if(code==50000051) iNfeeds++;
- if(charge!=0) iNmips++;
+ TParticle* particle = pStack->Particle(thecontrib);
+ if (particle) {
+ //PH the opposite should not happen
+ Int_t code = particle->GetPdgCode();
+ Double_t charge = 0;
+ if(particle->GetPDG()) charge=particle->GetPDG()->Charge();
+ AliDebug(1,Form(" charge of particle %f",charge));
+ if(code==50000050) iNckovs++;
+ if(code==50000051) iNfeeds++;
+ if(charge!=0) iNmips++;
+ }
+ }
}
pCluster->CFM(iNckovs,iNfeeds,iNmips);
void AliRICHClusterFinder::FormRawCluster(Int_t i, Int_t j)
{
//Builds the raw cluster (before deconvolution). Starts from the first pad (i,j) then calls itself recursevly for all neighbours.
- AliDebug(1,Form("Start with digit(%i,%i) Q=%f",i,j,((AliRICHdigit*)fDigitMap->GetHit(i,j))->Q()));
+ AliDebug(1,Form("Start with digit(%i,%i) Q=%f",i,j,((AliRICHDigit*)fDigitMap->GetHit(i,j))->Q()));
- fRawCluster.AddDigit((AliRICHdigit*) fDigitMap->GetHit(i,j));//take this pad in cluster
+ fRawCluster.AddDigit((AliRICHDigit*) fDigitMap->GetHit(i,j));//take this pad in cluster
fDigitMap->FlagHit(i,j);//flag this pad as taken
Int_t listX[4], listY[4]; // Now look recursively for all neighbours
AliDebug(1,Form("Cluster size of the Raw cluster ---> %i",fRawCluster.Size()));
for(Int_t iDig1=0;iDig1<fRawCluster.Size();iDig1++) {
Int_t iNotMax = 0;
- AliRICHdigit *pDig1 = (AliRICHdigit *)fRawCluster.Digits()->At(iDig1);
+ AliRICHDigit *pDig1 = (AliRICHDigit *)fRawCluster.Digits()->At(iDig1);
if(!pDig1) {fNlocals=0;return;}
TVector pad1 = pDig1->Pad();
Int_t padQ1 = (Int_t)(pDig1->Q()+0.1);
Int_t padC1 = pDig1->ChFbMi();
for(Int_t iDig2=0;iDig2<fRawCluster.Size();iDig2++) {
- AliRICHdigit *pDig2 = (AliRICHdigit *)fRawCluster.Digits()->At(iDig2);
+ AliRICHDigit *pDig2 = (AliRICHDigit *)fRawCluster.Digits()->At(iDig2);
if(!pDig2) {fNlocals=0;return;}
TVector pad2 = pDig2->Pad();
Int_t padQ2 = (Int_t)(pDig2->Q()+0.1);
}
}
if(iNotMax==0) {
- TVector2 x2=AliRICHParam::Pad2Loc(pad1);
- fLocalX[fNlocals]=x2.X();fLocalY[fNlocals]=x2.Y();
- fLocalQ[fNlocals] = (Double_t)padQ1;
- fLocalC[fNlocals] = padC1;
- fNlocals++;
+ if (fNlocals<100) {
+ TVector2 x2=AliRICHParam::Pad2Loc(pad1);
+ fLocalX[fNlocals]=x2.X();fLocalY[fNlocals]=x2.Y();
+ fLocalQ[fNlocals] = (Double_t)padQ1;
+ fLocalC[fNlocals] = padC1;
+ fNlocals++;
+ }
}
}
AliDebug(1,Form("Number of local maxima found ---> %i",fNlocals));
AliDebug(1,"Start.");
FindClusterContribs(&fRawCluster);
- R()->AddCluster(fRawCluster);
+ R()->ClusterAdd(fRawCluster);
ToAliDebug(1,fRawCluster.Print()); AliDebug(1,"Stop.");
}//WriteRawCluster()
AliDebug(1,"Start.");
FindClusterContribs(&fResolvedCluster);
- R()->AddCluster(fResolvedCluster);
+ R()->ClusterAdd(fResolvedCluster);
ToAliDebug(1,fResolvedCluster.Print()); AliDebug(1,"Stop.");
}//WriteResolvedCluster()
{
//Mathieson minimization function
- AliRICHcluster *pRawCluster = ((AliRICHClusterFinder*)gMinuit->GetObjectFit())->GetRawCluster();
+ AliRICHCluster *pRawCluster = ((AliRICHClusterFinder*)gMinuit->GetObjectFit())->GetRawCluster();
TVector2 centroid[50];
Double_t q[50];
chi2 = 0;
Int_t qtot = pRawCluster->Q();
for(Int_t i=0;i<pRawCluster->Size();i++) {
- TVector pad=((AliRICHdigit *)pRawCluster->Digits()->At(i))->Pad();
- Double_t padQ = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->Q();
+ TVector pad=((AliRICHDigit *)pRawCluster->Digits()->At(i))->Pad();
+ Double_t padQ = ((AliRICHDigit *)pRawCluster->Digits()->At(i))->Q();
Double_t qfracpar=0;
for(Int_t j=0;j<nFunctions;j++) {
qfracpar += q[j]*AliRICHParam::FracQdc(centroid[j],pad);