| b410dc6a |
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 | |
| 16 | //_________________________________________________________________________ |
| 17 | // Base class for the cluster unfolding algorithm |
| 18 | //*-- Author: Adam Matyja (SUBATECH) |
| 19 | // Based on unfolding in clusterizerv1 done by Cynthia Hadjidakis |
| 20 | //-- Unfolding for eta~0: Cynthia Hadjidakis - still in AliEMCALCLusterizerv1 |
| 21 | //-- Unfolding extension for whole EMCAL: Adam Matyja (SUBATECH & INP PAN) |
| 22 | // |
| 23 | // unfolds the clusters having several local maxima. |
| 24 | ////////////////////////////////////////////////////////////////////////////// |
| 25 | |
| 26 | // --- ROOT system --- |
| 27 | #include "TClonesArray.h" |
| 28 | #include <TMath.h> |
| 29 | #include <TMinuit.h> |
| 30 | |
| 31 | // --- Standard library --- |
| 32 | #include <cassert> |
| 33 | |
| 34 | // --- AliRoot header files --- |
| 35 | #include "AliEMCALUnfolding.h" |
| 36 | #include "AliEMCALGeometry.h" |
| 37 | #include "AliRunLoader.h" |
| 38 | #include "AliRun.h" |
| 39 | #include "AliEMCAL.h" |
| 40 | #include "AliEMCALRecParam.h" |
| 41 | #include "AliEMCALRecPoint.h" |
| 42 | #include "AliEMCALDigit.h" |
| 43 | #include "AliEMCALReconstructor.h" |
| 44 | |
| 45 | #include "AliLog.h" |
| 46 | #include "AliCDBManager.h" |
| 47 | class AliCDBStorage; |
| 48 | #include "AliCDBEntry.h" |
| 49 | |
| 50 | Double_t AliEMCALUnfolding::fgSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332}; |
| 51 | Double_t AliEMCALUnfolding::fgPar5[3]={12.31,-0.007381,-0.06936}; |
| 52 | Double_t AliEMCALUnfolding::fgPar6[3]={0.05452,0.0001228,0.001361}; |
| 53 | |
| 54 | ClassImp(AliEMCALUnfolding) |
| 55 | |
| 56 | //____________________________________________________________________________ |
| 57 | AliEMCALUnfolding::AliEMCALUnfolding(): |
| 58 | fNumberOfECAClusters(0), |
| 59 | fECALocMaxCut(0), |
| 60 | fThreshold(0.01),//10 MeV |
| 61 | fRejectBelowThreshold(0),//split |
| 62 | fGeom(NULL), |
| 63 | fRecPoints(NULL), |
| 64 | fDigitsArr(NULL) |
| 65 | { |
| 66 | // ctor with the indication of the file where header Tree and digits Tree are stored |
| 67 | Init() ; |
| 68 | } |
| 69 | |
| 70 | //____________________________________________________________________________ |
| 71 | AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry): |
| 72 | fNumberOfECAClusters(0), |
| 73 | fECALocMaxCut(0), |
| 74 | fThreshold(0.01),//10 MeV |
| 75 | fRejectBelowThreshold(0),//split |
| 76 | fGeom(geometry), |
| 77 | fRecPoints(NULL), |
| 78 | fDigitsArr(NULL) |
| 79 | { |
| 80 | // ctor with the indication of the file where header Tree and digits Tree are stored |
| 81 | // use this contructor to avoid usage of Init() which uses runloader |
| 82 | // change needed by HLT - MP |
| 83 | if (!fGeom) |
| 84 | { |
| 85 | AliFatal("AliEMCALUnfolding: Geometry not initialized."); |
| 86 | } |
| 87 | |
| 88 | } |
| 89 | |
| 90 | //____________________________________________________________________________ |
| 91 | AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry,Float_t ECALocMaxCut,Double_t *SSPars,Double_t *Par5,Double_t *Par6): |
| 92 | fNumberOfECAClusters(0), |
| 93 | fECALocMaxCut(ECALocMaxCut), |
| 94 | fThreshold(0.01),//10 MeV |
| 95 | fRejectBelowThreshold(0),//split |
| 96 | fGeom(geometry), |
| 97 | fRecPoints(NULL), |
| 98 | fDigitsArr(NULL) |
| 99 | { |
| 100 | // ctor with the indication of the file where header Tree and digits Tree are stored |
| 101 | // use this contructor to avoid usage of Init() which uses runloader |
| 102 | // change needed by HLT - MP |
| 103 | if (!fGeom) |
| 104 | { |
| 105 | AliFatal("AliEMCALUnfolding: Geometry not initialized."); |
| 106 | } |
| 107 | Int_t i=0; |
| 108 | for (i = 0; i < 8; i++) fgSSPars[i] = SSPars[i]; |
| 109 | for (i = 0; i < 3; i++) { |
| 110 | fgPar5[i] = Par5[i]; |
| 111 | fgPar6[i] = Par6[i]; |
| 112 | } |
| 113 | |
| 114 | } |
| 115 | |
| 116 | //____________________________________________________________________________ |
| 117 | void AliEMCALUnfolding::Init() |
| 118 | { |
| 119 | // Make all memory allocations which can not be done in default constructor. |
| 120 | // Attach the Clusterizer task to the list of EMCAL tasks |
| 121 | |
| 122 | AliRunLoader *rl = AliRunLoader::Instance(); |
| 123 | if (rl && rl->GetAliRun()){ |
| 124 | AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL")); |
| 125 | if(emcal)fGeom = emcal->GetGeometry(); |
| 126 | } |
| 127 | |
| 128 | if(!fGeom) |
| 129 | fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName()); |
| 130 | |
| 131 | AliDebug(1,Form("geom %p",fGeom)); |
| 132 | |
| 133 | if(!gMinuit) |
| 134 | // gMinuit = new TMinuit(100) ;//the same is in FindFitV2 |
| 135 | gMinuit = new TMinuit(30) ;//the same is in FindFitV2 |
| 136 | |
| 137 | } |
| 138 | |
| 139 | //____________________________________________________________________________ |
| 140 | AliEMCALUnfolding::~AliEMCALUnfolding() |
| 141 | { |
| 142 | // dtor |
| 143 | } |
| 144 | |
| 145 | //____________________________________________________________________________ |
| 146 | void AliEMCALUnfolding::SetInput(Int_t numberOfECAClusters,TObjArray *recPoints,TClonesArray *digitsArr) |
| 147 | { |
| 148 | // |
| 149 | //Set input for unfolding purposes |
| 150 | // |
| 151 | SetNumberOfECAClusters(numberOfECAClusters); |
| 152 | SetRecPoints(recPoints); |
| 153 | SetDigitsArr(digitsArr); |
| 154 | } |
| 155 | |
| 156 | //____________________________________________________________________________ |
| 157 | void AliEMCALUnfolding::MakeUnfolding() |
| 158 | { |
| 159 | // Unfolds clusters using the shape of an ElectroMagnetic shower |
| 160 | // Performs unfolding of all clusters |
| 161 | |
| 162 | AliDebug(4,Form(" V1: total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); |
| 163 | if(fNumberOfECAClusters > 0){ |
| 164 | if (fGeom==0) |
| 165 | AliFatal("Did not get geometry from EMCALLoader") ; |
| 166 | //Int_t nModulesToUnfold = fGeom->GetNCells(); |
| 167 | |
| 168 | Int_t numberOfClustersToUnfold=fNumberOfECAClusters; |
| 169 | //we unfold only clusters present in the array untill now |
| 170 | //fNumberOfECAClusters may change due to unfilded clusters |
| 171 | //so 0 to numberOfClustersToUnfold-1: clusters before unfolding |
| 172 | //numberOfClustersToUnfold to the end: new clusters from unfolding |
| 173 | //of course numberOfClustersToUnfold also is decreased but we don't loop over clusters added in UF |
| 174 | Int_t index ; |
| 175 | for(index = 0 ; index < numberOfClustersToUnfold ; index++){ |
| 176 | AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ; |
| 177 | if(recPoint){ |
| 178 | Int_t nMultipl = recPoint->GetMultiplicity() ; |
| 179 | AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ; |
| 180 | Float_t * maxAtEnergy = new Float_t[nMultipl] ; |
| 181 | Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ; |
| 182 | if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0 |
| 183 | AliDebug(4,Form(" *** V1+UNFOLD *** Cluster index before UF %d",fNumberOfECAClusters)); |
| 184 | if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){ |
| 185 | //if unfolding correct remove old recPoint |
| 186 | fRecPoints->Remove(recPoint); |
| 187 | fRecPoints->Compress() ;//is it really needed |
| 188 | index-- ; |
| 189 | fNumberOfECAClusters-- ; |
| 190 | numberOfClustersToUnfold--; |
| 191 | } |
| 192 | AliDebug(4,Form(" Cluster index after UF %d",fNumberOfECAClusters)); |
| 193 | } else{ |
| 194 | recPoint->SetNExMax(1) ; //Only one local maximum |
| 195 | } |
| 196 | |
| 197 | delete[] maxAt ; |
| 198 | delete[] maxAtEnergy ; |
| 199 | } else { |
| 200 | //AliError("RecPoint NULL"); //end of check if recPoint exist |
| 201 | Error("MakeUnfolding", "RecPoint NULL, index = %d, fNumberOfECAClusters = %d, numberOfClustersToUnfold = %d",index,fNumberOfECAClusters,numberOfClustersToUnfold) ; |
| 202 | } |
| 203 | } // rec point loop |
| 204 | }//end of check fNumberOfECAClusters |
| 205 | // End of Unfolding of clusters |
| 206 | |
| 207 | AliDebug(4,Form(" V1+UNFOLD: total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); |
| 208 | // for(Int_t i=0;i<fNumberOfECAClusters;i++){ |
| 209 | // AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(i)); |
| 210 | // Int_t nMultipl = recPoint->GetMultiplicity() ; |
| 211 | // Double_t energy=recPoint->GetEnergy(); |
| 212 | // Int_t absIdMaxDigit=recPoint->GetAbsIdMaxDigit(); |
| 213 | // Int_t sm=recPoint->GetSuperModuleNumber(); |
| 214 | // Double_t pointEne=recPoint->GetPointEnergy(); |
| 215 | // Float_t maxEne=recPoint->GetMaximalEnergy(); |
| 216 | // Int_t maxEneInd=recPoint->GetMaximalEnergyIndex(); |
| 217 | // printf(" cluster %d,ncells %d,ene %f,absIdMaxCell %d,sm %d,pointEne %f,maxEne %f,maxEneInd %d\n",i,nMultipl,energy,absIdMaxDigit,sm,pointEne,maxEne,maxEneInd); |
| 218 | // } |
| 219 | |
| 220 | } |
| 221 | |
| 222 | //____________________________________________________________________________ |
| 223 | Int_t AliEMCALUnfolding::UnfoldOneCluster(AliEMCALRecPoint * iniTower, |
| 224 | Int_t nMax, |
| 225 | AliEMCALDigit ** maxAt, |
| 226 | Float_t * maxAtEnergy, |
| 227 | TObjArray *list) |
| 228 | { |
| 229 | // Input one cluster |
| 230 | // Output list of clusters |
| 231 | // returns number of clusters |
| 232 | // if fit failed or unfolding is not applicable returns 0 and empty list |
| 233 | |
| 234 | //**************************** part 1 ******************************************* |
| 235 | // Performs the unfolding of a cluster with nMax overlapping showers |
| 236 | |
| 237 | //cout<<"unfolding check here part 1"<<endl; |
| 238 | AliDebug(5,Form(" Original cluster E %f, nMax = %d",iniTower->GetEnergy(),nMax )); |
| 239 | |
| 240 | Int_t nPar = 3 * nMax ; |
| 241 | Float_t * fitparameters = new Float_t[nPar] ; |
| 242 | |
| 243 | if (fGeom==0) |
| 244 | AliFatal("Did not get geometry from EMCALLoader") ; |
| 245 | |
| 246 | Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ; |
| 247 | if( !rv ) |
| 248 | { |
| 249 | // Fit failed, return (and remove cluster? - why? I leave the cluster) |
| 250 | iniTower->SetNExMax(-1) ; |
| 251 | delete[] fitparameters ; |
| 252 | return 0;//changed here |
| 253 | } |
| 254 | |
| 255 | //speed up solution for clusters with 2 maxima where one maximum is below threshold fThreshold |
| 256 | if(nMax==2){ |
| 257 | if(fitparameters[2]<fThreshold || fitparameters[5]<fThreshold){ |
| 258 | AliDebug(1,"One of fitted energy below threshold"); |
| 259 | iniTower->SetNExMax(1) ; |
| 260 | delete[] fitparameters ; |
| 261 | return 0;//changed here |
| 262 | } |
| 263 | } |
| 264 | |
| 265 | //**************************** part 2 ******************************************* |
| 266 | // create unfolded rec points and fill them with new energy lists |
| 267 | // First calculate energy deposited in each sell in accordance with |
| 268 | // fit (without fluctuations): efit[] |
| 269 | // and later correct this number in acordance with actual energy |
| 270 | // deposition |
| 271 | |
| 272 | // cout<<"unfolding check here part 2"<<endl; |
| 273 | Int_t nDigits = iniTower->GetMultiplicity() ; |
| 274 | Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells |
| 275 | Float_t xpar=0.,zpar=0.,epar=0. ;//center of gravity in cell units |
| 276 | |
| 277 | AliEMCALDigit * digit = 0 ; |
| 278 | Int_t * digitsList = iniTower->GetDigitsList() ; |
| 279 | |
| 280 | Int_t iSupMod = 0 ; |
| 281 | Int_t iTower = 0 ; |
| 282 | Int_t iIphi = 0 ; |
| 283 | Int_t iIeta = 0 ; |
| 284 | Int_t iphi = 0 ;//x direction |
| 285 | Int_t ieta = 0 ;//z direstion |
| 286 | |
| 287 | Int_t iparam = 0 ; |
| 288 | Int_t iDigit = 0 ; |
| 289 | |
| 290 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++) |
| 291 | { |
| 292 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ; |
| 293 | if(digit) |
| 294 | { |
| 295 | fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); |
| 296 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, |
| 297 | iIphi, iIeta,iphi,ieta); |
| 298 | EvalParsPhiDependence(digit->GetId(),fGeom); |
| 299 | |
| 300 | efit[iDigit] = 0.; |
| 301 | iparam = 0; |
| 302 | while(iparam < nPar ) |
| 303 | { |
| 304 | xpar = fitparameters[iparam] ; |
| 305 | zpar = fitparameters[iparam+1] ; |
| 306 | epar = fitparameters[iparam+2] ; |
| 307 | |
| 308 | efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; |
| 309 | iparam += 3 ; |
| 310 | } |
| 311 | |
| 312 | } else AliDebug(1,"Digit NULL part 2!"); |
| 313 | |
| 314 | }//digit loop |
| 315 | |
| 316 | //**************************** part 3 ******************************************* |
| 317 | // Now create new RecPoints and fill energy lists with efit corrected to fluctuations |
| 318 | // so that energy deposited in each cell is distributed between new clusters proportionally |
| 319 | // to its contribution to efit |
| 320 | |
| 321 | Float_t * energiesList = iniTower->GetEnergiesList() ; |
| 322 | Float_t ratio = 0. ; |
| 323 | Float_t eDigit = 0. ; |
| 324 | Int_t nSplittedClusters=(Int_t)nPar/3; |
| 325 | |
| 326 | Float_t * correctedEnergyList = new Float_t[nDigits*nSplittedClusters]; |
| 327 | //above - temporary table with energies after unfolding. |
| 328 | //the order is following: |
| 329 | //first cluster <first cell - last cell>, |
| 330 | //second cluster <first cell - last cell>, etc. |
| 331 | |
| 332 | //**************************** sub-part 3.1 ************************************* |
| 333 | //If not the energy from a given cell in the cluster is divided in correct proportions |
| 334 | //in accordance to the other clusters and added to them and set to 0. |
| 335 | |
| 336 | // cout<<"unfolding check here part 3.1"<<endl; |
| 337 | |
| 338 | iparam = 0 ; |
| 339 | while(iparam < nPar ) |
| 340 | { |
| 341 | xpar = fitparameters[iparam] ; |
| 342 | zpar = fitparameters[iparam+1] ; |
| 343 | epar = fitparameters[iparam+2] ; |
| 344 | |
| 345 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++) |
| 346 | { |
| 347 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ; |
| 348 | if(digit) |
| 349 | { |
| 350 | fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); |
| 351 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, |
| 352 | iIphi, iIeta,iphi,ieta); |
| 353 | |
| 354 | EvalParsPhiDependence(digit->GetId(),fGeom); |
| 355 | |
| 356 | if(efit[iDigit]==0) |
| 357 | {//just for sure |
| 358 | correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;//correction here |
| 359 | continue; |
| 360 | } |
| 361 | |
| 362 | ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ; |
| 363 | eDigit = energiesList[iDigit] * ratio ; |
| 364 | |
| 365 | //add energy to temporary matrix |
| 366 | correctedEnergyList[iparam/3*nDigits+iDigit] = eDigit; |
| 367 | |
| 368 | } else AliDebug(1,"NULL digit part 3"); |
| 369 | }//digit loop |
| 370 | iparam += 3 ; |
| 371 | }//while |
| 372 | |
| 373 | //**************************** sub-part 3.2 ************************************* |
| 374 | //here we check if energy of the cell in the cluster after unfolding is above threshold. |
| 375 | //here we correct energy for each cell and cluster |
| 376 | // cout<<"unfolding check here part 3.2"<<endl; |
| 377 | |
| 378 | |
| 379 | //here we have 3 possibilities |
| 380 | //when after UF cell energy in cluster is below threshold: |
| 381 | //1 - keep it associated to cluster - equivalent of threshold=0 |
| 382 | //2 - default - split (or add) energy of that cell into that cell in the other cluster(s) |
| 383 | //3 - reject that cell from cluster - fraction energy in cell=0 - breaks energy conservation |
| 384 | //Bool_t rejectBelowThreshold=kTRUE;//default option = 2 - split = kFALSE |
| 385 | |
| 386 | if(fThreshold > 0){//option 2 or 3 |
| 387 | if(fRejectBelowThreshold){//option 3 |
| 388 | for(iDigit = 0 ; iDigit < nDigits ; iDigit++){//digit loop |
| 389 | for(iparam = 0 ; iparam < nPar ; iparam+=3){//param0 loop = energy loop |
| 390 | if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold ) correctedEnergyList[iparam/3*nDigits+iDigit]=0.; |
| 391 | } |
| 392 | } |
| 393 | }else{//option 2 |
| 394 | Float_t maximumEne=0.; |
| 395 | Int_t maximumIndex=0; |
| 396 | Bool_t isAnyBelowThreshold=kFALSE; |
| 397 | // Float_t Threshold=0.01; |
| 398 | Float_t * energyFraction = new Float_t[nSplittedClusters]; |
| 399 | Int_t iparam2 = 0 ; |
| 400 | for(iDigit = 0 ; iDigit < nDigits ; iDigit++){ |
| 401 | isAnyBelowThreshold=kFALSE; |
| 402 | maximumEne=0.; |
| 403 | for(iparam = 0 ; iparam < nPar ; iparam+=3){ |
| 404 | if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold ) isAnyBelowThreshold = kTRUE; |
| 405 | if(correctedEnergyList[iparam/3*nDigits+iDigit] > maximumEne) |
| 406 | { |
| 407 | maximumEne = correctedEnergyList[iparam/3*nDigits+iDigit]; |
| 408 | maximumIndex = iparam; |
| 409 | } |
| 410 | }//end of loop over clusters after unfolding |
| 411 | |
| 412 | if(!isAnyBelowThreshold) continue; //no cluster-cell below threshold |
| 413 | |
| 414 | if(maximumEne < fThreshold) |
| 415 | {//add all cluster cells and put energy into max index, other set to 0 |
| 416 | maximumEne=0.; |
| 417 | for(iparam = 0 ; iparam < nPar ; iparam+=3) |
| 418 | { |
| 419 | maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit]; |
| 420 | correctedEnergyList[iparam/3*nDigits+iDigit]=0; |
| 421 | } |
| 422 | correctedEnergyList[maximumIndex/3*nDigits+iDigit]=maximumEne; |
| 423 | continue; |
| 424 | }//end if |
| 425 | |
| 426 | //divide energy of cell below threshold in the correct proportion and add to other cells |
| 427 | maximumEne=0.;//not used any more so use it for the energy sum |
| 428 | for(iparam = 0 ; iparam < nPar ; iparam+=3) |
| 429 | {//calculate energy sum |
| 430 | if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold) energyFraction[iparam/3]=0; |
| 431 | else |
| 432 | { |
| 433 | energyFraction[iparam/3]=1.; |
| 434 | maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit]; |
| 435 | } |
| 436 | }//end of loop over clusters after unfolding |
| 437 | if(maximumEne>0.) { |
| 438 | for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction |
| 439 | energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3*nDigits+iDigit] / maximumEne; |
| 440 | } |
| 441 | |
| 442 | for(iparam = 0 ; iparam < nPar ; iparam+=3) |
| 443 | {//add energy from cells below threshold to others |
| 444 | if(energyFraction[iparam/3]>0.) continue; |
| 445 | else |
| 446 | { |
| 447 | for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3) |
| 448 | { |
| 449 | correctedEnergyList[iparam2/3*nDigits+iDigit] += (energyFraction[iparam2/3] * |
| 450 | correctedEnergyList[iparam/3*nDigits+iDigit]) ; |
| 451 | }//inner loop |
| 452 | correctedEnergyList[iparam/3*nDigits+iDigit] = 0.; |
| 453 | } |
| 454 | } |
| 455 | } else { |
| 456 | //digit energy to be set to 0 |
| 457 | for(iparam = 0 ; iparam < nPar ; iparam+=3) |
| 458 | { |
| 459 | correctedEnergyList[iparam/3*nDigits+iDigit] = 0.; |
| 460 | } |
| 461 | }//correction for: is energy>0 |
| 462 | |
| 463 | }//end of loop over digits |
| 464 | delete[] energyFraction; |
| 465 | |
| 466 | }//end of option 2 or 3 |
| 467 | } else {//option 1 |
| 468 | //do nothing |
| 469 | } |
| 470 | |
| 471 | |
| 472 | //**************************** sub-part 3.3 ************************************* |
| 473 | //here we add digits to recpoints with corrected energy |
| 474 | // cout<<"unfolding check here part 3.3"<<endl; |
| 475 | |
| 476 | Int_t newClusterIndex=0; |
| 477 | iparam = 0 ; |
| 478 | while(iparam < nPar ) |
| 479 | { |
| 480 | AliEMCALRecPoint * recPoint = 0 ; |
| 481 | |
| 482 | if(nSplittedClusters >= list->GetSize()) |
| 483 | list->Expand(nSplittedClusters); |
| 484 | |
| 485 | //add recpoint |
| 486 | (*list)[newClusterIndex] = new AliEMCALRecPoint("") ; |
| 487 | recPoint = dynamic_cast<AliEMCALRecPoint *>( list->At(newClusterIndex) ) ; |
| 488 | |
| 489 | if(recPoint){//recPoint present -> good |
| 490 | recPoint->SetNExMax(nSplittedClusters) ;//can be wrong number, to be corrected in outer method |
| 491 | |
| 492 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++) { |
| 493 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ; |
| 494 | if(digit && correctedEnergyList[iparam/3*nDigits+iDigit]>0. ){ |
| 495 | //if(correctedEnergyList[iparam/3*nDigits+iDigit]<fThreshold) printf("Final E cell %f < %f\n",correctedEnergyList[iparam/3*nDigits+iDigit],fThreshold); |
| 496 | recPoint->AddDigit( *digit, correctedEnergyList[iparam/3*nDigits+iDigit], kFALSE ) ; //FIXME, need to study the shared case |
| 497 | } else { |
| 498 | AliDebug(1,Form("NULL digit part3.3 or NULL energy=%f",correctedEnergyList[iparam/3*nDigits+iDigit])); |
| 499 | } |
| 500 | }//digit loop |
| 501 | |
| 502 | if(recPoint->GetMultiplicity()==0){//recpoint exists but no digits associated -> remove from list |
| 503 | delete (*list)[newClusterIndex]; |
| 504 | list->RemoveAt(newClusterIndex); |
| 505 | nSplittedClusters--; |
| 506 | newClusterIndex--;//decrease cluster number |
| 507 | }else {//recPoint exists and has digits associated -> very good increase number of clusters |
| 508 | AliDebug(5,Form("cluster %d, digit no %d, energy %f",iparam/3,(recPoint->GetDigitsList())[0],(recPoint->GetEnergiesList())[0])); |
| 509 | } |
| 510 | |
| 511 | } else {//recPoint empty -> remove from list |
| 512 | AliError("NULL RecPoint"); |
| 513 | //protection from recpoint with no digits |
| 514 | delete (*list)[newClusterIndex]; |
| 515 | list->RemoveAt(newClusterIndex); |
| 516 | nSplittedClusters--; |
| 517 | newClusterIndex--;//decrease cluster number |
| 518 | } |
| 519 | |
| 520 | iparam += 3 ; |
| 521 | newClusterIndex++; |
| 522 | }//while |
| 523 | |
| 524 | delete[] fitparameters ; |
| 525 | delete[] efit ; |
| 526 | delete[] correctedEnergyList ; |
| 527 | |
| 528 | // print |
| 529 | AliDebug(5,Form(" nSplittedClusters %d, fNumberOfECAClusters %d, newClusterIndex %d,list->Entries() %d\n",nSplittedClusters,fNumberOfECAClusters,newClusterIndex,list->GetEntriesFast() )); |
| 530 | |
| 531 | // cout<<"end of unfolding check part 3.3"<<endl; |
| 532 | return nSplittedClusters; |
| 533 | } |
| 534 | |
| 535 | //____________________________________________________________________________ |
| 536 | Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower, |
| 537 | Int_t nMax, |
| 538 | AliEMCALDigit ** maxAt, |
| 539 | Float_t * maxAtEnergy) |
| 540 | { |
| 541 | // Extended to whole EMCAL |
| 542 | // Performs the unfolding of a cluster with nMax overlapping showers |
| 543 | // Returns true if success (1->several clusters), otherwise false (fit failed) |
| 544 | |
| 545 | TObjArray *list =new TObjArray(2);//temporary object |
| 546 | Int_t nUnfoldedClusters=UnfoldOneCluster(iniTower,nMax,maxAt,maxAtEnergy,list); |
| 547 | |
| 548 | // here we write new clusters from list to fRecPoints |
| 549 | AliDebug(5,Form("Number of clusters after unfolding %d",list->GetEntriesFast())); |
| 550 | Int_t iDigit=0; |
| 551 | AliEMCALDigit * digit = 0 ; |
| 552 | for(Int_t i=0;i<list->GetEntriesFast();i++) { |
| 553 | AliEMCALRecPoint * recPoint = 0 ; |
| 554 | |
| 555 | if(fNumberOfECAClusters >= fRecPoints->GetSize()) |
| 556 | fRecPoints->Expand(2*fNumberOfECAClusters) ; |
| 557 | |
| 558 | //add recpoint |
| 559 | (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;//fNumberOfECAClusters-1 is old cluster before unfolding |
| 560 | recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ; |
| 868ccccd |
561 | AliEMCALRecPoint * rpUFOne = dynamic_cast<AliEMCALRecPoint *>(list->At(i)) ; |
| 562 | |
| 563 | if( recPoint && rpUFOne ){//recPoint present -> good |
| 564 | |
| 565 | recPoint->SetNExMax(list->GetEntriesFast()) ; |
| 566 | |
| 567 | Int_t *digitsList = rpUFOne->GetDigitsList(); |
| 568 | Float_t *energyList = rpUFOne->GetEnergiesList(); |
| b410dc6a |
569 | |
| 78a3ea18 |
570 | if(!digitsList || ! energyList) |
| 571 | { |
| 572 | AliDebug(-1,"No digits index or energy available"); |
| 868ccccd |
573 | delete (*fRecPoints)[fNumberOfECAClusters]; |
| 574 | fRecPoints->RemoveAt(fNumberOfECAClusters); |
| 78a3ea18 |
575 | continue; |
| 576 | } |
| 577 | |
| b410dc6a |
578 | AliDebug(5,Form("cluster %d, digit no %d, energy %f\n",i,digitsList[0],energyList[0])); |
| 579 | |
| 868ccccd |
580 | for(iDigit = 0 ; iDigit < rpUFOne->GetMultiplicity(); iDigit ++) { |
| b410dc6a |
581 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ; |
| e4e78a4d |
582 | if(digit) recPoint->AddDigit( *digit, energyList[iDigit], kFALSE ) ; //FIXME, need to study the shared case |
| b410dc6a |
583 | }//digit loop |
| 584 | fNumberOfECAClusters++ ; |
| 585 | } else {//recPoint empty -> remove from list |
| 586 | AliError("NULL RecPoint"); |
| 587 | delete (*fRecPoints)[fNumberOfECAClusters]; |
| 588 | fRecPoints->RemoveAt(fNumberOfECAClusters); |
| 589 | } |
| 590 | |
| 591 | }//loop over unfolded clusters |
| 592 | |
| 593 | //print energy of new unfolded clusters |
| 594 | AliDebug(5,Form(" nUnfoldedClusters %d, fNumberOfECAClusters %d",nUnfoldedClusters,fNumberOfECAClusters )); |
| 595 | for(Int_t inewclus=0; inewclus<nUnfoldedClusters;inewclus++){ |
| e4e78a4d |
596 | AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters-1-inewclus)); |
| 597 | if(rp) AliDebug(5,Form(" Unfolded cluster %d E %f",inewclus, rp->GetEnergy() )); |
| b410dc6a |
598 | } |
| 599 | |
| 600 | //clear tables |
| 601 | list->SetOwner(kTRUE); |
| 602 | list->Delete(); |
| 603 | delete list; |
| 604 | if(nUnfoldedClusters>1) return kTRUE; |
| 605 | return kFALSE; |
| 606 | } |
| 607 | |
| 608 | |
| 609 | |
| 610 | //____________________________________________________________________________ |
| 611 | Bool_t AliEMCALUnfolding::UnfoldClusterV2old(AliEMCALRecPoint * iniTower, |
| 612 | Int_t nMax, |
| 613 | AliEMCALDigit ** maxAt, |
| 614 | Float_t * maxAtEnergy) |
| 615 | { |
| 616 | // Extended to whole EMCAL |
| 617 | // Performs the unfolding of a cluster with nMax overlapping showers |
| 618 | |
| 619 | Int_t nPar = 3 * nMax ; |
| 620 | Float_t * fitparameters = new Float_t[nPar] ; |
| 621 | |
| 622 | if (fGeom==0) |
| 623 | AliFatal("Did not get geometry from EMCALLoader") ; |
| 624 | |
| 625 | Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ; |
| 626 | if( !rv ) { |
| 627 | // Fit failed, return (and remove cluster? - why? I leave the cluster) |
| 628 | iniTower->SetNExMax(-1) ; |
| 629 | delete[] fitparameters ; |
| 630 | return kFALSE; |
| 631 | } |
| 632 | |
| 633 | // create unfolded rec points and fill them with new energy lists |
| 634 | // First calculate energy deposited in each sell in accordance with |
| 635 | // fit (without fluctuations): efit[] |
| 636 | // and later correct this number in acordance with actual energy |
| 637 | // deposition |
| 638 | |
| 639 | Int_t nDigits = iniTower->GetMultiplicity() ; |
| 640 | Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells |
| 641 | Float_t xpar=0.,zpar=0.,epar=0. ;//center of gravity in cell units |
| 642 | |
| 643 | AliEMCALDigit * digit = 0 ; |
| 644 | Int_t * digitsList = iniTower->GetDigitsList() ; |
| 645 | |
| 646 | Int_t iSupMod = 0 ; |
| 647 | Int_t iTower = 0 ; |
| 648 | Int_t iIphi = 0 ; |
| 649 | Int_t iIeta = 0 ; |
| 650 | Int_t iphi = 0 ;//x direction |
| 651 | Int_t ieta = 0 ;//z direstion |
| 652 | |
| 653 | Int_t iparam = 0 ; |
| 654 | Int_t iDigit = 0 ; |
| 655 | |
| 656 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ |
| 657 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ; |
| 658 | if(digit){ |
| 659 | fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); |
| 660 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, |
| 661 | iIphi, iIeta,iphi,ieta); |
| 662 | EvalParsPhiDependence(digit->GetId(),fGeom); |
| 663 | |
| 664 | efit[iDigit] = 0.; |
| 665 | iparam = 0; |
| 666 | while(iparam < nPar ){ |
| 667 | xpar = fitparameters[iparam] ; |
| 668 | zpar = fitparameters[iparam+1] ; |
| 669 | epar = fitparameters[iparam+2] ; |
| 670 | iparam += 3 ; |
| 671 | |
| 672 | efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; |
| 673 | } |
| 674 | } else AliError("Digit NULL!"); |
| 675 | |
| 676 | }//digit loop |
| 677 | |
| 678 | // Now create new RecPoints and fill energy lists with efit corrected to fluctuations |
| 679 | // so that energy deposited in each cell is distributed between new clusters proportionally |
| 680 | // to its contribution to efit |
| 681 | |
| 682 | Float_t * energiesList = iniTower->GetEnergiesList() ; |
| 683 | Float_t ratio = 0 ; |
| 684 | |
| 685 | iparam = 0 ; |
| 686 | while(iparam < nPar ){ |
| 687 | xpar = fitparameters[iparam] ; |
| 688 | zpar = fitparameters[iparam+1] ; |
| 689 | epar = fitparameters[iparam+2] ; |
| 690 | iparam += 3 ; |
| 691 | |
| 692 | AliEMCALRecPoint * recPoint = 0 ; |
| 693 | |
| 694 | if(fNumberOfECAClusters >= fRecPoints->GetSize()) |
| 695 | fRecPoints->Expand(2*fNumberOfECAClusters) ; |
| 696 | |
| 697 | //add recpoint |
| 698 | (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ; |
| 699 | recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ; |
| 700 | |
| 701 | if(recPoint){ |
| 702 | |
| 703 | fNumberOfECAClusters++ ; |
| 704 | recPoint->SetNExMax((Int_t)nPar/3) ; |
| 705 | |
| 706 | Float_t eDigit = 0. ; |
| 707 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ |
| 708 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ; |
| 709 | if(digit){ |
| 710 | fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); |
| 711 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, |
| 712 | iIphi, iIeta,iphi,ieta); |
| 713 | EvalParsPhiDependence(digit->GetId(),fGeom); |
| 714 | if(efit[iDigit]==0) continue;//just for sure |
| 715 | ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ; |
| 716 | eDigit = energiesList[iDigit] * ratio ; |
| 717 | recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case |
| 718 | } else AliError("NULL digit"); |
| 719 | }//digit loop |
| 720 | } else AliError("NULL RecPoint"); |
| 721 | }//while |
| 722 | |
| 723 | delete[] fitparameters ; |
| 724 | delete[] efit ; |
| 725 | |
| 726 | return kTRUE; |
| 727 | } |
| 728 | |
| 729 | |
| 730 | //____________________________________________________________________________ |
| 731 | Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt, |
| 732 | const Float_t* maxAtEnergy, |
| 733 | Int_t nPar, Float_t * fitparameters) const |
| 734 | { |
| 735 | // Calls TMinuit to fit the energy distribution of a cluster with several maxima |
| 736 | // The initial values for fitting procedure are set equal to the |
| 737 | // positions of local maxima. |
| 738 | // Cluster will be fitted as a superposition of nPar/3 |
| 739 | // electromagnetic showers |
| 740 | |
| 741 | if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader"); |
| 742 | |
| 743 | if(!gMinuit){ |
| 744 | // gMinuit = new TMinuit(100) ;//max 100 parameters |
| 745 | if(nPar<30) gMinuit = new TMinuit(30); |
| 746 | else gMinuit = new TMinuit(nPar) ;//max nPar parameters |
| 747 | // |
| 748 | } else { |
| 749 | if(gMinuit->fMaxpar < nPar) { |
| 750 | delete gMinuit; |
| 751 | gMinuit = new TMinuit(nPar); |
| 752 | } |
| 753 | } |
| 754 | |
| 755 | gMinuit->mncler(); // Reset Minuit's list of paramters |
| 756 | gMinuit->SetPrintLevel(-1) ; // No Printout |
| 757 | gMinuit->SetFCN(AliEMCALUnfolding::UnfoldingChiSquareV2) ; |
| 758 | // To set the address of the minimization function |
| 759 | TList * toMinuit = new TList(); |
| 760 | toMinuit->AddAt(recPoint,0) ; |
| 761 | toMinuit->AddAt(fDigitsArr,1) ; |
| 762 | toMinuit->AddAt(fGeom,2) ; |
| 763 | |
| 764 | gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare |
| 765 | |
| 766 | // filling initial values for fit parameters |
| 767 | AliEMCALDigit * digit ; |
| 768 | |
| 769 | Int_t ierflg = 0; |
| 770 | Int_t index = 0 ; |
| 771 | Int_t nDigits = (Int_t) nPar / 3 ; |
| 772 | |
| 773 | Int_t iDigit ; |
| 774 | |
| 775 | Int_t iSupMod = 0 ; |
| 776 | Int_t iTower = 0 ; |
| 777 | Int_t iIphi = 0 ; |
| 778 | Int_t iIeta = 0 ; |
| 779 | Int_t iphi = 0 ;//x direction |
| 780 | Int_t ieta = 0 ;//z direstion |
| 781 | |
| 782 | for(iDigit = 0; iDigit < nDigits; iDigit++){ |
| 783 | digit = maxAt[iDigit]; |
| 784 | if(digit==0) AliError("energy of digit = 0!"); |
| 785 | fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); |
| 786 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, |
| 787 | iIphi, iIeta,iphi,ieta); |
| 788 | |
| 789 | Float_t energy = maxAtEnergy[iDigit] ; |
| 790 | |
| 791 | //gMinuit->mnparm(index, "x", iphi, 0.1, 0, 0, ierflg) ;//original |
| 792 | gMinuit->mnparm(index, "x", iphi, 0.05, 0, 0, ierflg) ; |
| 793 | index++ ; |
| 794 | if(ierflg != 0){ |
| 795 | Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: x=%d, param.id=%d, nMaxima=%d",iphi,index-1,nPar/3 ) ; |
| 796 | toMinuit->Clear(); |
| 797 | delete toMinuit ; |
| 798 | return kFALSE; |
| 799 | } |
| 800 | //gMinuit->mnparm(index, "z", ieta, 0.1, 0, 0, ierflg) ;//original |
| 801 | gMinuit->mnparm(index, "z", ieta, 0.05, 0, 0, ierflg) ; |
| 802 | index++ ; |
| 803 | if(ierflg != 0){ |
| 804 | Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: z=%d, param.id=%d, nMaxima=%d", ieta, index-1,nPar/3) ; |
| 805 | toMinuit->Clear(); |
| 806 | delete toMinuit ; |
| 807 | return kFALSE; |
| 808 | } |
| 809 | //gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;//original |
| 810 | gMinuit->mnparm(index, "Energy", energy , 0.001*energy, 0., 5.*energy, ierflg) ;//was 0.05 |
| 811 | index++ ; |
| 812 | if(ierflg != 0){ |
| 813 | Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: energy = %f, param.id=%d, nMaxima=%d", energy, index-1, nPar/3) ; |
| 814 | toMinuit->Clear(); |
| 815 | delete toMinuit ; |
| 816 | return kFALSE; |
| 817 | } |
| 818 | } |
| 819 | |
| 820 | Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; |
| 821 | // The number of function call slightly depends on it. |
| 822 | // Double_t p1 = 1.0 ;// par to gradient |
| 823 | Double_t p2 = 0.0 ; |
| 824 | // Double_t p3 = 3.0 ; |
| 825 | gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls |
| 826 | // gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient |
| 827 | gMinuit->SetMaxIterations(5);//was 5 |
| 828 | gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings |
| 829 | //gMinuit->mnexcm("SET PRI", &p3 , 3, ierflg) ; // printouts |
| 830 | |
| 831 | gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize |
| 832 | //gMinuit->mnexcm("MINI", &p0, 0, ierflg) ; // minimize |
| 833 | if(ierflg == 4){ // Minimum not found |
| 834 | AliDebug(1,"EMCAL Unfolding Fit not converged, cluster abandoned " ) ; |
| 835 | toMinuit->Clear(); |
| 836 | delete toMinuit ; |
| 837 | return kFALSE ; |
| 838 | } |
| 839 | for(index = 0; index < nPar; index++){ |
| 840 | Double_t err = 0. ; |
| 841 | Double_t val = 0. ; |
| 842 | gMinuit->GetParameter(index, val, err) ; // Returns value and error ofOA parameter index |
| 843 | fitparameters[index] = val ; |
| 844 | } |
| 845 | |
| 846 | toMinuit->Clear(); |
| 847 | delete toMinuit ; |
| 848 | |
| 849 | if(gMinuit->fMaxpar>30) delete gMinuit; |
| 850 | |
| 851 | return kTRUE; |
| 852 | |
| 853 | } |
| 854 | |
| 855 | //____________________________________________________________________________ |
| 856 | Double_t AliEMCALUnfolding::ShowerShapeV2(Double_t x, Double_t y) |
| 857 | { |
| 858 | // extended to whole EMCAL |
| 859 | // Shape of the shower |
| 860 | // If you change this function, change also the gradient evaluation in ChiSquare() |
| 861 | |
| 862 | Double_t r = fgSSPars[7]*TMath::Sqrt(x*x+y*y); |
| 863 | Double_t rp1 = TMath::Power(r, fgSSPars[1]) ; |
| 864 | Double_t rp5 = TMath::Power(r, fgSSPars[5]) ; |
| 865 | Double_t shape = fgSSPars[0]*TMath::Exp( -rp1 * (1. / (fgSSPars[2] + fgSSPars[3] * rp1) + fgSSPars[4] / (1 + fgSSPars[6] * rp5) ) ) ; |
| 866 | return shape ; |
| 867 | } |
| 868 | |
| 869 | //____________________________________________________________________________ |
| 870 | void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad, |
| 871 | Double_t & fret, |
| 872 | Double_t * x, Int_t iflag) |
| 873 | { |
| 874 | // Calculates the Chi square for the cluster unfolding minimization |
| 875 | // Number of parameters, Gradient, Chi squared, parameters, what to do |
| 876 | |
| 877 | TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ; |
| 878 | if(toMinuit){ |
| 879 | AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) ) ; |
| 880 | TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ; |
| 881 | // A bit buggy way to get an access to the geometry |
| 882 | // To be revised! |
| 883 | AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2)); |
| 884 | |
| 885 | if(recPoint && digits && geom){ |
| 886 | |
| 887 | Int_t * digitsList = recPoint->GetDigitsList() ; |
| 888 | |
| 889 | Int_t nOdigits = recPoint->GetDigitsMultiplicity() ; |
| 890 | |
| 891 | Float_t * energiesList = recPoint->GetEnergiesList() ; |
| 892 | |
| 893 | fret = 0. ; |
| 894 | Int_t iparam = 0 ; |
| 895 | |
| 896 | if(iflag == 2) |
| 897 | for(iparam = 0 ; iparam < nPar ; iparam++) |
| 898 | Grad[iparam] = 0 ; // Will evaluate gradient |
| 899 | |
| 900 | Double_t efit = 0. ; |
| 901 | |
| 902 | AliEMCALDigit * digit ; |
| 903 | Int_t iDigit ; |
| 904 | |
| 905 | Int_t iSupMod = 0 ; |
| 906 | Int_t iTower = 0 ; |
| 907 | Int_t iIphi = 0 ; |
| 908 | Int_t iIeta = 0 ; |
| 909 | Int_t iphi = 0 ;//x direction |
| 910 | Int_t ieta = 0 ;//z direstion |
| 911 | |
| 912 | |
| 913 | for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) { |
| 914 | if(energiesList[iDigit]==0) continue; |
| 915 | |
| 916 | digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) ); |
| 917 | |
| 918 | if(digit){ |
| 919 | geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); |
| 920 | geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, |
| 921 | iIphi, iIeta,iphi,ieta); |
| 922 | EvalParsPhiDependence(digit->GetId(),geom); |
| 923 | |
| 924 | if(iflag == 2){ // calculate gradient |
| 925 | Int_t iParam = 0 ; |
| 926 | efit = 0. ; |
| 927 | while(iParam < nPar ){ |
| 928 | Double_t dx = ((Float_t)iphi - x[iParam]) ; |
| 929 | iParam++ ; |
| 930 | Double_t dz = ((Float_t)ieta - x[iParam]) ; |
| 931 | iParam++ ; |
| 932 | efit += x[iParam] * ShowerShapeV2(dx,dz) ; |
| 933 | iParam++ ; |
| 934 | } |
| 935 | |
| 936 | Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E) |
| 937 | iParam = 0 ; |
| 938 | while(iParam < nPar ){ |
| 939 | Double_t xpar = x[iParam] ; |
| 940 | Double_t zpar = x[iParam+1] ; |
| 941 | Double_t epar = x[iParam+2] ; |
| 942 | |
| 943 | Double_t dr = fgSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) ); |
| 944 | Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; |
| 945 | Double_t rp1 = TMath::Power(dr, fgSSPars[1]) ; |
| 946 | Double_t rp5 = TMath::Power(dr, fgSSPars[5]) ; |
| 947 | |
| 948 | Double_t deriv = -2 * TMath::Power(dr,fgSSPars[1]-2.) * fgSSPars[7] * fgSSPars[7] * |
| 949 | (fgSSPars[1] * ( 1/(fgSSPars[2]+fgSSPars[3]*rp1) + fgSSPars[4]/(1+fgSSPars[6]*rp5) ) - |
| 950 | (fgSSPars[1]*fgSSPars[3]*rp1/( (fgSSPars[2]+fgSSPars[3]*rp1)*(fgSSPars[2]+fgSSPars[3]*rp1) ) + |
| 951 | fgSSPars[4]*fgSSPars[5]*fgSSPars[6]*rp5/( (1+fgSSPars[6]*rp5)*(1+fgSSPars[6]*rp5) ) ) ); |
| 952 | |
| 953 | //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) ) |
| 954 | // - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ; |
| 955 | |
| 956 | Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ; // Derivative over x |
| 957 | iParam++ ; |
| 958 | Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ; // Derivative over z |
| 959 | iParam++ ; |
| 960 | Grad[iParam] += shape ; // Derivative over energy |
| 961 | iParam++ ; |
| 962 | } |
| 963 | } |
| 964 | efit = 0; |
| 965 | iparam = 0 ; |
| 966 | |
| 967 | while(iparam < nPar ){ |
| 968 | Double_t xpar = x[iparam] ; |
| 969 | Double_t zpar = x[iparam+1] ; |
| 970 | Double_t epar = x[iparam+2] ; |
| 971 | iparam += 3 ; |
| 972 | efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; |
| 973 | } |
| 974 | |
| 975 | fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ; |
| 976 | // Here we assume, that sigma = sqrt(E) |
| 977 | } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!, nPar %d \n", nPar); // put nPar here to cheat coverity and rule checker |
| 978 | } // digit loop |
| 979 | } // recpoint, digits and geom not NULL |
| 980 | }// List is not NULL |
| 981 | |
| 982 | } |
| 983 | |
| 984 | |
| 985 | //____________________________________________________________________________ |
| 986 | void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){ |
| 987 | for(UInt_t i=0;i<7;++i) |
| 988 | fgSSPars[i]=pars[i]; |
| 989 | if(pars[2]==0. && pars[3]==0.) fgSSPars[2]=1.;//to avoid dividing by 0 |
| 990 | } |
| 991 | |
| 992 | //____________________________________________________________________________ |
| 993 | void AliEMCALUnfolding::SetPar5(Double_t *pars){ |
| 994 | for(UInt_t i=0;i<3;++i) |
| 995 | fgPar5[i]=pars[i]; |
| 996 | } |
| 997 | |
| 998 | //____________________________________________________________________________ |
| 999 | void AliEMCALUnfolding::SetPar6(Double_t *pars){ |
| 1000 | for(UInt_t i=0;i<3;++i) |
| 1001 | fgPar6[i]=pars[i]; |
| 1002 | } |
| 1003 | |
| 1004 | //____________________________________________________________________________ |
| 1005 | void AliEMCALUnfolding::EvalPar5(Double_t phi){ |
| 1006 | // |
| 1007 | //Evaluate the 5th parameter of the shower shape function |
| 1008 | //phi in degrees range (-10,10) |
| 1009 | // |
| 1010 | //fSSPars[5] = 12.31 - phi*0.007381 - phi*phi*0.06936; |
| 1011 | fgSSPars[5] = fgPar5[0] + phi * fgPar5[1] + phi*phi * fgPar5[2]; |
| 1012 | } |
| 1013 | |
| 1014 | //____________________________________________________________________________ |
| 1015 | void AliEMCALUnfolding::EvalPar6(Double_t phi){ |
| 1016 | // |
| 1017 | //Evaluate the 6th parameter of the shower shape function |
| 1018 | //phi in degrees range (-10,10) |
| 1019 | // |
| 1020 | //fSSPars[6] = 0.05452 + phi*0.0001228 + phi*phi*0.001361; |
| 1021 | fgSSPars[6] = fgPar6[0] + phi * fgPar6[1] + phi*phi * fgPar6[2]; |
| 1022 | } |
| 1023 | |
| 1024 | //____________________________________________________________________________ |
| 1025 | void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, const AliEMCALGeometry *geom){ |
| 1026 | // |
| 1027 | // calculate params p5 and p6 depending on the phi angle in global coordinate |
| 1028 | // for the cell with given absId index |
| 1029 | // |
| 1030 | Double_t etaGlob = 0.;//eta in global c.s. - unused |
| 1031 | Double_t phiGlob = 0.;//phi in global c.s. in radians |
| 1032 | geom->EtaPhiFromIndex(absId, etaGlob, phiGlob); |
| 1033 | phiGlob*=180./TMath::Pi(); |
| 1034 | phiGlob-=90.; |
| 1035 | phiGlob-= (Double_t)((Int_t)geom->GetSuperModuleNumber(absId)/2 * 20); |
| 1036 | |
| 1037 | EvalPar5(phiGlob); |
| 1038 | EvalPar6(phiGlob); |
| 1039 | } |
| 1040 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame