Install macros
[u/mrichter/AliRoot.git] / TPC / TPCcalib / AliTPCcalibTracks.cxx
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 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
19 //     Class for calibration of the cluster properties:                       //
20 //     Cluster position resolution (RMS)  and short term distortions (within pad or within time bin)
21 // 
22 //  Algorithm:
23 //    1. Creation of the residual/properties  histograms
24 //    2. Filling of the histograms.
25 //       2.a Traklet fitting
26 //       2.b Resudual filling
27 //    3. Postprocessing: Creation of the tree with the mean values and RMS at different bins
28 //    4.               : Paramaterization fitting. 
29 //                       In old AliRoot the RMS paramterization was used to characterize cluster resolution
30 //                       Mean value /bias was neglected
31 //    5. To be implemented 
32 //          a.) lookup table for the distortion parmaterization - THn
33 //          b.) Usage in the reconstruction  
34 //                               
35 // 
36 //   1. Creation of the histograms:   MakeHistos() 
37 //
38 //         6 n dimensional histograms  are filled during the calibration:
39 //         cluster  performance bins
40 //         0 - variable of interest
41 //         1 - pad type   - 0- IROC Short 1- OCROC medium 2 - OROC long pads
42 //         2 - drift length - drift length -0-1
43 //         3 - Qmax         - Qmax  - 2- 400
44 //         4 - cog          - COG position - 0-1
45 //         5 - tan(phi)     - local  angle
46 //        Histograms:
47 //        THnSparse  *fHisDeltaY;    // THnSparse - delta Y between the cluster and tracklet interpolation (+-X(5?) rows) 
48 //        THnSparse  *fHisDeltaZ;    // THnSparse - delta Z 
49 //        THnSparse  *fHisRMSY;      // THnSparse - rms Y 
50 //        THnSparse  *fHisRMSZ;      // THnSparse - rms Z 
51 //        THnSparse  *fHisQmax;      // THnSparse - qmax 
52 //        THnSparse  *fHisQtot;      // THnSparse - qtot 
53 //
54 //
55 //
56
57
58                      //
59 //     The parameter 'clusterParam', a AliTPCClusterParam object             //
60 //      (needed for TPC cluster error and shape parameterization)            //
61 //
62 //     Normally you get this object out of the file 'TPCClusterParam.root'   //
63 //     In the parameter 'cuts' the cuts are specified, that decide           //
64 //     weather a track will be accepted for calibration or not.              //
65 //                                                                           //
66 //       
67 //     The data flow:
68 //     
69 /*
70    Raw Data -> Local Reconstruction -> Tracking ->     Calibration -> RefData (component itself)
71                Offline/HLT             Offline/HLT                    OCDB entries (AliTPCClusterParam) 
72 */            
73
74 /*
75   Example usage - creation of the residual trees:
76   TFile f("CalibObjects.root");
77   AliTPCcalibTracks *calibTracks = ( AliTPCcalibTracks *)TPCCalib->FindObject("calibTracks");
78   TH2 * his2 = calibTracks->fHisDeltaY->Projection(0,4);
79   his2->SetName("his2");
80   his2->FitSlicesY();
81
82   
83   TTreeSRedirector *pcstream = new TTreeSRedirector("clusterLookup.root");
84   AliTPCcalibTracks::MakeSummaryTree(calibTracks->fHisDeltaY, pcstream, 0);
85   AliTPCcalibTracks::MakeSummaryTree(calibTracks->fHisDeltaZ, pcstream, 1);
86   delete pcstream;
87   TFile fl("clusterLookup.root");
88   TCut cutNI="cogA==0&&angleA==0&&qmaxA==0&&zA==0&&ipadA==0"
89
90 */
91
92                                                                //
93 ///////////////////////////////////////////////////////////////////////////////
94
95 //
96 // ROOT includes 
97 //
98 #include <iostream>
99 #include <fstream>
100 using namespace std;
101
102 #include <TROOT.h>
103 #include <TChain.h>
104 #include <TFile.h>
105 #include <TH3F.h>
106 #include <TProfile.h>
107
108 //
109 //#include <TPDGCode.h>
110 #include <TStyle.h>
111 #include "TLinearFitter.h"
112 //#include "TMatrixD.h"
113 #include "TTreeStream.h"
114 #include "TF1.h"
115 #include <TCanvas.h>
116 #include <TGraph2DErrors.h>
117 #include "TPostScript.h"
118 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
119 #include "TCint.h"
120 #endif
121
122 #include <TH2D.h>
123 #include <TF2.h>
124 #include <TSystem.h>
125 #include <TCollection.h>
126 #include <iostream>
127 #include <TLinearFitter.h>
128 #include <TString.h>
129
130 //
131 // AliROOT includes 
132 //
133 #include "AliMagF.h"
134 #include "AliTracker.h"
135 #include "AliESD.h"
136 #include "AliESDtrack.h"
137 #include "AliESDfriend.h"
138 #include "AliESDfriendTrack.h" 
139 #include "AliTPCseed.h"
140 #include "AliTPCclusterMI.h"
141 #include "AliTPCROC.h"
142
143 #include "AliTPCParamSR.h"
144 #include "AliTrackPointArray.h"
145 #include "AliTPCcalibTracks.h"
146 #include "AliTPCClusterParam.h"
147 #include "AliTPCcalibTracksCuts.h"
148 #include "AliTPCCalPad.h"
149 #include "AliTPCCalROC.h"
150 #include "TText.h"
151 #include "TPaveText.h"
152 #include "TSystem.h"
153 #include "TStatToolkit.h"
154 #include "TCut.h"
155 #include "THnSparse.h"
156 #include "AliRieman.h"
157 #include "TRandom.h"
158
159
160 Double_t          AliTPCcalibTracks::fgkMergeEntriesCut=10000000.; //10**7 clusters
161
162 ClassImp(AliTPCcalibTracks)
163
164
165 AliTPCcalibTracks::AliTPCcalibTracks():
166   AliTPCcalibBase(),
167   fClusterParam(0),
168   fROC(0),
169   fHisDeltaY(0),    // THnSparse - delta Y 
170   fHisDeltaZ(0),    // THnSparse - delta Z 
171   fHisRMSY(0),      // THnSparse - rms Y 
172   fHisRMSZ(0),      // THnSparse - rms Z 
173   fHisQmax(0),      // THnSparse - qmax 
174   fHisQtot(0),      // THnSparse - qtot 
175   fPtDownscaleRatio(20),       // pt downscaling ratio (use subsample of data)
176   fQDownscaleRatio(20),        // Q downscaling ratio (use subsample of dta)
177   fArrayQDY(0), 
178   fArrayQDZ(0), 
179   fArrayQRMSY(0),
180   fArrayQRMSZ(0),
181   fResolY(0),
182   fResolZ(0),
183   fRMSY(0),
184   fRMSZ(0),
185   fCuts(0),
186   fRejectedTracksHisto(0),
187   fClusterCutHisto(0),
188   fCalPadClusterPerPad(0),
189   fCalPadClusterPerPadRaw(0)
190
191    // 
192    // AliTPCcalibTracks default constructor
193    //    
194   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl;  
195 }   
196
197
198
199 AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
200   AliTPCcalibBase(calibTracks),
201   fClusterParam(0),
202   fROC(0),
203   fHisDeltaY(0),    // THnSparse - delta Y 
204   fHisDeltaZ(0),    // THnSparse - delta Z 
205   fHisRMSY(0),      // THnSparse - rms Y 
206   fHisRMSZ(0),      // THnSparse - rms Z 
207   fHisQmax(0),      // THnSparse - qmax 
208   fHisQtot(0),      // THnSparse - qtot 
209   fPtDownscaleRatio(20),       // pt downscaling ratio (use subsample of data)
210   fQDownscaleRatio(20),        // Q downscaling ratio (use subsample of dta)
211   fArrayQDY(0), 
212   fArrayQDZ(0), 
213   fArrayQRMSY(0),
214   fArrayQRMSZ(0),
215   fResolY(0),
216   fResolZ(0),
217   fRMSY(0),
218   fRMSZ(0),
219   fCuts(0),
220   fRejectedTracksHisto(0),
221   fClusterCutHisto(0),
222   fCalPadClusterPerPad(0),
223   fCalPadClusterPerPadRaw(0)
224 {
225    // 
226    // AliTPCcalibTracks copy constructor
227    // 
228   if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl;
229    
230    Bool_t dirStatus = TH1::AddDirectoryStatus();
231    TH1::AddDirectory(kFALSE);
232    
233    Int_t length = -1;
234    
235    (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
236    fArrayQDY= new TObjArray(length);
237    fArrayQDZ= new TObjArray(length);
238    fArrayQRMSY= new TObjArray(length);
239    fArrayQRMSZ= new TObjArray(length);
240    for (Int_t i = 0; i < length; i++) {
241       fArrayQDY->AddAt( ((TH1F*)calibTracks.fArrayQDY->At(i)->Clone()), i);
242       fArrayQDZ->AddAt( ((TH1F*)calibTracks.fArrayQDZ->At(i)->Clone()), i);
243       fArrayQRMSY->AddAt( ((TH1F*)calibTracks.fArrayQRMSY->At(i)->Clone()), i);
244       fArrayQRMSZ->AddAt( ((TH1F*)calibTracks.fArrayQRMSZ->At(i)->Clone()), i);
245    }
246    
247    (calibTracks.fResolY) ? length = calibTracks.fResolY->GetEntriesFast() : length = -1;
248    fResolY = new TObjArray(length);
249    fResolZ = new TObjArray(length);
250    fRMSY = new TObjArray(length);
251    fRMSZ = new TObjArray(length);
252    for (Int_t i = 0; i < length; i++) {
253       fResolY->AddAt( ((TH1F*)calibTracks.fResolY->At(i)->Clone()), i);
254       fResolZ->AddAt( ((TH1F*)calibTracks.fResolZ->At(i)->Clone()), i);
255       fRMSY->AddAt( ((TH1F*)calibTracks.fRMSY->At(i)->Clone()), i);
256       fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
257    } 
258    
259    
260    fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
261    fRejectedTracksHisto    = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
262    fCalPadClusterPerPad    = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
263    fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
264
265    fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(), 
266       calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise());
267    SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle());
268    TH1::AddDirectory(dirStatus); // set status back to original status
269 //    cout << "+++++ end of copy constructor +++++" << endl;   // TO BE REMOVED
270 }
271
272
273 AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibTracks){
274   //
275   // assgnment operator
276   //
277   if (this != &calibTracks) {
278     new (this) AliTPCcalibTracks(calibTracks);
279   }
280   return *this;
281
282 }
283
284
285 AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam,  AliTPCcalibTracksCuts* cuts, Int_t logLevel) : 
286   AliTPCcalibBase(),
287   fClusterParam(0),
288   fROC(0),
289   fHisDeltaY(0),    // THnSparse - delta Y 
290   fHisDeltaZ(0),    // THnSparse - delta Z 
291   fHisRMSY(0),      // THnSparse - rms Y 
292   fHisRMSZ(0),      // THnSparse - rms Z 
293   fHisQmax(0),      // THnSparse - qmax 
294   fHisQtot(0),      // THnSparse - qtot 
295   fPtDownscaleRatio(20),       // pt downscaling ratio (use subsample of data)
296   fQDownscaleRatio(20),        // Q downscaling ratio (use subsample of dta)
297   fArrayQDY(0), 
298   fArrayQDZ(0), 
299   fArrayQRMSY(0),
300   fArrayQRMSZ(0),
301   fResolY(0),
302   fResolZ(0),
303   fRMSY(0),
304   fRMSZ(0),
305   fCuts(0),
306   fRejectedTracksHisto(0),
307   fClusterCutHisto(0),
308   fCalPadClusterPerPad(0),
309   fCalPadClusterPerPadRaw(0)
310  {
311    // 
312    // AliTPCcalibTracks constructor
313    // specify 'name' and 'title' of your object
314    // specify 'clusterParam', (needed for TPC cluster error and shape parameterization)
315    // In the parameter 'cuts' the cuts are specified, that decide           
316    // weather a track will be accepted for calibration or not.              
317    //
318    // fDebugLevel - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStreamer, 6: waste your screen
319    // 
320    // All histograms are instatiated in this constructor.
321    // 
322    this->SetName(name);
323    this->SetTitle(title);
324
325    if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl;
326 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
327    G__SetCatchException(0);     
328 #endif
329    
330    fClusterParam = clusterParam;
331    if (fClusterParam){
332      fClusterParam->SetInstance(fClusterParam);
333    }
334    else {
335      Error("AliTPCcalibTracks","No cluster parametrization found! A valid clusterParam object is needed in the constructor. (To be found in 'TPCClusterParam.root'.)");
336    } 
337    fCuts = cuts;
338    SetDebugLevel(logLevel);
339    MakeHistos();
340    
341    TH1::AddDirectory(kFALSE);
342    
343    fRejectedTracksHisto    = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
344    fCalPadClusterPerPad    = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
345    fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
346    fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
347    
348    
349    TH1::AddDirectory(kFALSE);
350    
351    
352    fResolY = new TObjArray(3);
353    fResolZ = new TObjArray(3);
354    fRMSY   = new TObjArray(3);
355    fRMSZ   = new TObjArray(3);
356    TH3F * his3D;
357    //
358    his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
359    fResolY->AddAt(his3D,0);     
360    his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
361    fResolY->AddAt(his3D,1);
362    his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
363    fResolY->AddAt(his3D,2);
364    //
365    his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
366    fResolZ->AddAt(his3D,0);
367    his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
368    fResolZ->AddAt(his3D,1);
369    his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
370    fResolZ->AddAt(his3D,2);
371    //
372    his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
373    fRMSY->AddAt(his3D,0);
374    his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
375    fRMSY->AddAt(his3D,1);
376    his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
377    fRMSY->AddAt(his3D,2);
378    //
379    his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
380    fRMSZ->AddAt(his3D,0);
381    his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
382    fRMSZ->AddAt(his3D,1);
383    his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
384    fRMSZ->AddAt(his3D,2);
385    //
386       
387    TH1::AddDirectory(kFALSE);
388    
389    fArrayQDY = new TObjArray(300);
390    fArrayQDZ = new TObjArray(300);
391    fArrayQRMSY = new TObjArray(300);
392    fArrayQRMSZ = new TObjArray(300);
393    for (Int_t iq = 0; iq <= 10; iq++){
394       for (Int_t ipad = 0; ipad < 3; ipad++){
395          Int_t   bin   = GetBin(iq, ipad);
396          Float_t qmean = GetQ(bin);
397          char hname[200];
398          snprintf(hname,100,"ResolY Pad%d Qmiddle%f",ipad, qmean);
399          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
400          fArrayQDY->AddAt(his3D, bin);
401          snprintf(hname,100,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
402          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
403          fArrayQDZ->AddAt(his3D, bin);
404          snprintf(hname,100,"RMSY Pad%d Qmiddle%f",ipad, qmean);
405          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
406          fArrayQRMSY->AddAt(his3D, bin);
407          snprintf(hname,100,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
408          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
409          fArrayQRMSZ->AddAt(his3D, bin);
410       }
411    }
412    
413    
414
415    if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl; 
416    cout << "end of main constructor" << endl; // TO BE REMOVED
417 }    
418
419
420 AliTPCcalibTracks::~AliTPCcalibTracks() {
421    // 
422    // AliTPCcalibTracks destructor
423    // 
424    
425   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
426    Int_t length = 0;
427    
428    
429    if (fResolY) {
430      length = fResolY->GetEntriesFast();
431      for (Int_t i = 0; i < length; i++){
432        delete fResolY->At(i);
433        delete fResolZ->At(i);
434        delete fRMSY->At(i);
435        delete fRMSZ->At(i);
436      }
437      delete fResolY;
438      delete fResolZ;
439      delete fRMSY;
440      delete fRMSZ;
441    }
442    
443    if (fArrayQDY) {
444      length = fArrayQDY->GetEntriesFast();
445      for (Int_t i = 0; i < length; i++){
446        delete fArrayQDY->At(i);
447        delete fArrayQDZ->At(i);
448        delete fArrayQRMSY->At(i);
449        delete fArrayQRMSZ->At(i);
450      }
451    }
452    
453    
454     
455    delete fArrayQDY;
456    delete fArrayQDZ;
457    delete fArrayQRMSY;
458    delete fArrayQRMSZ;
459    
460   delete fRejectedTracksHisto;
461   delete fClusterCutHisto;
462   if (fCalPadClusterPerPad)    delete fCalPadClusterPerPad;
463   if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
464   delete fHisDeltaY;    // THnSparse - delta Y 
465   delete fHisDeltaZ;    // THnSparse - delta Z 
466   delete fHisRMSY;      // THnSparse - rms Y 
467   delete fHisRMSZ;      // THnSparse - rms Z 
468   delete fHisQmax;      // THnSparse - qmax 
469   delete fHisQtot;      // THnSparse - qtot 
470
471 }
472    
473   
474
475 void AliTPCcalibTracks::Process(AliTPCseed *track){
476    // 
477    // To be called in the selector
478    // first AcceptTrack is evaluated, then calls all the following analyse functions: 
479    // FillResolutionHistoLocal(track)
480
481    // 
482   Double_t scalept= TMath::Min(1./TMath::Abs(track->GetParameter()[4]),2.)/0.5;
483   Bool_t   isSelected = (TMath::Exp(scalept)>fPtDownscaleRatio*gRandom->Rndm());
484   if (!isSelected) return;
485
486   if (GetDebugLevel() > 5) Info("Process","Starting to process the track...");
487    Int_t accpetStatus = AcceptTrack(track);
488    if (accpetStatus == 0) {
489       FillResolutionHistoLocal(track);
490    }
491    else fRejectedTracksHisto->Fill(accpetStatus);
492 }
493
494
495
496 Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){
497   //
498   // calculate bins for given q and pad type 
499   // used in TObjArray
500   //
501   Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 );  
502   res *= 3;
503   res += pad;
504   return res;
505 }
506
507
508 Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){
509   //
510   // calculate bins for given iq and pad type 
511   // used in TObjArray
512   //
513   return iq * 3 + pad;;
514 }
515
516
517 Float_t AliTPCcalibTracks::GetQ(Int_t bin){
518    // 
519    // returns to bin belonging charge
520    // (bin / 3 + 3)^2
521    // 
522    Int_t bin0 = bin / 3;
523    bin0 += 3;
524    return bin0 * bin0;
525 }
526
527
528 Float_t AliTPCcalibTracks::GetPad(Int_t bin){
529    // 
530    // returns to bin belonging pad
531    // bin % 3
532    // 
533    return bin % 3; 
534 }
535
536
537
538 Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
539   //
540   // Function, that decides wheather a given track is accepted for 
541   // the analysis or not. 
542   // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts'
543   // Returns 0 if a track is accepted or an integer different from 0 
544   // to indicate the failed cut
545   //
546   const Int_t   kMinClusters  = fCuts->GetMinClusters();
547   const Float_t kMinRatio     = fCuts->GetMinRatio();
548   const Float_t kMax1pt       = fCuts->GetMax1pt();
549   const Float_t kEdgeYXCutNoise    = fCuts->GetEdgeYXCutNoise();
550   const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise();
551   
552   //
553   // edge induced noise tracks - NEXT RELEASE will be removed during tracking
554   if ( TMath::Abs(track->GetY() / track->GetX()) > kEdgeYXCutNoise )
555     if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return 1;
556   if (track->GetNumberOfClusters() < kMinClusters) return 2;
557   Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.);
558   if (ratio < kMinRatio) return 3;
559 //   Float_t mpt = track->Get1Pt();       // Get1Pt() doesn't exist any more
560   Float_t mpt = track->GetSigned1Pt();
561   if (TMath::Abs(mpt) > kMax1pt) return 4;
562   //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
563   //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
564   //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
565   
566   if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");  
567   return 0;
568 }
569
570
571 void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
572    //
573    // fill resolution histograms - localy - tracklet in the neighborhood
574    // write debug information to 'TPCSelectorDebug.root'
575    // 
576    // _ the main function, called during track analysis _
577    // 
578    // loop over all padrows along the track
579    // fit tracklets (length: 13 clusters) calculate mean chi^2 for this track-fit in Y and Z direction
580    // 
581    // loop again over all padrows along the track
582    // fit tracklet (clusters in given padrow +- kDelta padrows) 
583    // with polynom of 2nd order and two polynoms of 1st order
584    // take both polynoms of 1st order, calculate difference of their parameters
585    // add covariance matrixes and calculate chi2 of this difference
586    // if this chi2 is bigger than a given threshold, assume that the current cluster is
587    // a kink an goto next padrow
588    // if not:
589    // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
590    // 
591    // write debug information to 'TPCSelectorDebug.root'
592    // only for every kDeltaWriteDebugStream'th padrow to reduce data volume 
593    // and to avoid redundant data
594    // 
595
596   /*  {//SG
597     static long n=0;
598     if( n==0 ){
599       if( fClusterParam && fClusterParam->GetWaveCorrectionMap() ){
600         cout<<"Map found "<<endl; 
601       }else{
602         cout<<"Can't find the map "<<endl;    
603       }
604     }
605     if( n%100 ==0 )cout<<n<<" Tracks processed"<<endl;
606     n++;
607   }*/
608   static TLinearFitter fFitterParY(3,"pol2");    // 
609   static TLinearFitter fFitterParZ(3,"pol2");    //
610   static TLinearFitter fFitterParYWeight(3,"pol2");    // 
611   static TLinearFitter fFitterParZWeight(3,"pol2");    //
612   fFitterParY.StoreData(kTRUE);
613   fFitterParZ.StoreData(kTRUE);
614   fFitterParYWeight.StoreData(kTRUE);
615   fFitterParZWeight.StoreData(kTRUE);
616   if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
617   const Int_t   kDelta     = 10;          // delta rows to fit
618   const Float_t kMinRatio  = 0.75;        // minimal ratio
619   const Float_t kChi2Cut   =  10.;          // cut chi2 - left right
620   const Float_t kSigmaCut  = 3.;        //sigma cut
621   const Float_t kdEdxCut=300;
622   const Float_t kPtCut=0.300;
623
624   if (track->GetTPCsignal()>kdEdxCut) return;  
625   if (TMath::Abs(AliTracker::GetBz())>0.1 &&TMath::Abs(track->Pt())<kPtCut) return;  
626
627   // estimate mean error
628   Int_t nTrackletsAll = 0;       // number of tracklets for given track
629   Float_t csigmaY     = 0;       // mean sigma for tracklet refit in Y direction
630   Float_t csigmaZ     = 0;       // mean sigma for tracklet refit in Z direction
631   Int_t nClusters     = 0;       // working variable, number of clusters per tracklet
632   Int_t sectorG       = -1;      // working variable, sector of tracklet, has to stay constant for one tracklet
633   Double_t refX=0;
634   // ---------------------------------------------------------------------
635   for (Int_t irow = 0; irow < 159; irow++){
636     // loop over all rows along the track
637     // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
638     // calculate mean chi^2 for this track-fit in Y and Z direction
639     AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
640     if (!cluster0) continue;  // no cluster found
641     Int_t sector = cluster0->GetDetector();
642     
643     Int_t ipad = TMath::Nint(cluster0->GetPad());
644     Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
645     fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
646     
647     if (sector != sectorG){
648       // track leaves sector before it crossed enough rows to fit / initialization
649       nClusters = 0;
650       fFitterParY.ClearPoints();
651       fFitterParZ.ClearPoints();
652       sectorG = sector;
653     }
654     else {
655       nClusters++;
656       if (refX<1) refX=cluster0->GetX()+kDelta*0.5;
657       Double_t x = cluster0->GetX()-refX;
658       fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
659       fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
660       //
661       if ( nClusters >= kDelta + 3 ){  
662         // if more than 13 (kDelta+3) clusters were added to the fitters
663         // fit the tracklet, increase trackletCounter
664         fFitterParY.Eval();
665         fFitterParZ.Eval();
666         nTrackletsAll++;
667         csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
668         csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
669         nClusters = -1;
670         fFitterParY.ClearPoints();
671         fFitterParZ.ClearPoints();
672         refX=0;
673       }
674     }
675   }      // for (Int_t irow = 0; irow < 159; irow++)
676   // mean chi^2 for all tracklet fits in Y and in Z direction: 
677   csigmaY = TMath::Sqrt(TMath::Abs(csigmaY) / (nTrackletsAll+0.1));
678   csigmaZ = TMath::Sqrt(TMath::Abs(csigmaZ) / (nTrackletsAll+0.1));
679   if (csigmaY<=TMath::KUncertainty() || csigmaZ<= TMath::KUncertainty()) return;
680   // ---------------------------------------------------------------------
681   //
682   //
683
684   for (Int_t irow = kDelta; irow < 159-kDelta; irow++){
685     // loop again over all rows along the track
686     // do analysis
687     // 
688     Int_t nclFound = 0;  // number of clusters in the neighborhood
689     Int_t ncl0 = 0;      // number of clusters in rows < rowOfCenterCluster
690     Int_t ncl1 = 0;      // number of clusters in rows > rowOfCenterCluster
691     AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
692     if (!cluster0) continue;
693     Int_t sector = cluster0->GetDetector();
694     Float_t xref = cluster0->GetX();
695     
696     // Make Fit
697     fFitterParY.ClearPoints();
698     fFitterParZ.ClearPoints();    
699     fFitterParYWeight.ClearPoints();
700     fFitterParZWeight.ClearPoints();    
701     // fit tracklet (clusters in given padrow +- kDelta padrows) 
702     // with polynom of 2nd order and two polynoms of 1st order
703     // take both polynoms of 1st order, calculate difference of their parameters
704     // add covariance matrixes and calculate chi2 of this difference
705     // if this chi2 is bigger than a given threshold, assume that the current cluster is
706     // a kink an goto next padrow    
707
708     // first fit the track angle for wave correction
709
710     AliRieman riemanFitAngle( 2*kDelta ); //SG
711
712     if( fClusterParam && fClusterParam->GetWaveCorrectionMap() ){
713       for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
714         // loop over irow +- kDelta rows (neighboured rows)
715         if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
716         AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
717         if (!currentCluster) continue;
718         if (currentCluster->GetType() < 0) continue;
719         if (currentCluster->GetDetector() != sector) continue;      
720         riemanFitAngle.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY,csigmaZ);
721       }
722       if( riemanFitAngle.GetN()>3 ) riemanFitAngle.Update();    
723     }
724
725     // do fit
726
727     AliRieman riemanFit(2*kDelta);
728     AliRieman riemanFitW(2*kDelta);
729     for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
730       // loop over irow +- kDelta rows (neighboured rows)
731       // 
732       // 
733       if (idelta == 0) continue;                                // don't use center cluster
734       if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
735       AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
736       if (!currentCluster) continue;
737       if (currentCluster->GetType() < 0) continue;
738       if (currentCluster->GetDetector() != sector) continue;
739       nclFound++;
740       if (idelta < 0){
741         ncl0++;
742       }
743       if (idelta > 0){
744         ncl1++;
745       }
746       //SG!!!
747       double dY = 0;
748       if( fClusterParam ){
749         Int_t padSize = 0;                          // short pads
750         if (currentCluster->GetDetector() >= 36) {
751           padSize = 1;                              // medium pads 
752           if (currentCluster->GetRow() > 63) padSize = 2; // long pads
753         }
754         dY = fClusterParam->GetWaveCorrection( padSize, 
755                                                currentCluster->GetZ(), 
756                                                currentCluster->GetMax(),
757                                                currentCluster->GetPad(),
758                                                riemanFitAngle.GetDYat(currentCluster->GetX())
759                                                );       
760       }
761       riemanFit.AddPoint(currentCluster->GetX(), currentCluster->GetY()-dY,currentCluster->GetZ(), csigmaY,csigmaZ);
762       riemanFitW.AddPoint(currentCluster->GetX(), currentCluster->GetY()-dY,currentCluster->GetZ(), TMath::Abs(csigmaY*TMath::Sqrt(1+TMath::Abs(idelta))),TMath::Abs(csigmaZ*TMath::Sqrt(1+TMath::Abs(idelta))));
763     }  // loop over neighbourhood for fitter filling 
764     if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
765     if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
766     if (nclFound < kDelta * kMinRatio) continue;    // if not enough clusters (7.5) found in neighbourhood goto next padrow
767     riemanFit.Update();
768     riemanFitW.Update();
769     Double_t chi2R=TMath::Sqrt(riemanFit.CalcChi2()/(2*nclFound-5));
770     Double_t chi2RW=TMath::Sqrt(riemanFitW.CalcChi2()/(2*nclFound-5));
771     Double_t paramYR[3], paramZR[3];
772     Double_t paramYRW[3], paramZRW[3];
773     //
774     paramYR[0]    = riemanFit.GetYat(xref);
775     paramZR[0]    = riemanFit.GetZat(xref);
776     paramYRW[0]   = riemanFitW.GetYat(xref);
777     paramZRW[0]   = riemanFitW.GetZat(xref);
778     //
779     paramYR[1]    = riemanFit.GetDYat(xref);
780     paramZR[1]    = riemanFit.GetDZat(xref);
781     paramYRW[1]   = riemanFitW.GetDYat(xref);
782     paramZRW[1]   = riemanFitW.GetDZat(xref);
783     //
784     Int_t reject=0;
785     if (chi2R>kChi2Cut) reject+=1;
786     if (chi2RW>kChi2Cut) reject+=2;
787     if (TMath::Abs(paramYR[0]-paramYRW[0])>kSigmaCut*csigmaY) reject+=4;
788     if (TMath::Abs(paramZR[0]-paramZRW[0])>kSigmaCut*csigmaZ) reject+=8;
789     if (TMath::Abs(paramYR[1]-paramYRW[1])>csigmaY) reject+=16;
790     if (TMath::Abs(paramZR[1]-paramZRW[1])>csigmaZ) reject+=32;
791     //
792     TTreeSRedirector *cstream = GetDebugStreamer();    
793     // get fit parameters from pol2 fit:     
794     Double_t tracky = paramYR[0];
795     Double_t trackz = paramZR[0];
796     Float_t  deltay = cluster0->GetY()-tracky;
797     Float_t  deltaz = cluster0->GetZ()-trackz;
798     Float_t  angley = paramYR[1];
799     Float_t  anglez = paramZR[1];
800     
801     Int_t padSize = 0;                          // short pads
802     if (cluster0->GetDetector() >= 36) {
803       padSize = 1;                              // medium pads 
804       if (cluster0->GetRow() > 63) padSize = 2; // long pads
805     }
806     Int_t ipad = TMath::Nint(cluster0->GetPad());
807     if (cstream){
808       Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
809       (*cstream)<<"Resol0"<<    
810         "run="<<fRun<<              //  run number
811         "event="<<fEvent<<          //  event number
812         "time="<<fTime<<            //  time stamp of event
813         "trigger="<<fTrigger<<      //  trigger
814         "mag="<<fMagF<<             //  magnetic field        
815         "padSize="<<padSize<<
816         //
817         "reject="<<reject<<
818         "cl.="<<cluster0<<          // cluster info
819         "xref="<<xref<<             // cluster reference X
820         //rieman fit
821         "yr="<<paramYR[0]<<         // track position y - rieman fit
822         "zr="<<paramZR[0]<<         // track position z - rieman fit 
823         "yrW="<<paramYRW[0]<<         // track position y - rieman fit - weight
824         "zrW="<<paramZRW[0]<<         // track position z - rieman fit - weight
825         "dyr="<<paramYR[1]<<         // track position y - rieman fit
826         "dzr="<<paramZR[1]<<         // track position z - rieman fit 
827         "dyrW="<<paramYRW[1]<<         // track position y - rieman fit - weight
828         "dzrW="<<paramZRW[1]<<         // track position z - rieman fit - weight
829         //
830         "angley="<<angley<<         // angle par fit
831         "anglez="<<anglez<<         // angle par fit
832         "zdr="<<zdrift<<            //
833         "dy="<<deltay<<
834         "dz="<<deltaz<<        
835         "sy="<<csigmaY<<
836         "sz="<<csigmaZ<<
837         "chi2R="<<chi2R<<
838         "chi2RW="<<chi2RW<<
839         "\n";
840     }    
841     if (reject) continue;
842
843     
844     // =========================================
845     // wirte collected information to histograms
846     // =========================================
847         
848     Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
849     fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
850     
851     
852     TH3F * his3 = 0;
853     his3 = (TH3F*)fRMSY->At(padSize);
854     if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) ));
855     his3 = (TH3F*)fRMSZ->At(padSize);
856     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) ));
857     
858     his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
859     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaY2()) ));
860     his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
861     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(TMath::Abs(cluster0->GetSigmaZ2()) ));
862     
863     
864     his3 = (TH3F*)fResolY->At(padSize);
865     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
866     his3 = (TH3F*)fResolZ->At(padSize);
867     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
868     his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
869     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
870     his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
871     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );        
872     //=============================================================================================    
873     //
874     // Fill THN histograms
875     //
876     Double_t scaleQ= TMath::Min(Double_t(cluster0->GetMax()),200.)/30.;
877     Bool_t   isSelected = (TMath::Exp(scaleQ)>fQDownscaleRatio*gRandom->Rndm());
878     if (!isSelected) continue;
879     Double_t xvar[9];
880     xvar[1]=padSize;   // pad type 
881     xvar[2]=cluster0->GetZ();  // 
882     xvar[3]=cluster0->GetMax();
883     
884     xvar[0]=deltay;
885     double cog = cluster0->GetPad()-Int_t(cluster0->GetPad());// distance to the center of the pad       
886     xvar[4] = cog;
887     xvar[5]=angley;    
888
889     if( TMath::Abs(cog-0.5)<1.e-8 ) xvar[4]=-1; // fill one-pad clusters in the underflow bin
890     fHisDeltaY->Fill(xvar);       
891
892     xvar[4]=cog;
893     xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
894     fHisRMSY->Fill(xvar);
895     
896     xvar[0]=deltaz;
897     xvar[4]=cluster0->GetTimeBin()-Int_t(cluster0->GetTimeBin()); // distance to the center of the time bin
898     xvar[5]=anglez;
899     fHisDeltaZ->Fill(xvar);
900     xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
901     fHisRMSZ->Fill(xvar);
902     
903   }    // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
904 }  // FillResolutionHistoLocal(...)
905
906
907
908
909
910
911
912
913 void  AliTPCcalibTracks::SetStyle() const {
914    // 
915    // set style, can be called by all draw functions
916    // 
917    gROOT->SetStyle("Plain");
918    gStyle->SetFillColor(10);
919    gStyle->SetPadColor(10);
920    gStyle->SetCanvasColor(10);
921    gStyle->SetStatColor(10);
922    gStyle->SetPalette(1,0);
923    gStyle->SetNumberContours(60);
924 }
925
926
927
928 void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){ 
929    // 
930    // all functions are called, that produce pictures
931    // the histograms are written to the directory 'pathName'
932    // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
933    // 'stat' is also the number of minEntries for MakeResPlotsQTree
934    // 
935
936   if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
937   MakeResPlotsQTree(stat, pathName);
938 }
939    
940
941
942
943 void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
944    //
945    // Make tree with resolution parameters
946    // the result is written to 'resol.root' in directory 'pathname'
947    // file information are available in fileInfo
948    // available variables in the tree 'Resol':
949    //  Entries: number of entries for this resolution point
950    //  nbins:   number of bins that were accumulated
951    //  Dim:     direction, Dim==0: y-direction, Dim==1: z-direction
952    //  Pad:     padSize; short, medium and long
953    //  Length:  pad length, 0.75, 1, 1.5
954    //  QMean:   mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
955    //  Zc:      center of middle bin in drift direction
956    //  Zm:      mean dirftlength for accumulated Delta-Histograms
957    //  Zs:      width of driftlength bin
958    //  AngleC:  center of middle bin in Angle-Direction
959    //  AngleM:  mean angle for accumulated Delta-Histograms
960    //  AngleS:  width of Angle-bin
961    //  Resol:   sigma for gaus fit through Delta-Histograms
962    //  Sigma:   error of sigma for gaus fit through Delta Histograms
963    //  MeanR:   mean of the Delta-Histogram
964    //  SigmaR:  rms of the Delta-Histogram
965    //  RMSm:    mean of the gaus fit through RMS-Histogram
966    //  RMS:     sigma of the gaus fit through RMS-Histogram
967    //  RMSe0:   error of mean of gaus fit in RMS-Histogram
968    //  RMSe1:   error of sigma of gaus fit in RMS-Histogram
969    //  
970       
971   if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
972   if (GetDebugLevel() > -1) cout << "    relax, the calculation will take a while..." << endl;
973   
974    gSystem->MakeDirectory(pathName);
975    gSystem->ChangeDirectory(pathName);
976    TString kFileName = "resol.root";
977    TTreeSRedirector fTreeResol(kFileName.Data());
978    
979    TH3F *resArray[2][3][11];
980    TH3F *rmsArray[2][3][11];
981   
982    // load histograms from fArrayQDY and fArrayQDZ 
983    // into resArray and rmsArray
984    // that is all we need here
985    for (Int_t idim = 0; idim < 2; idim++){
986       for (Int_t ipad = 0; ipad < 3; ipad++){
987          for (Int_t iq = 0; iq <= 10; iq++){
988             rmsArray[idim][ipad][iq]=0;
989             resArray[idim][ipad][iq]=0;
990             Int_t bin = GetBin(iq,ipad); 
991             TH3F *hresl = 0;
992             if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
993             if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
994             if (!hresl) continue;
995             resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
996             resArray[idim][ipad][iq]->SetDirectory(0);
997             TH3F * hreslRMS = 0;
998             if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
999             if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
1000             if (!hreslRMS) continue;
1001             rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
1002             rmsArray[idim][ipad][iq]->SetDirectory(0);
1003          }
1004       }
1005    }
1006    if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
1007    
1008    //--------------------------------------------------------------------------------------------
1009    
1010    char name[200];
1011    Double_t qMean;
1012    Double_t zMean, angleMean, zCenter, angleCenter;
1013    Double_t zSigma, angleSigma;
1014    TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
1015    TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
1016    TF1 *fitFunction = new TF1("fitFunction", "gaus");
1017    Float_t entriesQ = 0;
1018    Int_t loopCounter = 1;
1019   
1020    for (Int_t idim = 0; idim < 2; idim++){
1021       // Loop y-z corrdinate
1022       for (Int_t ipad = 0; ipad < 3; ipad++){
1023          // loop pad type
1024          for (Int_t iq = -1; iq < 10; iq++){
1025             // LOOP Q
1026            if (GetDebugLevel() > -1) 
1027                cout << "Loop-counter, this is loop " << loopCounter << " of 66, (" 
1028                   << (Int_t)((loopCounter)/66.*100) << "% done), " 
1029                   << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << "  \r" << std::flush;
1030             loopCounter++;
1031             TH3F *hres = 0;
1032             TH3F *hrms = 0;
1033             qMean    = 0;
1034             entriesQ = 0;
1035             
1036             // calculate qMean
1037             if (iq == -1){
1038                // integrated spectra
1039                for (Int_t iql = 0; iql < 10; iql++){    
1040                   Int_t bin = GetBin(iql,ipad); 
1041                   TH3F *hresl = resArray[idim][ipad][iql];
1042                   TH3F *hrmsl = rmsArray[idim][ipad][iql];
1043                   if (!hresl) continue;
1044                   if (!hrmsl) continue;     
1045                   entriesQ += hresl->GetEntries();
1046                   qMean += hresl->GetEntries() * GetQ(bin);      
1047                   if (!hres) {
1048                      hres = (TH3F*)hresl->Clone();
1049                      hrms = (TH3F*)hrmsl->Clone();
1050                   }
1051                   else{
1052                      hres->Add(hresl);
1053                      hrms->Add(hrmsl);
1054                   }
1055                }
1056                qMean /= entriesQ;
1057                qMean *= -1.;  // integral mean charge
1058             }
1059             else {
1060                // loop over neighboured Q-bins 
1061                // accumulate entries from neighboured Q-bins
1062                for (Int_t iql = iq - 1; iql <= iq + 1; iql++){              
1063                   if (iql < 0) continue;
1064                   Int_t bin = GetBin(iql,ipad);
1065                   TH3F * hresl = resArray[idim][ipad][iql];
1066                   TH3F * hrmsl = rmsArray[idim][ipad][iql];
1067                   if (!hresl) continue;
1068                   if (!hrmsl) continue;
1069                   entriesQ += hresl->GetEntries(); 
1070                   qMean += hresl->GetEntries() * GetQ(bin);      
1071                   if (!hres) {
1072                      hres = (TH3F*) hresl->Clone();
1073                      hrms = (TH3F*) hrmsl->Clone();
1074                   }
1075                   else{
1076                      hres->Add(hresl);
1077                      hrms->Add(hrmsl);
1078                   }
1079                }
1080                qMean/=entriesQ;
1081             }
1082             if (!hres) continue;
1083             if (!hrms) continue;
1084
1085             TAxis *xAxisDriftLength = hres->GetXaxis();   // driftlength / z - axis
1086             TAxis *yAxisAngle       = hres->GetYaxis();   // angle axis
1087             TAxis *zAxisDelta       = hres->GetZaxis();   // delta axis
1088             TAxis *zAxisRms         = hrms->GetZaxis();   // rms axis
1089             
1090             // loop over all angle bins
1091             for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
1092                angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
1093                // loop over all driftlength bins
1094                for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
1095                   zCenter    = xAxisDriftLength->GetBinCenter(ibinxDL);
1096                   zSigma     = xAxisDriftLength->GetBinWidth(ibinxDL);
1097                   angleSigma = yAxisAngle->GetBinWidth(ibinyAngle); 
1098                   zMean      = zCenter;      // changens, when more statistic is accumulated
1099                   angleMean  = angleCenter;  // changens, when more statistic is accumulated
1100                   
1101                   // create 2 1D-Histograms, projectionRes and projectionRms
1102                   // these histograms are delta histograms for given direction, padSize, chargeBin,
1103                   // angleBin and driftLengthBin
1104                   // later on they will be fitted with a gausian, its sigma is the resoltuion...
1105                   snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
1106                   // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1107                   projectionRes->SetNameTitle(name, name);
1108                   snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
1109                   // TH1D * projectionRms =  new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1110                   projectionRms->SetNameTitle(name, name);
1111                   
1112                   projectionRes->Reset();
1113                   projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1114                   projectionRms->Reset();
1115                   projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1116                   projectionRes->SetDirectory(0);
1117                   projectionRms->SetDirectory(0);
1118                   
1119                   Double_t entries = 0;
1120                   Int_t    nbins   = 0;   // counts, how many bins were accumulated
1121                   zMean     = 0;
1122                   angleMean = 0;
1123                   entries   = 0;
1124                   
1125                   // fill projectionRes and projectionRms for given dim, ipad and iq, 
1126                   // as well as for given angleBin and driftlengthBin
1127                   // if this gives not enough statistic, include neighbourhood 
1128                   // (angle and driftlength) successifely
1129                   for (Int_t dbin = 0; dbin <= 8; dbin++){              // delta-bins around centered angleBin and driftlengthBin
1130                      for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) {   // delta-bins in angle direction
1131                         for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
1132                            if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue;   // add each bin only one time !
1133                            Int_t binx2 = ibinxDL + dbinx2;                       // position variable in x (driftlength) direction
1134                            Int_t biny2 = ibinyAngle + dbiny2;                    // position variable in y (angle)  direction
1135                            if (binx2 < 1 || biny2 < 1) continue;                 // don't go out of the histogram!
1136                            if (binx2 >= xAxisDriftLength->GetNbins()) continue;  // don't go out of the histogram!
1137                            if (biny2 >= yAxisAngle->GetNbins()) continue;        // don't go out of the histogram!
1138                            nbins++;                                              // count the number of accumulated bins
1139                            // Fill resolution histo
1140                            for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
1141                               // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3);     // unused variable
1142                               projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
1143                               entries   += hres->GetBinContent(binx2, biny2, ibin3);
1144                               zMean     += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
1145                               angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
1146                            }  // ibin3 loop
1147                            // fill RMS histo
1148                            for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
1149                               projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
1150                            }
1151                         }  //dbinx2 loop
1152                         if (entries > minEntries) break; // enough statistic accumulated
1153                      }  // dbiny2 loop
1154                      if (entries > minEntries) break;    // enough statistic accumulated
1155                   }  // dbin loop
1156                   if ( entries< minEntries) continue;  // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree  
1157                   zMean /= entries;
1158                   angleMean /= entries;
1159                   
1160                   if (entries > minEntries) {
1161                      //  when enough statistic is accumulated
1162                      //  fit Delta histograms with a gausian
1163                      //  of the gausian is the resolution (resol), its fit error is sigma
1164                      //  store also mean and RMS of the histogram
1165                      Float_t xmin     = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1166                      Float_t xmax     = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1167                      
1168 //                      projectionRes->Fit("gaus", "q0", "", xmin, xmax);
1169 //                      Float_t resol    = projectionRes->GetFunction("gaus")->GetParameter(2);
1170 //                      Float_t sigma    = projectionRes->GetFunction("gaus")->GetParError(2);
1171                      fitFunction->SetMaximum(xmax);
1172                      fitFunction->SetMinimum(xmin);
1173                      projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
1174                      Float_t resol    = fitFunction->GetParameter(2);
1175                      Float_t sigma    = fitFunction->GetParError(2);
1176                      
1177                      Float_t meanR    = projectionRes->GetMean();
1178                      Float_t sigmaR   = projectionRes->GetRMS();
1179                      // fit also RMS histograms with a gausian
1180                      // store mean and sigma of the gausian in rmsMean and rmsSigma
1181                      // store also the fit errors in errorRMS and errorSigma
1182                      xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1183                      xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1184                      
1185 //                      projectionRms->Fit("gaus","q0","",xmin,xmax);
1186 //                      Float_t rmsMean    = projectionRms->GetFunction("gaus")->GetParameter(1);
1187 //                      Float_t rmsSigma   = projectionRms->GetFunction("gaus")->GetParameter(2);
1188 //                      Float_t errorRMS   = projectionRms->GetFunction("gaus")->GetParError(1);
1189 //                      Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
1190                      projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
1191                      Float_t rmsMean    = fitFunction->GetParameter(1);
1192                      Float_t rmsSigma   = fitFunction->GetParameter(2);
1193                      Float_t errorRMS   = fitFunction->GetParError(1);
1194                      Float_t errorSigma = fitFunction->GetParError(2);
1195                     
1196                      Float_t length = 0.75;
1197                      if (ipad == 1) length = 1;
1198                      if (ipad == 2) length = 1.5;
1199                      
1200                      fTreeResol<<"Resol"<<
1201                         "Entries="<<entries<<      // number of entries for this resolution point
1202                         "nbins="<<nbins<<          // number of bins that were accumulated
1203                         "Dim="<<idim<<             // direction, Dim==0: y-direction, Dim==1: z-direction
1204                         "Pad="<<ipad<<             // padSize; short, medium and long
1205                         "Length="<<length<<        // pad length, 0.75, 1, 1.5
1206                         "QMean="<<qMean<<          // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1207                         "Zc="<<zCenter<<           // center of middle bin in drift direction
1208                         "Zm="<<zMean<<             // mean dirftlength for accumulated Delta-Histograms
1209                         "Zs="<<zSigma<<            // width of driftlength bin
1210                         "AngleC="<<angleCenter<<   // center of middle bin in Angle-Direction
1211                         "AngleM="<<angleMean<<     // mean angle for accumulated Delta-Histograms
1212                         "AngleS="<<angleSigma<<    // width of Angle-bin
1213                         "Resol="<<resol<<          // sigma for gaus fit through Delta-Histograms
1214                         "Sigma="<<sigma<<          // error of sigma for gaus fit through Delta Histograms
1215                         "MeanR="<<meanR<<          // mean of the Delta-Histogram
1216                         "SigmaR="<<sigmaR<<        // rms of the Delta-Histogram
1217                         "RMSm="<<rmsMean<<         // mean of the gaus fit through RMS-Histogram
1218                         "RMSs="<<rmsSigma<<        // sigma of the gaus fit through RMS-Histogram
1219                         "RMSe0="<<errorRMS<<       // error of mean of gaus fit in RMS-Histogram
1220                         "RMSe1="<<errorSigma<<     // error of sigma of gaus fit in RMS-Histogram
1221                         "\n";
1222                      if (GetDebugLevel() > 5) {
1223                         projectionRes->SetDirectory(fTreeResol.GetFile());
1224                         projectionRes->Write(projectionRes->GetName());
1225                         projectionRes->SetDirectory(0);
1226                         projectionRms->SetDirectory(fTreeResol.GetFile());
1227                         projectionRms->Write(projectionRms->GetName());
1228                         projectionRes->SetDirectory(0);
1229                      }
1230                   }  // if (projectionRes->GetSum() > minEntries)
1231                }  // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
1232             }  // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
1233             
1234          }  // iq-loop
1235       }  // ipad-loop
1236    }  // idim-loop
1237    delete projectionRes;
1238    delete projectionRms;
1239    
1240 //    TFile resolFile(fTreeResol.GetFile());
1241    TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
1242    fileInfo.Write("fileInfo");
1243 //    resolFile.Close();
1244 //    fTreeResol.GetFile()->Close();
1245    if (GetDebugLevel() > -1) cout << endl;
1246    if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
1247    gSystem->ChangeDirectory("..");
1248 }
1249
1250
1251
1252
1253
1254 Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
1255    // 
1256    // function to merge several AliTPCcalibTracks objects after PROOF calculation
1257    // The object's histograms are merged via their merge functions
1258    // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
1259    // 
1260    
1261   if (GetDebugLevel() > 0) cout << " *****  this is AliTPCcalibTracks::Merge(TCollection *collectionList)  *****"<< endl;  
1262    if (!collectionList) return 0;
1263    if (collectionList->IsEmpty()) return -1;
1264    
1265    if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl;     //    REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
1266    if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
1267    collectionList->Print();
1268    
1269    // create a list for each data member
1270    TList* deltaYList = new TList;
1271    TList* deltaZList = new TList;
1272    TList* arrayAmpRowList = new TList;
1273    TList* rejectedTracksList = new TList;
1274    TList* clusterCutHistoList = new TList;
1275    TList* arrayAmpList = new TList;
1276    TList* arrayQDYList = new TList;
1277    TList* arrayQDZList = new TList;
1278    TList* arrayQRMSYList = new TList;
1279    TList* arrayQRMSZList = new TList;
1280    TList* resolYList = new TList;
1281    TList* resolZList = new TList;
1282    TList* rMSYList = new TList;
1283    TList* rMSZList = new TList;
1284    
1285 //    TList* nRowsList = new TList;
1286 //    TList* nSectList = new TList;
1287 //    TList* fileNoList = new TList;
1288    
1289    TIterator *listIterator = collectionList->MakeIterator();
1290    AliTPCcalibTracks *calibTracks = 0;
1291    if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;    
1292    Int_t counter = 0;
1293    while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
1294       // loop over all entries in the collectionList and get dataMembers into lists
1295       
1296       arrayQDYList->Add(calibTracks->GetfArrayQDY());
1297       arrayQDZList->Add(calibTracks->GetfArrayQDZ());
1298       arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
1299       arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
1300       resolYList->Add(calibTracks->GetfResolY());
1301       resolZList->Add(calibTracks->GetfResolZ());
1302       rMSYList->Add(calibTracks->GetfRMSY());
1303       rMSZList->Add(calibTracks->GetfRMSZ());
1304       rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
1305       clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
1306       //
1307       if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
1308         fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());      
1309       //      fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
1310       counter++;
1311       if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
1312       AddHistos(calibTracks);
1313    }
1314    
1315    
1316    // merge data members
1317    if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl; 
1318    fClusterCutHisto->Merge(clusterCutHistoList);
1319    fRejectedTracksHisto->Merge(rejectedTracksList);
1320    
1321    TObjArray* objarray = 0;
1322    TH1* hist = 0;
1323    TList* histList = 0;
1324    TIterator *objListIterator = 0;
1325    
1326       
1327    if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
1328    // merge fArrayQDY
1329    for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
1330       objListIterator = arrayQDYList->MakeIterator();
1331       histList = new TList;
1332       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1333          // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
1334          hist = (TH3F*)(objarray->At(i));
1335          histList->Add(hist);
1336       }
1337       ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
1338       delete histList;
1339       delete objListIterator;
1340    }
1341
1342    if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
1343    // merge fArrayQDZ
1344    for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1345       objListIterator = arrayQDZList->MakeIterator();
1346       histList = new TList;
1347       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1348          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1349          hist = (TH3F*)(objarray->At(i));
1350          histList->Add(hist);
1351       }
1352       ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
1353       delete histList;
1354       delete objListIterator;
1355    }
1356
1357    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
1358    // merge fArrayQRMSY
1359    for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
1360       objListIterator = arrayQRMSYList->MakeIterator();
1361       histList = new TList;
1362       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1363          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1364          //   if (!objarray) continue; // removed for coverity -> JMT
1365          hist = (TH3F*)(objarray->At(i));
1366          histList->Add(hist);
1367       }
1368       ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
1369       delete histList;
1370       delete objListIterator;
1371    }   
1372
1373    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
1374    // merge fArrayQRMSZ
1375    for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1376       objListIterator = arrayQRMSZList->MakeIterator();
1377       histList = new TList;
1378       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1379          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1380          hist = (TH3F*)(objarray->At(i));
1381          histList->Add(hist);
1382       }
1383       ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
1384       delete histList;
1385       delete objListIterator;
1386    } 
1387    
1388    
1389   
1390    
1391         
1392    
1393    if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
1394    // merge fResolY
1395    for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
1396       objListIterator = resolYList->MakeIterator();
1397       histList = new TList;
1398       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1399          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1400          hist = (TH3F*)(objarray->At(i));
1401          histList->Add(hist);
1402       }
1403       ((TH3F*)(fResolY->At(i)))->Merge(histList);
1404       delete histList;
1405       delete objListIterator;
1406    }
1407    
1408    // merge fResolZ
1409    for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1410       objListIterator = resolZList->MakeIterator();
1411       histList = new TList;
1412       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1413          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1414          hist = (TH3F*)(objarray->At(i));
1415          histList->Add(hist);
1416       }
1417       ((TH3F*)(fResolZ->At(i)))->Merge(histList);
1418       delete histList;
1419       delete objListIterator;
1420    }
1421
1422    // merge fRMSY
1423    for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
1424       objListIterator = rMSYList->MakeIterator();
1425       histList = new TList;
1426       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1427          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1428          hist = (TH3F*)(objarray->At(i));
1429          histList->Add(hist);
1430       }
1431       ((TH3F*)(fRMSY->At(i)))->Merge(histList);
1432       delete histList;
1433       delete objListIterator;
1434    }
1435          
1436    // merge fRMSZ
1437    for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1438       objListIterator = rMSZList->MakeIterator();
1439       histList = new TList;
1440       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1441          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1442          hist = (TH3F*)(objarray->At(i));
1443          histList->Add(hist);
1444       }
1445       ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
1446       delete histList;
1447       delete objListIterator;
1448    }
1449    
1450    delete deltaYList;
1451    delete deltaZList;
1452    delete arrayAmpRowList;
1453    delete arrayAmpList;
1454    delete arrayQDYList;
1455    delete arrayQDZList;
1456    delete arrayQRMSYList;
1457    delete arrayQRMSZList;
1458    delete resolYList;
1459    delete resolZList;
1460    delete rMSYList;
1461    delete rMSZList;
1462    delete listIterator;
1463    
1464    if (GetDebugLevel() > 0) cout << "merging done!" << endl;
1465    
1466    return 1;
1467 }
1468
1469
1470
1471
1472 void AliTPCcalibTracks::MakeHistos(){
1473   //
1474   ////make THnSparse
1475   //
1476   //THnSparse  *fHisDeltaY;    // THnSparse - delta Y 
1477   //THnSparse  *fHisDeltaZ;    // THnSparse - delta Z 
1478   //THnSparse  *fHisRMSY;      // THnSparse - rms Y 
1479   //THnSparse  *fHisRMSZ;      // THnSparse - rms Z 
1480   //THnSparse  *fHisQmax;      // THnSparse - qmax 
1481   //THnSparse  *fHisQtot;      // THnSparse - qtot 
1482   // cluster  performance bins
1483   // 0 - variable of interest
1484   // 1 - pad type   - 0- short 1-medium 2-long pads
1485   // 2 - drift length - drift length -0-1
1486   // 3 - Qmax         - Qmax  - 2- 400
1487   // 4 - cog          - COG position - 0-1
1488   // 5 - tan(phi)     - local y angle
1489   // 6 - tan(theta)   - local z angle
1490   // 7 - sector       - sector number
1491   Double_t xminTrack[8], xmaxTrack[8];
1492   Int_t binsTrack[8];
1493   TString axisName[8];
1494   
1495   //
1496   binsTrack[0]=200;
1497   axisName[0]  ="var";
1498
1499   binsTrack[1] =3;
1500   xminTrack[1] =0; xmaxTrack[1]=3;
1501   axisName[1]  ="pad type";
1502   //
1503   binsTrack[2] =20;
1504   xminTrack[2] =-250; xmaxTrack[2]=250;
1505   axisName[2]  ="z";
1506   //
1507   binsTrack[3] =10;
1508   xminTrack[3] =1; xmaxTrack[3]=400;
1509   axisName[3]  ="Qmax";
1510   //
1511   binsTrack[4] =20;
1512   xminTrack[4] =0; xmaxTrack[4]=1;
1513   axisName[4]  ="cog";
1514   //
1515   binsTrack[5] =15;
1516   xminTrack[5] =-1.5; xmaxTrack[5]=1.5;
1517   axisName[5]  ="tan(angle)";
1518   //
1519   //
1520   xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1521   fHisDeltaY=new THnSparseF("#Delta_{y} (cm)","#Delta_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1522   xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1523   fHisDeltaZ=new THnSparseF("#Delta_{z} (cm)","#Delta_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1524   xminTrack[0] =0.; xmaxTrack[0]=0.5;
1525   fHisRMSY=new THnSparseF("#RMS_{y} (cm)","#RMS_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1526   xminTrack[0] =0.; xmaxTrack[0]=0.5;
1527   fHisRMSZ=new THnSparseF("#RMS_{z} (cm)","#RMS_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1528   xminTrack[0] =0.; xmaxTrack[0]=100;
1529   fHisQmax=new THnSparseF("Qmax (ADC)","Qmax (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
1530
1531   xminTrack[0] =0.; xmaxTrack[0]=250;
1532   fHisQtot=new THnSparseF("Qtot (ADC)","Qtot (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
1533
1534
1535   for (Int_t ivar=0;ivar<6;ivar++){
1536     fHisDeltaY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1537     fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1538     fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1539     fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1540     fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1541     fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1542     fHisDeltaY->GetAxis(ivar)->SetTitle(axisName[ivar].Data());
1543     fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1544     fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1545     fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1546     fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1547     fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1548   }
1549
1550
1551   BinLogX(fHisDeltaY,3);
1552   BinLogX(fHisDeltaZ,3);
1553   BinLogX(fHisRMSY,3);
1554   BinLogX(fHisRMSZ,3);
1555   BinLogX(fHisQmax,3);
1556   BinLogX(fHisQtot,3);
1557
1558 }  
1559
1560 void    AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){
1561   //
1562   // Add histograms
1563   //
1564   if (!calib->fHisDeltaY) return;
1565   if (calib->fHisDeltaY->GetEntries()> fgkMergeEntriesCut) return; 
1566   if (calib->fHisDeltaY) fHisDeltaY->Add(calib->fHisDeltaY);
1567   if (calib->fHisDeltaZ) fHisDeltaZ->Add(calib->fHisDeltaZ);
1568   if (calib->fHisRMSY)   fHisRMSY->Add(calib->fHisRMSY);
1569   if (calib->fHisRMSZ)   fHisRMSZ->Add(calib->fHisRMSZ);
1570 }
1571
1572
1573
1574 void AliTPCcalibTracks::MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype){
1575   //
1576   // Dump summary info
1577   //
1578   //  0.OBJ: TAxis     var
1579   //  1.OBJ: TAxis     pad type
1580   //  2.OBJ: TAxis     z
1581   //  3.OBJ: TAxis     Qmax
1582   //  4.OBJ: TAxis     cog
1583   //  5.OBJ: TAxis     tan(angle)
1584   //
1585   if (ptype>3) return;
1586   Int_t idim[6]={0,1,2,3,4,5};
1587   TString hname[4]={"dy","dz","rmsy","rmsz"};
1588   //
1589   Int_t nbins5=hisInput->GetAxis(5)->GetNbins();
1590   Int_t first5=hisInput->GetAxis(5)->GetFirst();
1591   Int_t last5 =hisInput->GetAxis(5)->GetLast();
1592   //
1593   for (Int_t ibin5=first5-1; ibin5<=last5; ibin5+=1){   // axis 5 - local angle
1594     Bool_t bin5All=kTRUE;
1595     if (ibin5>=first5){
1596       hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5-1,first5),TMath::Min(ibin5+1,nbins5));
1597       if (ibin5==first5) hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5,first5),TMath::Min(ibin5,nbins5));
1598       bin5All=kFALSE;
1599     }
1600     Double_t      x5= hisInput->GetAxis(5)->GetBinCenter(ibin5);
1601     THnSparse * his5= hisInput->Projection(5,idim);         //projected histogram according selection 5    
1602     Int_t nbins4=his5->GetAxis(4)->GetNbins();
1603     Int_t first4=his5->GetAxis(4)->GetFirst();
1604     Int_t last4 =his5->GetAxis(4)->GetLast();
1605     
1606     for (Int_t ibin4=first4-1; ibin4<=last4; ibin4+=1){   // axis 4 - cog
1607       Bool_t bin4All=kTRUE;
1608       if (ibin4>=first4){
1609         his5->GetAxis(4)->SetRange(TMath::Max(ibin4+1,first4),TMath::Min(ibin4-1,nbins4));
1610         if (ibin4==first4||ibin4==last4) his5->GetAxis(4)->SetRange(TMath::Max(ibin4,first4),TMath::Min(ibin4,nbins4));
1611         bin4All=kFALSE;
1612       }
1613       Double_t      x4= his5->GetAxis(4)->GetBinCenter(ibin4);
1614       THnSparse * his4= his5->Projection(4,idim);         //projected histogram according selection 4
1615       //
1616       Int_t nbins3=his4->GetAxis(3)->GetNbins();
1617       Int_t first3=his4->GetAxis(3)->GetFirst();
1618       Int_t last3 =his4->GetAxis(3)->GetLast();
1619       //
1620       for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){   // axis 3 - Qmax
1621         Bool_t bin3All=kTRUE;
1622         if (ibin3>=first3){
1623           his4->GetAxis(3)->SetRange(TMath::Max(ibin3,first3),TMath::Min(ibin3,nbins3));
1624           bin3All=kFALSE;
1625         }
1626         Double_t      x3= his4->GetAxis(3)->GetBinCenter(ibin3);
1627         THnSparse * his3= his4->Projection(3,idim);         //projected histogram according selection 3
1628         //
1629         Int_t nbins2    = his3->GetAxis(2)->GetNbins();
1630         Int_t first2    = his3->GetAxis(2)->GetFirst();
1631         Int_t last2     = his3->GetAxis(2)->GetLast();
1632         //
1633         for (Int_t ibin2=first2-1; ibin2<=last2; ibin2+=1){   // axis 2 - z     
1634           Bool_t bin2All=kTRUE;
1635           Double_t      x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1636           if (ibin2>=first2){
1637             his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,first2),TMath::Min(ibin2+1,nbins2));
1638             if (ibin2==first2||ibin2==last2||TMath::Abs(x2)<20) his3->GetAxis(2)->SetRange(TMath::Max(ibin2,first2),TMath::Min(ibin2,nbins2));
1639             bin2All=kFALSE;
1640           }
1641           THnSparse * his2= his3->Projection(2,idim);         //projected histogram according selection 2
1642           //
1643           Int_t nbins1     = his2->GetAxis(1)->GetNbins();
1644           Int_t first1     = his2->GetAxis(1)->GetFirst();
1645           Int_t last1      = his2->GetAxis(1)->GetLast();
1646           for (Int_t ibin1=first1-1; ibin1<=last1; ibin1++){   //axis 1 - pad type 
1647             Bool_t bin1All=kTRUE;
1648             if (ibin1>=first1){
1649               his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));           
1650               bin1All=kFALSE;
1651             }
1652             Double_t       x1= TMath::Nint(his2->GetAxis(1)->GetBinCenter(ibin1)-0.5);
1653             TH1 * hisDelta = his2->Projection(0);
1654             Double_t entries = hisDelta->GetEntries();
1655             Double_t mean=0, rms=0;
1656             if (entries>10){
1657               mean    = hisDelta->GetMean();
1658               rms = hisDelta->GetRMS();   
1659               hisDelta->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
1660               mean    = hisDelta->GetMean();
1661               rms = hisDelta->GetRMS();
1662             }
1663             //
1664             (*pcstream)<<hname[ptype].Data()<<
1665               // flag - intgrated values for given bin
1666               "angleA="<<bin5All<<   
1667               "cogA="<<bin4All<<
1668               "qmaxA="<<bin3All<<
1669               "zA="<<bin2All<<
1670               "ipadA="<<bin1All<<
1671               // center of bin value
1672               "angle="<<x5<<       // local angle
1673               "cog="<<x4<<         // distance cluster to center
1674               "qmax="<<x3<<        // qmax
1675               "z="<<x2<<           // z of the cluster
1676               "ipad="<<x1<<        // type of the region
1677               // mean values
1678               //
1679               "entries="<<entries<<
1680               "mean="<<mean<<
1681               "rms="<<rms<<
1682               "ptype="<<ptype<<   //
1683               "\n";
1684             delete hisDelta;
1685             printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n",x5,x4,x3,x2,x1, entries,mean);          
1686           }//loop z
1687           delete his2;
1688         } // loop Q max
1689         delete his3;
1690       } // loop COG
1691       delete his4;
1692     }//loop local angle
1693     delete his5;
1694   }
1695 }
1696
1697
1698 int AliTPCcalibTracks::GetTHnStat(const  THnBase *H, THnBase *&Mean, THnBase *&Sigma, THnBase *&Entr )
1699 {
1700   // H is THnSparse:
1701   //
1702   // OBJ: TAxis     var     var
1703   // OBJ: TAxis     axis 1
1704   // OBJ: TAxis     axis 2
1705   // ...  
1706
1707   // Output: 
1708
1709   // Mean, Sigma and Entr is THnF of mean, RMS and entries of var:
1710   //
1711   // OBJ: TAxis     axis 1
1712   // OBJ: TAxis     axis 2
1713   // ..
1714     
1715   Mean = 0;
1716   Sigma = 0;
1717   Entr = 0;
1718   if( !H ) return -1;
1719
1720   Bool_t ok = kTRUE;
1721   
1722   int nDim = H->GetNdimensions();
1723   Long64_t nFilledBins = H->GetNbins();
1724   Long64_t nStatBins = 0;
1725
1726   Int_t *nBins = new Int_t [nDim];
1727   Double_t *xMin = new Double_t [nDim]; 
1728   Double_t *xMax = new Double_t [nDim];
1729   Int_t *idx = new Int_t [nDim];
1730
1731   TString nameMean = H->GetName();
1732   TString nameSigma = H->GetName();
1733   TString nameEntr = H->GetName();
1734   nameMean+="_Mean";
1735   nameSigma+="_Sigma";
1736   nameEntr+="_Entr";
1737
1738   ok = ok && ( nBins && xMin && xMax && idx );
1739
1740   if( ok ){
1741     for( int i=0; i<nDim; i++){
1742       xMin[i] = H->GetAxis(i)->GetXmin();
1743       xMax[i] = H->GetAxis(i)->GetXmax();
1744       nBins[i] = H->GetAxis(i)->GetNbins();
1745     }
1746   
1747     Mean  = new THnSparseF( nameMean.Data(),  nameMean.Data(),  nDim-1, nBins+1, xMin+1, xMax+1 );
1748     Sigma = new THnSparseF( nameSigma.Data(), nameSigma.Data(), nDim-1, nBins+1, xMin+1, xMax+1 );
1749     Entr  = new THnSparseF( nameEntr.Data(),  nameEntr.Data(),  nDim-1, nBins+1, xMin+1, xMax+1 );
1750   }
1751  
1752   ok = ok && ( Mean && Sigma && Entr );
1753
1754   if( ok ){
1755     for( int i=0; i<nDim-1; i++){
1756       const TAxis *ax = H->GetAxis(i+1);
1757       Mean->GetAxis(i)->SetName(ax->GetName());
1758       Mean->GetAxis(i)->SetTitle(ax->GetTitle());
1759       Sigma->GetAxis(i)->SetName(ax->GetName());
1760       Sigma->GetAxis(i)->SetTitle(ax->GetTitle());
1761       Entr->GetAxis(i)->SetName(ax->GetName());
1762       Entr->GetAxis(i)->SetTitle(ax->GetTitle());
1763       if( ax->GetXbins()->fN>0 ){
1764         Mean->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1765         Sigma->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1766         Entr->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1767       }
1768     } 
1769  
1770     // Allocate bins
1771  
1772     for( Long64_t binS=0; binS<nFilledBins; binS++){   
1773       H->GetBinContent(binS,idx);
1774       Mean->GetBin(idx+1,kTRUE);
1775       Sigma->GetBin(idx+1,kTRUE);
1776       Entr->GetBin(idx+1,kTRUE);
1777     }
1778   
1779     nStatBins = Mean->GetNbins();
1780       
1781   }
1782
1783  
1784   Long64_t *hMap = new Long64_t[nFilledBins];
1785   Double_t *hVal  = new Double_t[nFilledBins];
1786   Double_t *hEntr = new Double_t[nFilledBins];
1787   Double_t *meanD = new Double_t[nStatBins];
1788   Double_t *sigmaD = new Double_t[nStatBins];
1789
1790   ok = ok && ( hMap && hVal && hEntr && meanD && sigmaD );
1791
1792   // Map bins to Stat bins
1793
1794   if( ok ){ 
1795     for( Long64_t binS=0; binS<nFilledBins; binS++){   
1796       double val = H->GetBinContent(binS,idx);
1797       Long64_t bin = Mean->GetBin(idx+1,kFALSE);
1798       ok = ok && ( bin>=0 && bin<nStatBins && bin==Sigma->GetBin(idx+1,kFALSE) && bin== Entr->GetBin(idx+1,kFALSE) );
1799       if( !ok ) break;      
1800       if( val < 0 ){
1801         cout << "AliTPCcalibTracks: GetSparseStat : Unexpected content of the input THn"<<endl;
1802         bin = -1;
1803       }
1804       if( idx[0]<1 || idx[0]>nBins[0] ) bin = -1;
1805       hMap[binS] = bin;      
1806       hVal[binS] = idx[0];
1807       hEntr[binS] = val;
1808     }
1809   }
1810   
1811   // do 2 iteration with cut at 4 Sigma
1812
1813   for( int iter=0; ok && iter<2; iter++ ){
1814     
1815     // clean up entries
1816     
1817     for( Long64_t bin=0; bin < nStatBins; bin++){
1818       Entr->SetBinContent(bin, 0); 
1819       meanD[bin]=0;
1820       sigmaD[bin]=0;
1821     }
1822
1823     // get Entries and Mean
1824
1825     for( Long64_t binS=0; binS<nFilledBins; binS++){   
1826       Long64_t bin = hMap[binS];
1827       Double_t val  = hVal[binS];
1828       Double_t entr = hEntr[binS];
1829       if( bin < 0 ) continue;
1830       if( iter!=0 ){ // cut
1831         double s2 = Sigma->GetBinContent(bin);
1832         double d = val - Mean->GetBinContent(bin);
1833         if( d*d>s2*16 ) continue;
1834       } 
1835       meanD[bin]+= entr*val;
1836       Entr->AddBinContent(bin,entr);    
1837     }
1838   
1839     for( Long64_t bin=0; bin<nStatBins; bin++){
1840       double val = meanD[bin];
1841       double sum = Entr->GetBinContent(bin);
1842       if( sum>=1 ){
1843         Mean->SetBinContent(bin, val/sum );
1844       }
1845       else Mean->SetBinContent(bin, 0);
1846       Entr->SetBinContent(bin, 0); 
1847     }
1848
1849     // get RMS
1850
1851     for( Long64_t binS=0; binS<nFilledBins; binS++){   
1852       Long64_t bin = hMap[binS];
1853       Double_t val  = hVal[binS];
1854       Double_t entr = hEntr[binS];
1855       if( bin < 0 ) continue;  
1856       double d2 = val - Mean->GetBinContent(bin);
1857       d2 *= d2;
1858       if( iter!=0 ){ // cut
1859         double s2 = Sigma->GetBinContent(bin);
1860         if( d2>s2*16 ) continue;
1861       }
1862       sigmaD[bin] += entr*d2;
1863       Entr->AddBinContent(bin,entr);    
1864     }
1865
1866     for( Long64_t bin=0; bin<nStatBins; bin++ ){
1867       double val = sigmaD[bin];
1868       double sum = Entr->GetBinContent(bin);
1869       if( sum>=1 && val>=0 ){
1870         Sigma->SetBinContent(bin, val/sum );
1871       }
1872       else Sigma->SetBinContent(bin, 0);    
1873     }
1874   }
1875   
1876   // scale the Mean and the Sigma
1877
1878   if( ok ){
1879     double v0 = H->GetAxis(0)->GetBinCenter(1);
1880     double vscale = H->GetAxis(0)->GetBinWidth(1);
1881   
1882     for( Long64_t bin=0; bin<nStatBins; bin++){
1883       double m = Mean->GetBinContent(bin);
1884       double s2 = Sigma->GetBinContent(bin);
1885       double sum = Entr->GetBinContent(bin);
1886       if( sum>0 && s2>=0 ){
1887         Mean->SetBinContent(bin, v0 + (m-1)*vscale );
1888         Sigma->SetBinContent(bin, sqrt(s2)*vscale );
1889       }
1890       else{
1891         Mean->SetBinContent(bin, 0);
1892         Sigma->SetBinContent(bin, 0);
1893         Entr->SetBinContent(bin, 0);
1894       }
1895     }
1896   }
1897
1898   delete[] nBins;
1899   delete[] xMin;
1900   delete[] xMax;
1901   delete[] idx; 
1902   delete[] hMap;
1903   delete[] hVal;
1904   delete[] hEntr;
1905   delete[] meanD;
1906   delete[] sigmaD;
1907
1908   if( !ok ){
1909     cout << "AliTPCcalibTracks: GetSparseStat : No memory or internal Error "<<endl;
1910     delete Mean;
1911     delete Sigma;
1912     delete Entr;
1913     Mean =0;
1914     Sigma = 0;
1915     Entr = 0;
1916     return -1;
1917   }
1918   
1919   return 0;
1920 }
1921
1922
1923
1924 int AliTPCcalibTracks::CreateWaveCorrection(const  THnBase *DeltaY, THnBase *&MeanY, THnBase *&SigmaY, THnBase *&EntrY, 
1925                                             Bool_t MirrorZ, Bool_t MirrorPad, Bool_t MirrorAngle, Int_t MinStat  )
1926 {
1927   // DeltaY is THnSparse:
1928   //
1929   // OBJ: TAxis 0    var     var
1930   // OBJ: TAxis 1    pad type        pad type
1931   // OBJ: TAxis 2    z       z
1932   // OBJ: TAxis 3    Qmax    Qmax
1933   // OBJ: TAxis 4    cog     cog
1934   // OBJ: TAxis 5    tan(angle)      tan(angle)
1935   // ...  
1936
1937   MeanY = 0;
1938   SigmaY = 0;
1939   EntrY = 0;
1940
1941   int nDim = DeltaY->GetNdimensions();
1942   if( nDim<6 ){
1943     cout << "AliTPCcalibTracks: CreateWaveCorrection: Unknown input"<<endl;
1944     return -1;
1945   }
1946
1947   Int_t *nBins = new Int_t[nDim];
1948   Int_t *nBinsNew = new Int_t[nDim];
1949   Double_t *xMin = new Double_t[nDim]; 
1950   Double_t *xMax = new Double_t[nDim];
1951   Int_t *idx = new Int_t[nDim];
1952   THnSparseD *mergedDeltaY = 0;
1953
1954   int ret = 0;
1955   
1956   if( !nBins || !nBinsNew || !xMin || !xMax || !idx ){
1957     ret = -1;
1958     cout << "AliTPCcalibTracks: CreateWaveCorrection: Out of memory"<<endl;
1959   }
1960
1961   if( ret==0 ){
1962     
1963     for( int i=0; i<nDim; i++){
1964       xMin[i] = DeltaY->GetAxis(i)->GetXmin();
1965       xMax[i] = DeltaY->GetAxis(i)->GetXmax();
1966       nBins[i] = DeltaY->GetAxis(i)->GetNbins();
1967       nBinsNew[i] = nBins[i];
1968     }
1969   
1970     // Merge cog axis
1971     if( MirrorPad ){ 
1972       Int_t centralBin = DeltaY->GetAxis(4)->FindFixBin(0.5);
1973       xMin[4] = DeltaY->GetAxis(4)->GetBinLowEdge(centralBin);
1974       nBinsNew[4] = nBinsNew[4]-centralBin+1;   
1975     }
1976
1977     // Merge Z axis
1978     if( MirrorZ ){ 
1979       Int_t centralBin = DeltaY->GetAxis(2)->FindFixBin(0.0);
1980       xMin[2] = DeltaY->GetAxis(2)->GetBinLowEdge(centralBin)-0.0;
1981       nBinsNew[2] = nBinsNew[2]-centralBin+1;
1982     }
1983
1984     // Merge Angle axis
1985     if( MirrorAngle ){ 
1986       Int_t centralBin = DeltaY->GetAxis(5)->FindFixBin(0.0);
1987       xMin[5] = DeltaY->GetAxis(5)->GetBinLowEdge(centralBin)-0.0;
1988       nBinsNew[5] = nBinsNew[5]-centralBin+1;
1989     }
1990
1991     // Merge Sparse 
1992     
1993     mergedDeltaY = new THnSparseD("mergedDeltaY", "mergedDeltaY",nDim, nBinsNew, xMin, xMax );
1994     
1995     if( !mergedDeltaY ){
1996       cout << "AliTPCcalibTracks: CreateWaveCorrection: Can not copy a Sparse"<<endl;
1997       ret = -1;
1998     }
1999   }
2000   
2001   if( ret == 0 ){
2002     
2003     for( int i=0; i<nDim; i++){
2004       const TAxis *ax = DeltaY->GetAxis(i);
2005       mergedDeltaY->GetAxis(i)->SetName(ax->GetName());
2006       mergedDeltaY->GetAxis(i)->SetTitle(ax->GetTitle());
2007       if( ax->GetXbins()->fN>0 ){
2008         mergedDeltaY->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
2009       }
2010     }
2011
2012     for( Long64_t binS=0; binS<DeltaY->GetNbins(); binS++){   
2013       Double_t stat = DeltaY->GetBinContent(binS,idx); 
2014       if( stat < 1 ) continue;
2015       bool swap0=0;
2016
2017       if( MirrorPad && idx[4]>0){ // underflow reserved for contains one-pad clusters, don't mirror
2018         Double_t v = DeltaY->GetAxis(4)->GetBinCenter(idx[4]);
2019         if( v < 0.5 ) swap0 = !swap0;
2020         idx[4] = mergedDeltaY->GetAxis(4)->FindFixBin( 0.5 + TMath::Abs(0.5 - v) );
2021       }
2022       
2023       if( MirrorZ ){
2024         Double_t v = DeltaY->GetAxis(2)->GetBinCenter(idx[2]);
2025         if( v < 0.0 ) swap0 = !swap0;
2026         if( idx[2]<=0 ) idx[2] = nBinsNew[2]+1;
2027         else idx[2] = mergedDeltaY->GetAxis(2)->FindFixBin( TMath::Abs(v) );
2028       }
2029       
2030       if( MirrorAngle ){   
2031         Double_t v = DeltaY->GetAxis(5)->GetBinCenter(idx[5]);
2032         if( idx[5]<=0 ) idx[5] = nBinsNew[5]+1;
2033         else idx[5] = mergedDeltaY->GetAxis(5)->FindFixBin( TMath::Abs(v) );
2034       }
2035    
2036       if( swap0 ){
2037         if( idx[0]<=0 ) idx[0] = nBinsNew[0]+1;
2038         else if( idx[0] >= nBins[0]+1 ) idx[0] = 0;
2039         else {
2040           Double_t v = DeltaY->GetAxis(0)->GetBinCenter(idx[0]); 
2041           idx[0] = mergedDeltaY->GetAxis(0)->FindFixBin(-v);
2042         }
2043       }
2044       
2045       Long64_t bin = mergedDeltaY->GetBin(idx,kTRUE);
2046       if( bin<0 ){
2047         cout << "AliTPCcalibTracks: CreateWaveCorrection : wrong bining"<<endl;
2048         continue;
2049       }
2050       mergedDeltaY->AddBinContent(bin,stat);
2051     }
2052
2053     ret = GetTHnStat(mergedDeltaY, MeanY, SigmaY, EntrY );
2054   }
2055
2056   if( ret==0 ){
2057     
2058     MeanY->SetName("TPCWaveCorrectionMap");
2059     MeanY->SetTitle("TPCWaveCorrectionMap");
2060     SigmaY->SetName("TPCResolutionMap");
2061     SigmaY->SetTitle("TPCResolutionMap");
2062     EntrY->SetName("TPCWaveCorrectionEntr");
2063     EntrY->SetTitle("TPCWaveCorrectionEntr");
2064  
2065     for( Long64_t bin=0; bin<MeanY->GetNbins(); bin++){
2066       Double_t stat = EntrY->GetBinContent(bin,idx);
2067       
2068       // Normalize : Set no correction for one pad clusters
2069       
2070       if( idx[3]<=0 ) MeanY->SetBinContent(bin,0);
2071       
2072       // Suppress bins with low statistic
2073       
2074       if( stat<MinStat ){
2075         EntrY->SetBinContent(bin,0);
2076         MeanY->SetBinContent(bin,0);
2077         SigmaY->SetBinContent(bin,-1);
2078       }
2079     }
2080     
2081   }
2082
2083   delete[] nBins;
2084   delete[] nBinsNew;
2085   delete[] xMin;
2086   delete[] xMax;
2087   delete[] idx;
2088   delete mergedDeltaY;
2089
2090   if( ret!=0 ){
2091     delete MeanY;
2092     delete SigmaY;
2093     delete EntrY;
2094     MeanY = 0;
2095     SigmaY = 0;
2096     EntrY = 0;
2097   }
2098  
2099   return ret;
2100 }
2101
2102 int AliTPCcalibTracks::UpdateClusterParam( AliTPCClusterParam *cParam, Bool_t MirrorZ, Bool_t MirrorPad, Bool_t MirrorAngle, Int_t MinStat )
2103 {
2104   if( !cParam || !fHisDeltaY) return -1;
2105
2106   THnBase *meanY  = 0; 
2107   THnBase *sigmaY = 0;
2108   THnBase *entrY  = 0; 
2109   int ret = CreateWaveCorrection(fHisDeltaY, meanY, sigmaY, entrY,  MirrorZ, MirrorAngle, MirrorPad, MinStat  );
2110   if( ret<0 ) return ret;
2111   cParam->SetWaveCorrectionMap(meanY);
2112   cParam->SetResolutionYMap(sigmaY);
2113   delete meanY;
2114   delete sigmaY;
2115   delete entrY;
2116   return 0;
2117 }