// Unfolds clusters using the shape of an ElectroMagnetic shower\r
// Performs unfolding of all clusters\r
\r
+ //cout<<"fNumberOfECAClusters "<<fNumberOfECAClusters<<endl;\r
if(fNumberOfECAClusters > 0){\r
if (fGeom==0)\r
AliFatal("Did not get geometry from EMCALLoader") ;\r
AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;\r
Float_t * maxAtEnergy = new Float_t[nMultipl] ;\r
Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;\r
- \r
+ //cout<<"nMax "<<nMax<<endl;\r
if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0\r
if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){\r
fRecPoints->Remove(recPoint);\r
\r
delete[] maxAt ;\r
delete[] maxAtEnergy ;\r
- } else AliError("RecPoint NULL");\r
+ } else AliError("RecPoint NULL"); //end of check if recPoint exist\r
} // rec point loop\r
- }\r
+ }//end of check fNumberOfECAClusters\r
// End of Unfolding of clusters\r
}\r
\r
//**************************** part 1 *******************************************\r
// Performs the unfolding of a cluster with nMax overlapping showers \r
\r
+ //cout<<"unfolding check here part 1"<<endl;\r
//printf("Original cluster E %f\n",iniTower->GetEnergy());\r
\r
Int_t nPar = 3 * nMax ;\r
Float_t * fitparameters = new Float_t[nPar] ;\r
\r
+ //cout<<"number of parameters "<<nPar<<" nMax "<<nMax<<endl;\r
if (fGeom==0)\r
AliFatal("Did not get geometry from EMCALLoader") ;\r
\r
return kFALSE;\r
}\r
\r
+// for(Int_t iii=0;iii<nPar;iii++){\r
+// cout<<" param "<<iii<<" from fit "<<fitparameters[iii];\r
+// if((iii+1)%3==0)cout<<endl;\r
+// }\r
+\r
+//speed up solution for clusters with 2 maxima where one maximum is below threshold fThreshold\r
+ if(nMax==2){\r
+ if(fitparameters[2]<fThreshold || fitparameters[5]<fThreshold){\r
+ //cout<<"one of fitted energy below threshold"<<endl;\r
+ AliDebug(1,"One of fitted energy below threshold");\r
+ iniTower->SetNExMax(1) ;\r
+ delete[] fitparameters ;\r
+ return kFALSE;\r
+ }\r
+ }\r
+\r
//**************************** part 2 *******************************************\r
// create unfolded rec points and fill them with new energy lists\r
// First calculate energy deposited in each sell in accordance with\r
// and later correct this number in acordance with actual energy\r
// deposition\r
\r
+ // cout<<"unfolding check here part 2"<<endl;\r
Int_t nDigits = iniTower->GetMultiplicity() ;\r
+ // cout<<"cluster multiplicity "<<nDigits<<" "<< iniTower->GetMultiplicity() <<endl;\r
Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells\r
Float_t xpar=0.,zpar=0.,epar=0. ;//center of gravity in cell units\r
\r
xpar = fitparameters[iparam] ;\r
zpar = fitparameters[iparam+1] ;\r
epar = fitparameters[iparam+2] ;\r
+// cout<<" xpar from fit "<<xpar<<" "<< fitparameters[iparam]<<endl; \r
+// cout<<" zpar from fit "<<zpar<<" "<< fitparameters[iparam+1]<<endl; \r
+// cout<<" epar from fit "<<epar<<" "<< fitparameters[iparam+2]<<endl; \r
+\r
\r
efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
+ // cout<<"idigit "<<iDigit<<" efit[idigit]="<<efit[iDigit]<<" iparam "<<iparam<<" "<< (Float_t)iphi - xpar <<" "<< (Float_t)ieta - zpar<<endl;\r
\r
iparam += 3 ;\r
}\r
\r
+ // cout<<"idigit "<<iDigit<<" total efit[idigit]="<<efit[iDigit]<<endl;\r
} else AliDebug(1,"Digit NULL part 2!");\r
\r
}//digit loop\r
// to its contribution to efit\r
\r
Float_t * energiesList = iniTower->GetEnergiesList() ;\r
- Float_t ratio = 0 ;\r
+ Float_t ratio = 0. ;//0 -> 0. changed\r
Float_t eDigit = 0. ;\r
Int_t nSplittedClusters=(Int_t)nPar/3;\r
\r
//If not the energy from a given cell in the cluster is divided in correct proportions \r
//in accordance to the other clusters and added to them and set to 0.\r
\r
+ // cout<<"unfolding check here part 3.1"<<endl;\r
+\r
iparam = 0 ;\r
while(iparam < nPar )\r
{\r
\r
EvalParsPhiDependence(digit->GetId(),fGeom);\r
\r
+ // cout<<"iparam "<<iparam<<" iDigit "<<iDigit<<" efit[iDigit] "<<efit[iDigit]<<endl;\r
\r
if(efit[iDigit]==0) \r
{//just for sure\r
- correctedEnergyList[iparam/3*nDigits+iDigit] = 0;//correction here\r
+ // cout<<"energy =0"<<endl;\r
+ correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;//correction here\r
continue;\r
}\r
\r
//add energy to temporary matrix\r
correctedEnergyList[iparam/3*nDigits+iDigit] = eDigit;\r
\r
+ // cout<<"iDigit "<<iDigit<<" ratio "<<ratio<<" efit[iDigit] "<<efit[iDigit]<<" iparam "<<iparam<< "iparam/3*nDigits+iDigit "<<iparam/3*nDigits+iDigit<<" eDigit "<<eDigit<<" correctedEnergyList[] "<< correctedEnergyList[iparam/3*nDigits+iDigit]<<endl;\r
+\r
} else AliDebug(1,"NULL digit part 3");\r
}//digit loop \r
iparam += 3 ;\r
\r
//**************************** sub-part 3.2 *************************************\r
//here we correct energy for each cell and cluster\r
- Float_t maximumEne=0;\r
+ // cout<<"unfolding check here part 3.2"<<endl;\r
+\r
+ Float_t maximumEne=0.;//0 -> 0. changed\r
Int_t maximumIndex=0;\r
Bool_t isAnyBelowThreshold=kFALSE;\r
// Float_t Threshold=0.01;\r
for(iDigit = 0 ; iDigit < nDigits ; iDigit++)\r
{\r
isAnyBelowThreshold=kFALSE;\r
- maximumEne=0;\r
+ maximumEne=0.;//0 -> 0. changed\r
for(iparam = 0 ; iparam < nPar ; iparam+=3)\r
{\r
\r
}//end of loop over clusters after unfolding\r
\r
if(!isAnyBelowThreshold) continue; //no cluster-cell below threshold \r
- \r
- // printf("Correct E cell %f < %f, org Digit index %d, e = %f\n",correctedEnergyList[iparam/3*nDigits+iDigit],fThreshold,iDigit, energiesList[iDigit]);\r
+\r
+\r
+ //printf("Correct E cell %f < %f, org Digit index %d, e = %f\n",correctedEnergyList[iparam/3*nDigits+iDigit],fThreshold,iDigit, energiesList[iDigit]);\r
//if( energiesList[iDigit] < correctedEnergyList[iparam/3*nDigits+iDigit]) printf("\t What? \n");\r
+ //cout<<" Correct E cell "<<correctedEnergyList[iparam/3*nDigits+iDigit]<<" < "<<fThreshold<<", org Digit index "<<iDigit<<", e = "<<energiesList[iDigit]<<endl;\r
+\r
\r
if(maximumEne < fThreshold) \r
{//add all cluster cells and put energy into max index, other set to 0\r
}//end if\r
\r
//divide energy of cell below threshold in the correct proportion and add to other cells\r
- maximumEne=0;//not used any more so use it for the energy sum\r
+ maximumEne=0.;//not used any more so use it for the energy sum //0 -> 0. changed\r
for(iparam = 0 ; iparam < nPar ; iparam+=3)\r
{//calculate energy sum\r
if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold) energyFraction[iparam/3]=0;\r
else \r
{\r
- energyFraction[iparam/3]=1;\r
+ energyFraction[iparam/3]=1.;//1 -> 1. changed\r
maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];\r
}\r
}//end of loop over clusters after unfolding\r
- if(maximumEne>0)\r
+ if(maximumEne>0.)//0 -> 0. changed\r
{\r
for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction\r
energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3*nDigits+iDigit] / maximumEne;\r
\r
for(iparam = 0 ; iparam < nPar ; iparam+=3)\r
{//add energy from cells below threshold to others\r
- if(energyFraction[iparam/3]>0) continue;\r
+ if(energyFraction[iparam/3]>0.) continue;//0 -> 0. changed\r
else\r
{\r
for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3)\r
correctedEnergyList[iparam2/3*nDigits+iDigit] += (energyFraction[iparam2/3] * \r
correctedEnergyList[iparam/3*nDigits+iDigit]) ;\r
}//inner loop\r
- correctedEnergyList[iparam/3*nDigits+iDigit] = 0;\r
+ correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;//0 -> 0. changed\r
}\r
}\r
}\r
//digit energy to be set to 0\r
for(iparam = 0 ; iparam < nPar ; iparam+=3)\r
{\r
- correctedEnergyList[iparam/3*nDigits+iDigit] = 0;\r
+ correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;//0 -> 0. changed\r
}\r
}//new adam correction for is energy>0\r
\r
}//end of loop over digits\r
delete[] energyFraction;\r
\r
+ //cout section\r
+// cout<<"nDigits "<<nDigits<<endl;\r
+// for(Int_t iii=0;iii<nDigits*nPar/3;++iii){\r
+// cout<<"correctedEnergyList["<<iii<<"]="<<correctedEnergyList[iii]<<endl;\r
+// }\r
+ //end of cout section\r
+\r
//**************************** sub-part 3.3 *************************************\r
//here we add digits to recpoints with corrected energy\r
+ // cout<<"unfolding check here part 3.3"<<endl;\r
+\r
iparam = 0 ;\r
while(iparam < nPar )\r
{\r
(*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;\r
recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;\r
\r
- if(recPoint)\r
- {\r
- fNumberOfECAClusters++ ;\r
+ if(recPoint){//recPoint preseent -> good\r
recPoint->SetNExMax(nSplittedClusters) ;\r
\r
for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)\r
{\r
digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r
\r
- \r
- if(digit && correctedEnergyList[iparam/3*nDigits+iDigit]>0. )\r
- {\r
+ if(digit && correctedEnergyList[iparam/3*nDigits+iDigit]>0. ){\r
+ // cout<<"idigit,nMax"<<iDigit<<","<<iparam/3<<" correctedEnergyList["<<iparam/3*nDigits+iDigit<<"]="<<correctedEnergyList[iparam/3*nDigits+iDigit]<<endl;\r
if(correctedEnergyList[iparam/3*nDigits+iDigit]<fThreshold) printf("Final E cell %f < %f\n",correctedEnergyList[iparam/3*nDigits+iDigit],fThreshold);\r
recPoint->AddDigit( *digit, correctedEnergyList[iparam/3*nDigits+iDigit], kFALSE ) ; //FIXME, need to study the shared case\r
- } else \r
- {\r
+ } else {\r
AliDebug(1,Form("NULL digit part3.3 or NULL energy=%f",correctedEnergyList[iparam/3*nDigits+iDigit]));\r
+ //cout<<"nDigits "<<nDigits<<" iParam/3 "<<iparam/3<< endl;\r
}\r
- }//digit loop \r
- } else AliError("NULL RecPoint");\r
- \r
- //protection from recpoint with no digits\r
- if(recPoint->GetMultiplicity()==0)\r
- {\r
- delete (*fRecPoints)[fNumberOfECAClusters];\r
- //cout<<"size fRecPoints before "<<fRecPoints->GetSize()<<endl;\r
- fRecPoints->RemoveAt(fNumberOfECAClusters);\r
- //cout<<"size fRecPoints after "<<fRecPoints->GetSize()<<endl;\r
- fNumberOfECAClusters--;\r
- nSplittedClusters--;\r
+\r
+ }//digit loop\r
+\r
+ if(recPoint->GetMultiplicity()==0){//recpoint exists but no digits associated -> remove from list\r
+ delete (*fRecPoints)[fNumberOfECAClusters];\r
+ //cout<<"size fRecPoints before "<<fRecPoints->GetSize()<<endl;\r
+ fRecPoints->RemoveAt(fNumberOfECAClusters);\r
+ //cout<<"size fRecPoints after "<<fRecPoints->GetSize()<<endl;\r
+ fNumberOfECAClusters--;\r
+ nSplittedClusters--;\r
+ } else {//recPoint exists and has digits associated -> very good increase number of clusters \r
+ fNumberOfECAClusters++ ; \r
+ }\r
\r
- }\r
+ } else {//recPoint empty -> remove from list\r
+ AliError("NULL RecPoint");\r
+ \r
+ //protection from recpoint with no digits\r
+ // cout<<"multi rec "<<recPoint->GetMultiplicity()<<endl;\r
+ if(recPoint->GetMultiplicity()==0)\r
+ {\r
+ delete (*fRecPoints)[fNumberOfECAClusters];\r
+ //cout<<"size fRecPoints before "<<fRecPoints->GetSize()<<endl;\r
+ fRecPoints->RemoveAt(fNumberOfECAClusters);\r
+ //cout<<"size fRecPoints after "<<fRecPoints->GetSize()<<endl;\r
+ fNumberOfECAClusters--;\r
+ nSplittedClusters--;\r
+ \r
+ }\r
\r
+ }\r
+\r
iparam += 3 ;\r
}//while\r
\r
delete[] efit ;\r
delete[] correctedEnergyList ;\r
\r
+ // cout<<"end of unfolding check part 3.3"<<endl;\r
return kTRUE;\r
}\r
\r