New classes and test macro for significance maximisation in multi-dimensional cuts...
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliMultiDimVector.cxx
1 /**************************************************************************
2  * Copyright(c) 2007-2009, 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 /* $Id: $ */
17
18 ///////////////////////////////////////////////////////////////////
19 //                                                               //
20 // Implementation of the class to store the number of signal     //
21 // and background events in bins of the cut values               //
22 // Origin:       Elena Bruna (bruna@to.infn.it)                  //
23 // Updated:      Sergey Senyukov (senyukov@to.infn.it)           //
24 // Last updated: Francesco Prino (prino@to.infn.it)              //
25 //                                                               //
26 ///////////////////////////////////////////////////////////////////
27
28
29 #include "TH2.h"
30 #include "AliMultiDimVector.h"
31 #include "AliLog.h"
32 #include "TMath.h"
33 #include "TString.h"
34
35 ClassImp(AliMultiDimVector)
36 //___________________________________________________________________________
37 AliMultiDimVector::AliMultiDimVector():TNamed("AliMultiDimVector","default"),
38 fNVariables(0),
39 fNPtBins(0),
40 fVett(0),
41 fNTotCells(0),
42 fIsIntegrated(0)
43 {
44   // default constructor
45 }
46 //___________________________________________________________________________
47 AliMultiDimVector::AliMultiDimVector(const char *name,const char *title, const Int_t npars, Int_t nptbins, Int_t *nofcells, Float_t *loosecuts, Float_t *tightcuts, TString *axisTitles):TNamed(name,title),
48 fNVariables(npars),
49 fNPtBins(nptbins),
50 fVett(0),
51 fNTotCells(0),
52 fIsIntegrated(0){
53 // standard constructor
54   ULong64_t ntot=1; 
55   for(Int_t i=0;i<fNVariables;i++){
56     ntot*=nofcells[i];
57     fNCutSteps[i]=nofcells[i];
58     if(loosecuts[i] <= tightcuts[i]){
59       fMinLimits[i]=loosecuts[i];
60       fMaxLimits[i]=tightcuts[i];
61       fGreaterThan[i]=kTRUE;
62     }else{
63       fMinLimits[i]=tightcuts[i];
64       fMaxLimits[i]=loosecuts[i];
65       fGreaterThan[i]=kFALSE;
66     }
67     fAxisTitles[i]=axisTitles[i].Data();
68   }
69   fNTotCells=ntot*fNPtBins;
70   fVett.Set(fNTotCells); 
71   for (ULong64_t j=0;j<fNTotCells;j++) fVett.AddAt(0,j);
72 }
73 //___________________________________________________________________________
74 AliMultiDimVector::AliMultiDimVector(const AliMultiDimVector &mv):TNamed(mv.GetName(),mv.GetTitle()),
75 fNVariables(mv.fNVariables),
76 fNPtBins(mv.fNPtBins),
77 fVett(0),
78 fNTotCells(mv.fNTotCells),
79 fIsIntegrated(mv.fIsIntegrated)
80 {
81   // copy constructor
82   for(Int_t i=0;i<fNVariables;i++){
83     fNCutSteps[i]=mv.GetNCutSteps(i);
84     fMinLimits[i]=mv.GetMinLimit(i);
85     fMaxLimits[i]=mv.GetMaxLimit(i);
86     fGreaterThan[i]=mv.GetGreaterThan(i);
87     fAxisTitles[i]=mv.GetAxisTitle(i);
88   }
89   fVett.Set(fNTotCells); 
90   for(ULong64_t i=0;i<fNTotCells;i++) fVett[i]=mv.GetElement(i);
91 }
92 //___________________________________________________________________________
93 void AliMultiDimVector::CopyStructure(const AliMultiDimVector* mv){
94 // Sets dimensions and limit from mv
95   fNVariables=mv->GetNVariables();
96   fNPtBins=mv->GetNPtBins();
97   fNTotCells=mv->GetNTotCells();
98   fIsIntegrated=mv->IsIntegrated();
99   for(Int_t i=0;i<fNVariables;i++){
100     fNCutSteps[i]=mv->GetNCutSteps(i);
101     fMinLimits[i]=mv->GetMinLimit(i);
102     fMaxLimits[i]=mv->GetMaxLimit(i);
103     fGreaterThan[i]=mv->GetGreaterThan(i);
104     fAxisTitles[i]=mv->GetAxisTitle(i);
105   }
106   fVett.Set(fNTotCells);
107   
108 }
109 //______________________________________________________________________
110 Bool_t AliMultiDimVector::GetIndicesFromGlobalAddress(ULong64_t globadd, Int_t *ind, Int_t &ptbin) const {
111 // returns matrix element indices and Pt bin from global index
112   if(globadd>=fNTotCells) return kFALSE;
113   ULong64_t r=globadd;
114   Int_t prod=1;
115   Int_t nOfCellsPlusLevel[fgkMaxNVariables+1];
116   for(Int_t k=0;k<fNVariables;k++) nOfCellsPlusLevel[k]=fNCutSteps[k];
117   nOfCellsPlusLevel[fNVariables]=fNPtBins;
118         
119   for(Int_t i=0;i<fNVariables+1;i++) prod*=nOfCellsPlusLevel[i];
120   for(Int_t i=0;i<fNVariables+1;i++){
121     prod/=nOfCellsPlusLevel[i];
122     if(i<fNVariables) ind[i]=r/prod;
123     else ptbin=r/prod;
124     r=globadd%prod;
125   }
126   return kTRUE;
127 }
128 //______________________________________________________________________
129 Bool_t AliMultiDimVector::GetCutValuesFromGlobalAddress(ULong64_t globadd, Float_t *cuts, Int_t &ptbin) const {
130   Int_t ind[fgkMaxNVariables];
131   Bool_t retcode=GetIndicesFromGlobalAddress(globadd,ind,ptbin);
132   if(!retcode) return kFALSE;
133   for(Int_t i=0;i<fNVariables;i++) cuts[i]=GetCutValue(i,ind[i]);
134   return kTRUE;
135 }
136 //______________________________________________________________________
137 ULong64_t AliMultiDimVector::GetGlobalAddressFromIndices(Int_t *ind, Int_t ptbin) const {
138   // Returns the global index of the cell in the matrix
139   Int_t prod=1;
140   ULong64_t elem=0;
141   Int_t indexPlusLevel[fgkMaxNVariables+1];
142   Int_t nOfCellsPlusLevel[fgkMaxNVariables+1];
143   for(Int_t i=0;i<fNVariables;i++){
144     indexPlusLevel[i]=ind[i];
145     nOfCellsPlusLevel[i]=fNCutSteps[i];
146   }
147   indexPlusLevel[fNVariables]=ptbin;
148   nOfCellsPlusLevel[fNVariables]=fNPtBins;
149         
150   for(Int_t i=0;i<fNVariables+1;i++){
151     prod=indexPlusLevel[i];
152     if(i<fNVariables){
153       for(Int_t j=i+1;j<fNVariables+1;j++){
154         prod*=nOfCellsPlusLevel[j];
155       }
156     }
157     elem+=prod;
158   }
159   return elem;
160 }
161 //______________________________________________________________________
162 Bool_t AliMultiDimVector::GetIndicesFromValues(Float_t *values, Int_t *ind) const {
163   // Fills the array of matrix indices strating from variable values
164   for(Int_t i=0;i<fNVariables;i++){
165     if(fGreaterThan[i]){ 
166       if(values[i]<GetMinLimit(i)) return kFALSE;
167       ind[i]=(Int_t)((values[i]-fMinLimits[i])/GetCutStep(i));
168       if(ind[i]>=GetNCutSteps(i)) ind[i]=GetNCutSteps(i)-1;
169     }else{
170       if(values[i]>GetMaxLimit(i)) return kFALSE;
171       ind[i]=(Int_t)((fMaxLimits[i]-values[i])/GetCutStep(i));
172       if(ind[i]>=GetNCutSteps(i)) ind[i]=GetNCutSteps(i)-1;
173     }
174   }
175   return kTRUE;
176 }
177 //______________________________________________________________________
178 ULong64_t AliMultiDimVector::GetGlobalAddressFromValues(Float_t *values, Int_t ptbin) const {
179   // Returns the global index of the cell in the matrix
180    Int_t ind[fgkMaxNVariables];
181    Bool_t retcode=GetIndicesFromValues(values,ind);
182    if(retcode) return GetGlobalAddressFromIndices(ind,ptbin);
183    else{
184      AliError("Values out of range");
185      return fNTotCells+999;
186    }
187 }
188 //_____________________________________________________________________________
189 void AliMultiDimVector::MultiplyBy(Float_t factor){
190   // multiply the AliMultiDimVector by a constant factor
191   for(ULong64_t i=0;i<fNTotCells;i++){
192     if(fVett.At(i)!=-1)
193       fVett.AddAt(fVett.At(i)*factor,i);
194     else fVett.AddAt(-1,i);
195   }
196   
197 }
198 //_____________________________________________________________________________
199 void AliMultiDimVector::Multiply(AliMultiDimVector* mv,Float_t factor){
200   //  Sets AliMultiDimVector=mv*constant factor
201   for(ULong64_t i=0;i<fNTotCells;i++){
202     if(mv->GetElement(i)!=-1)
203       fVett.AddAt(mv->GetElement(i)*factor,i);
204     else fVett.AddAt(-1,i); 
205   }
206 }
207 //_____________________________________________________________________________
208 void AliMultiDimVector::Multiply(AliMultiDimVector* mv1,AliMultiDimVector* mv2){
209   //  Sets AliMultiDimVector=mv1*mv2
210   for(ULong64_t i=0;i<fNTotCells;i++){
211     if(mv1->GetElement(i)!=-1 && mv2->GetElement(i)!=-1)
212       fVett.AddAt(mv1->GetElement(i)*mv2->GetElement(i),i);
213     else fVett.AddAt(-1,i); 
214   }
215 }
216 //_____________________________________________________________________________
217 void AliMultiDimVector::Add(AliMultiDimVector* mv){
218   // Sums contents of mv to AliMultiDimVector
219   if (mv->GetNTotCells()!=fNTotCells){ 
220     AliError("Different dimension of the vectors!!");
221   }else{
222     for(ULong64_t i=0;i<fNTotCells;i++) 
223       if(mv->GetElement(i)!=-1 && fVett.At(i)!=-1)
224         fVett.AddAt(fVett.At(i)+mv->GetElement(i),i);
225       else fVett.AddAt(-1,i); 
226   }
227 }
228 //_____________________________________________________________________________
229 void AliMultiDimVector::Sum(AliMultiDimVector* mv1,AliMultiDimVector* mv2){
230   // Sets AliMultiDimVector=mv1+mv2
231   if (fNTotCells!=mv1->GetNTotCells()&&mv1->GetNTotCells()!=mv2->GetNTotCells()) {
232     AliError("Different dimension of the vectors!!");
233   }
234   else{
235     for(ULong64_t i=0;i<mv1->GetNTotCells();i++) {
236       if(mv1->GetElement(i)!=-1 && mv2->GetElement(i)!=-1)
237         fVett.AddAt(mv1->GetElement(i)+mv2->GetElement(i),i); 
238       else fVett.AddAt(-1,i); 
239     }
240   }
241 }
242 //_____________________________________________________________________________
243 void AliMultiDimVector::LinearComb(AliMultiDimVector* mv1, Float_t norm1, AliMultiDimVector* mv2, Float_t norm2){
244   // Sets AliMultiDimVector=n1*mv1+n2*mv2
245   if (fNTotCells!=mv1->GetNTotCells()&&mv1->GetNTotCells()!=mv2->GetNTotCells()) {
246     AliError("Different dimension of the vectors!!");
247   }
248   else{
249     for(ULong64_t i=0;i<mv1->GetNTotCells();i++) {
250       if(mv1->GetElement(i)!=-1 && mv2->GetElement(i)!=-1)
251         fVett.AddAt(norm1*mv1->GetElement(i)+norm2*mv2->GetElement(i),i); 
252       else fVett.AddAt(-1,i); 
253     }
254   }
255 }
256 //_____________________________________________________________________________
257 void AliMultiDimVector::DivideBy(AliMultiDimVector* mv){
258   // Divide AliMulivector by mv
259   if (mv->GetNTotCells()!=fNTotCells) {
260     AliError("Different dimension of the vectors!!");
261   }
262   else{
263     for(ULong64_t i=0;i<fNTotCells;i++) 
264       if(mv->GetElement(i)!=0 &&mv->GetElement(i)!=-1 && fVett.At(i)!=-1)
265         fVett.AddAt(fVett.At(i)/mv->GetElement(i),i);
266       else fVett.AddAt(-1,i);
267   }
268
269 }
270 //_____________________________________________________________________________
271 void AliMultiDimVector::Divide(AliMultiDimVector* mv1,AliMultiDimVector* mv2){
272   // Sets AliMultiDimVector=mv1/mv2
273   if (fNTotCells!=mv1->GetNTotCells()&&mv1->GetNTotCells()!=mv2->GetNTotCells()) {
274     AliError("Different dimension of the vectors!!");
275   }
276   else{
277     for(ULong64_t i=0;i<mv1->GetNTotCells();i++) 
278       if(mv2->GetElement(i)!=0&& mv2->GetElement(i)!=-1&& mv1->GetElement(i)!=-1)
279         {
280           fVett.AddAt(mv1->GetElement(i)/mv2->GetElement(i),i);
281         }
282       else fVett.AddAt(-1,i);
283   }
284 }
285 //_____________________________________________________________________________
286 void AliMultiDimVector::Sqrt(){
287   // Sqrt of elements of AliMultiDimVector
288   for(ULong64_t i=0;i<fNTotCells;i++) {
289     if(fVett.At(i)>=0) fVett.AddAt(TMath::Sqrt(fVett.At(i)),i);
290     else {
291       fVett.AddAt(-1,i);
292     }
293   }
294 }
295 //_____________________________________________________________________________
296 void AliMultiDimVector::Sqrt(AliMultiDimVector* mv){
297   // Sets AliMultiDimVector=sqrt(mv)
298   for(ULong64_t i=0;i<fNTotCells;i++) 
299     if(mv->GetElement(i)>=0) fVett.AddAt(TMath::Sqrt(mv->GetElement(i)),i);
300     else fVett.AddAt(-1,i);
301 }
302 //_____________________________________________________________________________
303 void AliMultiDimVector::FindMaximum(Float_t& maxValue, Int_t *ind , Int_t ptbin){
304   // finds the element with maximum contents
305   const ULong64_t nelem=fNTotCells/fNPtBins;
306   TArrayF vett;
307   vett.Set(nelem);
308   ULong64_t runningAddress;
309   for(ULong64_t i=0;i<nelem;i++){
310     runningAddress=ptbin+i*fNPtBins;
311     vett.AddAt(fVett[runningAddress],i);
312   }
313   maxValue=TMath::MaxElement(nelem,vett.GetArray());
314   ULong64_t maxAddress=TMath::LocMax(nelem,vett.GetArray());
315   ULong64_t maxGlobalAddress=ptbin+maxAddress*fNPtBins;
316   Int_t checkedptbin;
317   GetIndicesFromGlobalAddress(maxGlobalAddress,ind,checkedptbin);
318 }
319
320 //_____________________________________________________________________________
321 TH2F*  AliMultiDimVector::Project(Int_t firstVar, Int_t secondVar, Int_t* fixedVars, Int_t ptbin, Float_t norm){
322   // Project the AliMultiDimVector on a 2D histogram
323
324   TString hisName=Form("hproj%s%dv%d",GetName(),secondVar,firstVar);
325   TString hisTit=Form("%s vs. %s",fAxisTitles[secondVar].Data(),fAxisTitles[firstVar].Data());
326   TH2F* h2=new TH2F(hisName.Data(),hisTit.Data(),fNCutSteps[firstVar],fMinLimits[firstVar],fMaxLimits[firstVar],fNCutSteps[secondVar],fMinLimits[secondVar],fMaxLimits[secondVar]);
327         
328   Int_t index[fgkMaxNVariables];
329   for(Int_t i=0;i<fNVariables;i++){
330     index[i]=fixedVars[i];
331   }
332   
333   for(Int_t i=0;i<fNCutSteps[firstVar];i++){
334     for(Int_t j=0;j<fNCutSteps[secondVar];j++){
335       index[firstVar]=i;
336       index[secondVar]=j;
337       Float_t cont=GetElement(index,ptbin)/norm;
338       Int_t bin1=i+1;
339       if(!fGreaterThan[firstVar]) bin1=fNCutSteps[firstVar]-i;
340       Int_t bin2=j+1;
341       if(!fGreaterThan[secondVar]) bin2=fNCutSteps[secondVar]-j;
342       h2->SetBinContent(bin1,bin2,cont);
343     }
344   }
345   return h2;
346 }
347 //_____________________________________________________________________________ 
348 void  AliMultiDimVector::GetIntegrationLimits(Int_t iVar, Int_t iCell, Int_t& minbin, Int_t& maxbin) const {
349   // computes bin limits for integrating the AliMultiDimVector
350   minbin=0;
351   maxbin=0;
352   if(iVar<fNVariables){
353     minbin=iCell;
354     maxbin=fNCutSteps[iVar]-1;
355   }
356 }
357 //_____________________________________________________________________________ 
358 void  AliMultiDimVector::GetFillRange(Int_t iVar, Int_t iCell, Int_t& minbin, Int_t& maxbin) const {
359   // computes range of cells passing the cuts for FillAndIntegrate
360   minbin=0;
361   maxbin=0;
362   if(iVar<fNVariables){
363     minbin=0; // bin 0 corresponds to loose cuts
364     maxbin=iCell;
365   }
366 }
367 //_____________________________________________________________________________ 
368 void AliMultiDimVector::Integrate(){
369   // integrates the matrix
370   if(fIsIntegrated){
371     AliError("MultiDimVector already integrated");
372     return;
373   }
374   TArrayF integral(fNTotCells);
375   for(ULong64_t i=0;i<fNTotCells;i++) integral[i]=CountsAboveCell(i);
376   for(ULong64_t i=0;i<fNTotCells;i++) fVett[i]= integral[i];
377   fIsIntegrated=kTRUE;
378 }
379 //_____________________________________________________________________________ 
380 Float_t AliMultiDimVector::CountsAboveCell(ULong64_t globadd) const{
381   // integrates the counts of cells above cell with address globadd
382   Int_t ind[fgkMaxNVariables];
383   Int_t ptbin;
384   GetIndicesFromGlobalAddress(globadd,ind,ptbin);
385   for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
386   Int_t mink[fgkMaxNVariables];
387   Int_t maxk[fgkMaxNVariables];
388   for(Int_t i=0;i<fgkMaxNVariables;i++){
389     GetIntegrationLimits(i,ind[i],mink[i],maxk[i]);
390   }
391   Float_t sumcont=0.;
392   for(Int_t k0=mink[0]; k0<=maxk[0]; k0++){
393     for(Int_t k1=mink[1]; k1<=maxk[1]; k1++){
394       for(Int_t k2=mink[2]; k2<=maxk[2]; k2++){
395         for(Int_t k3=mink[3]; k3<=maxk[3]; k3++){
396           for(Int_t k4=mink[4]; k4<=maxk[4]; k4++){
397             for(Int_t k5=mink[5]; k5<=maxk[5]; k5++){
398               for(Int_t k6=mink[6]; k6<=maxk[6]; k6++){
399                 for(Int_t k7=mink[7]; k7<=maxk[7]; k7++){
400                   for(Int_t k8=mink[8]; k8<=maxk[8]; k8++){
401                     for(Int_t k9=mink[9]; k9<=maxk[9]; k9++){
402                       Int_t currentBin[fgkMaxNVariables]={k0,k1,k2,k3,k4,k5,k6,k7,k8,k9};
403                       sumcont+=GetElement(currentBin,ptbin);
404                     }
405                   }
406                 }
407               }
408             }
409           }
410         }
411       }
412     }
413   }
414   return sumcont;
415 }
416 //_____________________________________________________________________________ 
417 void AliMultiDimVector::Fill(Float_t* values, Int_t ptbin){
418   // fills the cells of AliMultiDimVector corresponding to values
419   if(fIsIntegrated){
420     AliError("MultiDimVector already integrated -- Use FillAndIntegrate");
421     return;
422   }
423   Int_t ind[fgkMaxNVariables];
424   Bool_t retcode=GetIndicesFromValues(values,ind);
425   for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
426   if(retcode) IncrementElement(ind,ptbin);
427 }
428 //_____________________________________________________________________________ 
429 void AliMultiDimVector::FillAndIntegrate(Float_t* values, Int_t ptbin){
430   // fills the cells of AliMultiDimVector passing the cuts
431   // The number of nested loops must match fgkMaxNVariables!!!!
432   fIsIntegrated=kTRUE;
433   Int_t ind[fgkMaxNVariables];
434   Bool_t retcode=GetIndicesFromValues(values,ind);
435   if(!retcode) return;
436   for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
437   Int_t mink[fgkMaxNVariables];
438   Int_t maxk[fgkMaxNVariables];
439   for(Int_t i=0;i<fgkMaxNVariables;i++){
440     GetFillRange(i,ind[i],mink[i],maxk[i]);
441   }
442   for(Int_t k0=mink[0]; k0<=maxk[0]; k0++){
443     for(Int_t k1=mink[1]; k1<=maxk[1]; k1++){
444       for(Int_t k2=mink[2]; k2<=maxk[2]; k2++){
445         for(Int_t k3=mink[3]; k3<=maxk[3]; k3++){
446           for(Int_t k4=mink[4]; k4<=maxk[4]; k4++){
447             for(Int_t k5=mink[5]; k5<=maxk[5]; k5++){
448               for(Int_t k6=mink[6]; k6<=maxk[6]; k6++){
449                 for(Int_t k7=mink[7]; k7<=maxk[7]; k7++){
450                   for(Int_t k8=mink[8]; k8<=maxk[8]; k8++){
451                     for(Int_t k9=mink[9]; k9<=maxk[9]; k9++){
452                       Int_t currentBin[fgkMaxNVariables]={k0,k1,k2,k3,k4,k5,k6,k7,k8,k9};
453                       IncrementElement(currentBin,ptbin);
454                     }
455                   }
456                 }
457               }
458             }
459           }
460         }
461       }
462     }
463   }
464
465 }
466 //_____________________________________________________________________________ 
467 void AliMultiDimVector::SuppressZeroBKGEffect(AliMultiDimVector* mvBKG){
468   // Sets to zero elements for which mvBKG=0
469   for(ULong64_t i=0;i<fNTotCells;i++)
470     if(mvBKG->GetElement(i)==0) fVett.AddAt(0,i);
471 }
472 //_____________________________________________________________________________ 
473 AliMultiDimVector* AliMultiDimVector:: ShrinkPtBins(Int_t firstBin, Int_t lastBin){
474   // sums the elements of pt bins between firstBin and lastBin
475   if(firstBin<0 || lastBin>=fNPtBins || firstBin>=lastBin){
476     AliError("Bad numbers of Pt bins to be shrinked");
477     return 0;
478   }
479   Int_t nofcells[fgkMaxNVariables];
480   Float_t loosecuts[fgkMaxNVariables];
481   Float_t tightcuts[fgkMaxNVariables];
482   TString axisTitles[fgkMaxNVariables];
483   for(Int_t i=0;i<fNVariables;i++){
484     nofcells[i]=fNCutSteps[i];
485     if(fGreaterThan[i]){
486       loosecuts[i]=fMinLimits[i];
487       tightcuts[i]=fMaxLimits[i];
488     }else{
489       loosecuts[i]=fMaxLimits[i];
490       tightcuts[i]=fMinLimits[i];
491     }
492     axisTitles[i]=fAxisTitles[i];
493   }
494   Int_t newNptbins=fNPtBins-(lastBin-firstBin);
495   AliMultiDimVector* shrinkedMV=new AliMultiDimVector(GetName(),GetTitle(),fNVariables,newNptbins,nofcells,loosecuts,tightcuts,axisTitles);
496   
497   ULong64_t nOfPointsPerPtbin=fNTotCells/fNPtBins;
498   ULong64_t addressOld,addressNew;
499   Int_t npb,opb;
500   for(npb=0;npb<firstBin;npb++){
501     opb=npb;
502     for(ULong64_t k=0;k<nOfPointsPerPtbin;k++){
503       addressOld=opb+k*fNPtBins;
504       addressNew=npb+k*newNptbins;
505       shrinkedMV->SetElement(addressNew,fVett[addressOld]);
506     }
507   }
508   npb=firstBin;
509   for(ULong64_t k=0;k<nOfPointsPerPtbin;k++){
510     Float_t summedValue=0.;
511     for(opb=firstBin;opb<=lastBin;opb++){
512       addressOld=opb+k*fNPtBins;
513       summedValue+=fVett[addressOld];
514     }
515     addressNew=npb+k*newNptbins;
516     shrinkedMV->SetElement(addressNew,summedValue);
517   }
518   for(npb=firstBin+1;npb<newNptbins;npb++){
519     opb=npb+(lastBin-firstBin);
520     for(ULong64_t k=0;k<nOfPointsPerPtbin;k++){
521       addressOld=opb+k*fNPtBins;
522       addressNew=npb+k*newNptbins;
523       shrinkedMV->SetElement(addressNew,fVett[addressOld]);
524     }
525   }
526   return shrinkedMV;
527 }