do not keep digits with null energy, speed up in case of 2 maxima clusters with...

diff --git a/EMCAL/AliEMCALUnfolding.cxx b/EMCAL/AliEMCALUnfolding.cxx
index 8abbc4c585c66b7dc2aca25245f1ad17467c3f21..6fb88cd477d77e63bf451c3cdb68eb1aeda793d3 100644 (file)
--- a/EMCAL/AliEMCALUnfolding.cxx
+++ b/EMCAL/AliEMCALUnfolding.cxx
@@ -155,6 +155,7 @@ void AliEMCALUnfolding::MakeUnfolding()
   // Unfolds clusters using the shape of an ElectroMagnetic shower\r
   // Performs unfolding of all clusters\r
   \r
   // 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
   if(fNumberOfECAClusters > 0){\r
     if (fGeom==0)\r
       AliFatal("Did not get geometry from EMCALLoader") ;\r


@@ -170,7 +171,7 @@ void AliEMCALUnfolding::MakeUnfolding()
         AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;\r
         Float_t * maxAtEnergy = new Float_t[nMultipl] ;\r
         Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;\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
         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


@@ -186,9 +187,9 @@ void AliEMCALUnfolding::MakeUnfolding()
         \r
         delete[] maxAt ;\r
         delete[] maxAtEnergy ;\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
     } // rec point loop\r

-  }\r
+  }//end of check fNumberOfECAClusters\r

   // End of Unfolding of clusters\r
 }\r
 \r
   // End of Unfolding of clusters\r
 }\r
 \r


@@ -203,11 +204,13 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
   //**************************** part 1 *******************************************\r
   // Performs the unfolding of a cluster with nMax overlapping showers \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
   //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
   if (fGeom==0)\r
     AliFatal("Did not get geometry from EMCALLoader") ;\r
   \r


@@ -220,6 +223,22 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
     return kFALSE;\r
   }\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
   //**************************** 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


@@ -227,7 +246,9 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
   // and later correct this number in acordance with actual energy\r
   // deposition\r
   \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
   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
   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


@@ -261,12 +282,18 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
         xpar = fitparameters[iparam] ;\r
         zpar = fitparameters[iparam+1] ;\r
         epar = fitparameters[iparam+2] ;\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
 \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
 \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
     } else AliDebug(1,"Digit NULL part 2!");\r
     \r
   }//digit loop\r


@@ -277,7 +304,7 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
   // to its contribution to efit\r
   \r
   Float_t * energiesList = iniTower->GetEnergiesList() ;\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
   Float_t eDigit = 0. ;\r
   Int_t nSplittedClusters=(Int_t)nPar/3;\r
   \r


@@ -292,6 +319,8 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
   //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
   //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
   iparam = 0 ;\r
   while(iparam < nPar )\r
   {\r


@@ -310,10 +339,12 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
         \r
         EvalParsPhiDependence(digit->GetId(),fGeom);\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
  \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
           continue;\r
         }\r
         \r


@@ -323,6 +354,8 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
         //add energy to temporary matrix\r
         correctedEnergyList[iparam/3*nDigits+iDigit] = eDigit;\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
       } else AliDebug(1,"NULL digit part 3");\r
     }//digit loop \r
     iparam += 3 ;\r


@@ -330,7 +363,9 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
   \r
   //**************************** sub-part 3.2 *************************************\r
   //here we correct energy for each cell and cluster\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
   Int_t maximumIndex=0;\r
   Bool_t isAnyBelowThreshold=kFALSE;\r
   //  Float_t Threshold=0.01;\r


@@ -339,7 +374,7 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
   for(iDigit = 0 ; iDigit < nDigits ; iDigit++)\r
   {\r
     isAnyBelowThreshold=kFALSE;\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
     for(iparam = 0 ; iparam < nPar ; iparam+=3)\r
     {\r
 \r


@@ -352,9 +387,12 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
     }//end of loop over clusters after unfolding\r
     \r
     if(!isAnyBelowThreshold) continue; //no cluster-cell below threshold \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
     //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
 \r
     if(maximumEne < fThreshold) \r
     {//add all cluster cells and put energy into max index, other set to 0\r


@@ -369,17 +407,17 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
     }//end if\r
     \r
     //divide energy of cell below threshold in the correct proportion and add to other cells\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
     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
         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){//calculate fraction\r
         energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3*nDigits+iDigit] / maximumEne;\r


@@ -387,7 +425,7 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
       \r
       for(iparam = 0 ; iparam < nPar ; iparam+=3)\r
       {//add energy from cells below threshold to others\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
         else\r
         {\r
           for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3)\r


@@ -395,7 +433,7 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
             correctedEnergyList[iparam2/3*nDigits+iDigit] += (energyFraction[iparam2/3] * \r
                                                       correctedEnergyList[iparam/3*nDigits+iDigit]) ;\r
           }//inner loop\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
         }\r
       }\r
     }\r


@@ -404,15 +442,24 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
       //digit energy to be set to 0\r
       for(iparam = 0 ; iparam < nPar ; iparam+=3)\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
       }\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
   //**************************** 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
   iparam = 0 ;\r
   while(iparam < nPar )\r
   {\r


@@ -425,39 +472,53 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
     (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;\r
     recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;\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
       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
           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
           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
         }\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

-    }\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
+\r

     iparam += 3 ;\r
   }//while\r
   \r
     iparam += 3 ;\r
   }//while\r
   \r


@@ -465,6 +526,7 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
   delete[] efit ;\r
   delete[] correctedEnergyList ;\r
   \r
   delete[] efit ;\r
   delete[] correctedEnergyList ;\r
   \r

+  //  cout<<"end of unfolding check part 3.3"<<endl;\r

   return kTRUE;\r
 }\r
 \r
   return kTRUE;\r
 }\r
 \r