Some of the coding violations corrected
[u/mrichter/AliRoot.git] / TPC / 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 to analyse tracks for calibration                               //
20 //     to be used as a component in AliTPCSelectorTracks                     //
21 //     In the constructor you have to specify name and title                 //
22 //     to get the Object out of a file.                                      //
23 //     The parameter 'clusterParam', a AliTPCClusterParam object             //
24 //      (needed for TPC cluster error and shape parameterization)            //
25 //     Normally you get this object out of the file 'TPCClusterParam.root'   //
26 //     In the parameter 'cuts' the cuts are specified, that decide           //
27 //     weather a track will be accepted for calibration or not.              //
28 //                                                                           //
29 //       
30 //     The data flow:
31 //     
32 /*
33    Raw Data -> Local Reconstruction -> Tracking ->     Calibration -> RefData (component itself)
34                Offline/HLT             Offline/HLT                    OCDB entries (AliTPCClusterParam) 
35 */            
36
37
38                                                                //
39 ///////////////////////////////////////////////////////////////////////////////
40
41 //
42 // ROOT includes 
43 //
44 #include <iostream>
45 #include <fstream>
46 using namespace std;
47
48 #include <TROOT.h>
49 #include <TChain.h>
50 #include <TFile.h>
51 #include <TH3F.h>
52 #include <TProfile.h>
53
54 //
55 //#include <TPDGCode.h>
56 #include <TStyle.h>
57 #include "TLinearFitter.h"
58 //#include "TMatrixD.h"
59 #include "TTreeStream.h"
60 #include "TF1.h"
61 #include <TCanvas.h>
62 #include <TGraph2DErrors.h>
63 #include "TPostScript.h"
64 #include "TCint.h"
65
66 #include <TH2D.h>
67 #include <TF2.h>
68 #include <TSystem.h>
69 #include <TCollection.h>
70 #include <iostream>
71 #include <TLinearFitter.h>
72 #include <TString.h>
73
74 //
75 // AliROOT includes 
76 //
77 #include "AliMagF.h"
78 #include "AliTracker.h"
79 #include "AliESD.h"
80 #include "AliESDtrack.h"
81 #include "AliESDfriend.h"
82 #include "AliESDfriendTrack.h" 
83 #include "AliTPCseed.h"
84 #include "AliTPCclusterMI.h"
85 #include "AliTPCROC.h"
86
87 #include "AliTPCParamSR.h"
88 #include "AliTrackPointArray.h"
89 #include "AliTPCcalibTracks.h"
90 #include "AliTPCClusterParam.h"
91 #include "AliTPCcalibTracksCuts.h"
92 #include "AliTPCCalPad.h"
93 #include "AliTPCCalROC.h"
94 #include "TText.h"
95 #include "TPaveText.h"
96 #include "TSystem.h"
97 #include "TStatToolkit.h"
98 #include "TCut.h"
99 #include "THnSparse.h"
100 #include "AliRieman.h"
101
102
103
104 ClassImp(AliTPCcalibTracks)
105
106
107 AliTPCcalibTracks::AliTPCcalibTracks():
108   AliTPCcalibBase(),
109   fClusterParam(0),
110   fROC(0),
111   fHisDeltaY(0),    // THnSparse - delta Y 
112   fHisDeltaZ(0),    // THnSparse - delta Z 
113   fHisRMSY(0),      // THnSparse - rms Y 
114   fHisRMSZ(0),      // THnSparse - rms Z 
115   fHisQmax(0),      // THnSparse - qmax 
116   fHisQtot(0),      // THnSparse - qtot 
117   fArrayQDY(0), 
118   fArrayQDZ(0), 
119   fArrayQRMSY(0),
120   fArrayQRMSZ(0),
121   fResolY(0),
122   fResolZ(0),
123   fRMSY(0),
124   fRMSZ(0),
125   fCuts(0),
126   fRejectedTracksHisto(0),
127   fClusterCutHisto(0),
128   fCalPadClusterPerPad(0),
129   fCalPadClusterPerPadRaw(0)
130
131    // 
132    // AliTPCcalibTracks default constructor
133    //    
134   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl;  
135 }   
136
137
138
139 AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
140   AliTPCcalibBase(calibTracks),
141   fClusterParam(0),
142   fROC(0),
143   fHisDeltaY(0),    // THnSparse - delta Y 
144   fHisDeltaZ(0),    // THnSparse - delta Z 
145   fHisRMSY(0),      // THnSparse - rms Y 
146   fHisRMSZ(0),      // THnSparse - rms Z 
147   fHisQmax(0),      // THnSparse - qmax 
148   fHisQtot(0),      // THnSparse - qtot 
149   fArrayQDY(0), 
150   fArrayQDZ(0), 
151   fArrayQRMSY(0),
152   fArrayQRMSZ(0),
153   fResolY(0),
154   fResolZ(0),
155   fRMSY(0),
156   fRMSZ(0),
157   fCuts(0),
158   fRejectedTracksHisto(0),
159   fClusterCutHisto(0),
160   fCalPadClusterPerPad(0),
161   fCalPadClusterPerPadRaw(0)
162 {
163    // 
164    // AliTPCcalibTracks copy constructor
165    // 
166   if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl;
167    
168    Bool_t dirStatus = TH1::AddDirectoryStatus();
169    TH1::AddDirectory(kFALSE);
170    
171    Int_t length = -1;
172    
173    (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
174    fArrayQDY= new TObjArray(length);
175    fArrayQDZ= new TObjArray(length);
176    fArrayQRMSY= new TObjArray(length);
177    fArrayQRMSZ= new TObjArray(length);
178    for (Int_t i = 0; i < length; i++) {
179       fArrayQDY->AddAt( ((TH1F*)calibTracks.fArrayQDY->At(i)->Clone()), i);
180       fArrayQDZ->AddAt( ((TH1F*)calibTracks.fArrayQDZ->At(i)->Clone()), i);
181       fArrayQRMSY->AddAt( ((TH1F*)calibTracks.fArrayQRMSY->At(i)->Clone()), i);
182       fArrayQRMSZ->AddAt( ((TH1F*)calibTracks.fArrayQRMSZ->At(i)->Clone()), i);
183    }
184    
185    (calibTracks.fResolY) ? length = calibTracks.fResolY->GetEntriesFast() : length = -1;
186    fResolY = new TObjArray(length);
187    fResolZ = new TObjArray(length);
188    fRMSY = new TObjArray(length);
189    fRMSZ = new TObjArray(length);
190    for (Int_t i = 0; i < length; i++) {
191       fResolY->AddAt( ((TH1F*)calibTracks.fResolY->At(i)->Clone()), i);
192       fResolZ->AddAt( ((TH1F*)calibTracks.fResolZ->At(i)->Clone()), i);
193       fRMSY->AddAt( ((TH1F*)calibTracks.fRMSY->At(i)->Clone()), i);
194       fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
195    } 
196    
197    
198    fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
199    fRejectedTracksHisto    = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
200    fCalPadClusterPerPad    = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
201    fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
202
203    fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(), 
204       calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise());
205    SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle());
206    TH1::AddDirectory(dirStatus); // set status back to original status
207 //    cout << "+++++ end of copy constructor +++++" << endl;   // TO BE REMOVED
208 }
209
210
211 AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibTracks){
212   //
213   // assgnment operator
214   //
215   if (this != &calibTracks) {
216     new (this) AliTPCcalibTracks(calibTracks);
217   }
218   return *this;
219
220 }
221
222
223 AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam,  AliTPCcalibTracksCuts* cuts, Int_t logLevel) : 
224   AliTPCcalibBase(),
225   fClusterParam(0),
226   fROC(0),
227   fHisDeltaY(0),    // THnSparse - delta Y 
228   fHisDeltaZ(0),    // THnSparse - delta Z 
229   fHisRMSY(0),      // THnSparse - rms Y 
230   fHisRMSZ(0),      // THnSparse - rms Z 
231   fHisQmax(0),      // THnSparse - qmax 
232   fHisQtot(0),      // THnSparse - qtot 
233   fArrayQDY(0), 
234   fArrayQDZ(0), 
235   fArrayQRMSY(0),
236   fArrayQRMSZ(0),
237   fResolY(0),
238   fResolZ(0),
239   fRMSY(0),
240   fRMSZ(0),
241   fCuts(0),
242   fRejectedTracksHisto(0),
243   fClusterCutHisto(0),
244   fCalPadClusterPerPad(0),
245   fCalPadClusterPerPadRaw(0)
246  {
247    // 
248    // AliTPCcalibTracks constructor
249    // specify 'name' and 'title' of your object
250    // specify 'clusterParam', (needed for TPC cluster error and shape parameterization)
251    // In the parameter 'cuts' the cuts are specified, that decide           
252    // weather a track will be accepted for calibration or not.              
253    //
254    // fDebugLevel - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStreamer, 6: waste your screen
255    // 
256    // All histograms are instatiated in this constructor.
257    // 
258    this->SetName(name);
259    this->SetTitle(title);
260
261    if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl;
262    G__SetCatchException(0);     
263    
264    fClusterParam = clusterParam;
265    if (fClusterParam){
266      fClusterParam->SetInstance(fClusterParam);
267    }
268    else {
269      Error("AliTPCcalibTracks","No cluster parametrization found! A valid clusterParam object is needed in the constructor. (To be found in 'TPCClusterParam.root'.)");
270    } 
271    fCuts = cuts;
272    SetDebugLevel(logLevel);
273    MakeHistos();
274    
275    TH1::AddDirectory(kFALSE);
276    
277    fRejectedTracksHisto    = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
278    fCalPadClusterPerPad    = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
279    fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
280    fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
281    
282    
283    TH1::AddDirectory(kFALSE);
284    
285    
286    fResolY = new TObjArray(3);
287    fResolZ = new TObjArray(3);
288    fRMSY   = new TObjArray(3);
289    fRMSZ   = new TObjArray(3);
290    TH3F * his3D;
291    //
292    his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
293    fResolY->AddAt(his3D,0);     
294    his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
295    fResolY->AddAt(his3D,1);
296    his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
297    fResolY->AddAt(his3D,2);
298    //
299    his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
300    fResolZ->AddAt(his3D,0);
301    his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
302    fResolZ->AddAt(his3D,1);
303    his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
304    fResolZ->AddAt(his3D,2);
305    //
306    his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
307    fRMSY->AddAt(his3D,0);
308    his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
309    fRMSY->AddAt(his3D,1);
310    his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
311    fRMSY->AddAt(his3D,2);
312    //
313    his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
314    fRMSZ->AddAt(his3D,0);
315    his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
316    fRMSZ->AddAt(his3D,1);
317    his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
318    fRMSZ->AddAt(his3D,2);
319    //
320       
321    TH1::AddDirectory(kFALSE);
322    
323    fArrayQDY = new TObjArray(300);
324    fArrayQDZ = new TObjArray(300);
325    fArrayQRMSY = new TObjArray(300);
326    fArrayQRMSZ = new TObjArray(300);
327    for (Int_t iq = 0; iq <= 10; iq++){
328       for (Int_t ipad = 0; ipad < 3; ipad++){
329          Int_t   bin   = GetBin(iq, ipad);
330          Float_t qmean = GetQ(bin);
331          char hname[200];
332          snprintf(hname,100,"ResolY Pad%d Qmiddle%f",ipad, qmean);
333          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
334          fArrayQDY->AddAt(his3D, bin);
335          snprintf(hname,100,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
336          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
337          fArrayQDZ->AddAt(his3D, bin);
338          snprintf(hname,100,"RMSY Pad%d Qmiddle%f",ipad, qmean);
339          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
340          fArrayQRMSY->AddAt(his3D, bin);
341          snprintf(hname,100,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
342          his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
343          fArrayQRMSZ->AddAt(his3D, bin);
344       }
345    }
346    
347    
348
349    if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl; 
350    cout << "end of main constructor" << endl; // TO BE REMOVED
351 }    
352
353
354 AliTPCcalibTracks::~AliTPCcalibTracks() {
355    // 
356    // AliTPCcalibTracks destructor
357    // 
358    
359   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
360    Int_t length = 0;
361    
362    
363    if (fResolY) {
364      length = fResolY->GetEntriesFast();
365      for (Int_t i = 0; i < length; i++){
366        delete fResolY->At(i);
367        delete fResolZ->At(i);
368        delete fRMSY->At(i);
369        delete fRMSZ->At(i);
370      }
371      delete fResolY;
372      delete fResolZ;
373      delete fRMSY;
374      delete fRMSZ;
375    }
376    
377    if (fArrayQDY) {
378      length = fArrayQDY->GetEntriesFast();
379      for (Int_t i = 0; i < length; i++){
380        delete fArrayQDY->At(i);
381        delete fArrayQDZ->At(i);
382        delete fArrayQRMSY->At(i);
383        delete fArrayQRMSZ->At(i);
384      }
385    }
386    
387    
388     
389    delete fArrayQDY;
390    delete fArrayQDZ;
391    delete fArrayQRMSY;
392    delete fArrayQRMSZ;
393    
394   delete fRejectedTracksHisto;
395   delete fClusterCutHisto;
396   if (fCalPadClusterPerPad)    delete fCalPadClusterPerPad;
397   if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
398   delete fHisDeltaY;    // THnSparse - delta Y 
399   delete fHisDeltaZ;    // THnSparse - delta Z 
400   delete fHisRMSY;      // THnSparse - rms Y 
401   delete fHisRMSZ;      // THnSparse - rms Z 
402   delete fHisQmax;      // THnSparse - qmax 
403   delete fHisQtot;      // THnSparse - qtot 
404
405 }
406    
407   
408
409 void AliTPCcalibTracks::Process(AliTPCseed *track){
410    // 
411    // To be called in the selector
412    // first AcceptTrack is evaluated, then calls all the following analyse functions: 
413    // FillResolutionHistoLocal(track)
414
415    // 
416   if (GetDebugLevel() > 5) Info("Process","Starting to process the track...");
417    Int_t accpetStatus = AcceptTrack(track);
418    if (accpetStatus == 0) {
419       FillResolutionHistoLocal(track);
420    }
421    else fRejectedTracksHisto->Fill(accpetStatus);
422 }
423
424
425
426 Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){
427   //
428   // calculate bins for given q and pad type 
429   // used in TObjArray
430   //
431   Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 );  
432   res *= 3;
433   res += pad;
434   return res;
435 }
436
437
438 Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){
439   //
440   // calculate bins for given iq and pad type 
441   // used in TObjArray
442   //
443   return iq * 3 + pad;;
444 }
445
446
447 Float_t AliTPCcalibTracks::GetQ(Int_t bin){
448    // 
449    // returns to bin belonging charge
450    // (bin / 3 + 3)^2
451    // 
452    Int_t bin0 = bin / 3;
453    bin0 += 3;
454    return bin0 * bin0;
455 }
456
457
458 Float_t AliTPCcalibTracks::GetPad(Int_t bin){
459    // 
460    // returns to bin belonging pad
461    // bin % 3
462    // 
463    return bin % 3; 
464 }
465
466
467
468 Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
469   //
470   // Function, that decides wheather a given track is accepted for 
471   // the analysis or not. 
472   // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts'
473   // Returns 0 if a track is accepted or an integer different from 0 
474   // to indicate the failed cut
475   //
476   const Int_t   kMinClusters  = fCuts->GetMinClusters();
477   const Float_t kMinRatio     = fCuts->GetMinRatio();
478   const Float_t kMax1pt       = fCuts->GetMax1pt();
479   const Float_t kEdgeYXCutNoise    = fCuts->GetEdgeYXCutNoise();
480   const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise();
481   
482   //
483   // edge induced noise tracks - NEXT RELEASE will be removed during tracking
484   if ( TMath::Abs(track->GetY() / track->GetX()) > kEdgeYXCutNoise )
485     if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return 1;
486   if (track->GetNumberOfClusters() < kMinClusters) return 2;
487   Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.);
488   if (ratio < kMinRatio) return 3;
489 //   Float_t mpt = track->Get1Pt();       // Get1Pt() doesn't exist any more
490   Float_t mpt = track->GetSigned1Pt();
491   if (TMath::Abs(mpt) > kMax1pt) return 4;
492   //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
493   //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
494   //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
495   
496   if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");  
497   return 0;
498 }
499
500
501 void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
502    //
503    // fill resolution histograms - localy - tracklet in the neighborhood
504    // write debug information to 'TPCSelectorDebug.root'
505    // 
506    // _ the main function, called during track analysis _
507    // 
508    // loop over all padrows along the track
509    // fit tracklets (length: 13 clusters) calculate mean chi^2 for this track-fit in Y and Z direction
510    // 
511    // loop again over all padrows along the track
512    // fit tracklet (clusters in given padrow +- kDelta padrows) 
513    // with polynom of 2nd order and two polynoms of 1st order
514    // take both polynoms of 1st order, calculate difference of their parameters
515    // add covariance matrixes and calculate chi2 of this difference
516    // if this chi2 is bigger than a given threshold, assume that the current cluster is
517    // a kink an goto next padrow
518    // if not:
519    // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
520    // 
521    // write debug information to 'TPCSelectorDebug.root'
522    // only for every kDeltaWriteDebugStream'th padrow to reduce data volume 
523    // and to avoid redundant data
524    // 
525
526   static TLinearFitter fFitterParY(3,"pol2");    // 
527   static TLinearFitter fFitterParZ(3,"pol2");    //
528   static TLinearFitter fFitterParYWeight(3,"pol2");    // 
529   static TLinearFitter fFitterParZWeight(3,"pol2");    //
530   fFitterParY.StoreData(kTRUE);
531   fFitterParZ.StoreData(kTRUE);
532   fFitterParYWeight.StoreData(kTRUE);
533   fFitterParZWeight.StoreData(kTRUE);
534   if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
535   const Int_t   kDelta     = 10;          // delta rows to fit
536   const Float_t kMinRatio  = 0.75;        // minimal ratio
537   const Float_t kChi2Cut   =  10.;          // cut chi2 - left right
538   const Float_t kSigmaCut  = 3.;        //sigma cut
539   const Float_t kdEdxCut=300;
540   const Float_t kPtCut=0.300;
541
542   if (track->GetTPCsignal()>kdEdxCut) return;  
543   if (TMath::Abs(AliTracker::GetBz())>0.1 &&TMath::Abs(track->Pt())<kPtCut) return;  
544
545   // estimate mean error
546   Int_t nTrackletsAll = 0;       // number of tracklets for given track
547   Float_t csigmaY     = 0;       // mean sigma for tracklet refit in Y direction
548   Float_t csigmaZ     = 0;       // mean sigma for tracklet refit in Z direction
549   Int_t nClusters     = 0;       // working variable, number of clusters per tracklet
550   Int_t sectorG       = -1;      // working variable, sector of tracklet, has to stay constant for one tracklet
551   Double_t refX=0;
552   // ---------------------------------------------------------------------
553   for (Int_t irow = 0; irow < 159; irow++){
554     // loop over all rows along the track
555     // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
556     // calculate mean chi^2 for this track-fit in Y and Z direction
557     AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
558     if (!cluster0) continue;  // no cluster found
559     Int_t sector = cluster0->GetDetector();
560     
561     Int_t ipad = TMath::Nint(cluster0->GetPad());
562     Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
563     fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
564     
565     if (sector != sectorG){
566       // track leaves sector before it crossed enough rows to fit / initialization
567       nClusters = 0;
568       fFitterParY.ClearPoints();
569       fFitterParZ.ClearPoints();
570       sectorG = sector;
571     }
572     else {
573       nClusters++;
574       if (refX<1) refX=cluster0->GetX()+kDelta*0.5;
575       Double_t x = cluster0->GetX()-refX;
576       fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
577       fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
578       //
579       if ( nClusters >= kDelta + 3 ){  
580         // if more than 13 (kDelta+3) clusters were added to the fitters
581         // fit the tracklet, increase trackletCounter
582         fFitterParY.Eval();
583         fFitterParZ.Eval();
584         nTrackletsAll++;
585         csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
586         csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
587         nClusters = -1;
588         fFitterParY.ClearPoints();
589         fFitterParZ.ClearPoints();
590         refX=0;
591       }
592     }
593   }      // for (Int_t irow = 0; irow < 159; irow++)
594   // mean chi^2 for all tracklet fits in Y and in Z direction: 
595   csigmaY = TMath::Sqrt(TMath::Abs(csigmaY) / (nTrackletsAll+0.1));
596   csigmaZ = TMath::Sqrt(TMath::Abs(csigmaZ) / (nTrackletsAll+0.1));
597   // ---------------------------------------------------------------------
598   //
599   //
600
601   for (Int_t irow = kDelta; irow < 159-kDelta; irow++){
602     // loop again over all rows along the track
603     // do analysis
604     // 
605     Int_t nclFound = 0;  // number of clusters in the neighborhood
606     Int_t ncl0 = 0;      // number of clusters in rows < rowOfCenterCluster
607     Int_t ncl1 = 0;      // number of clusters in rows > rowOfCenterCluster
608     AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
609     if (!cluster0) continue;
610     Int_t sector = cluster0->GetDetector();
611     Float_t xref = cluster0->GetX();
612     
613     // Make Fit
614     fFitterParY.ClearPoints();
615     fFitterParZ.ClearPoints();    
616     fFitterParYWeight.ClearPoints();
617     fFitterParZWeight.ClearPoints();    
618     // fit tracklet (clusters in given padrow +- kDelta padrows) 
619     // with polynom of 2nd order and two polynoms of 1st order
620     // take both polynoms of 1st order, calculate difference of their parameters
621     // add covariance matrixes and calculate chi2 of this difference
622     // if this chi2 is bigger than a given threshold, assume that the current cluster is
623     // a kink an goto next padrow    
624     AliRieman riemanFit(2*kDelta);
625     AliRieman riemanFitW(2*kDelta);
626     for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
627       // loop over irow +- kDelta rows (neighboured rows)
628       // 
629       // 
630       if (idelta == 0) continue;                                // don't use center cluster
631       if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
632       AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
633       if (!currentCluster) continue;
634       if (currentCluster->GetType() < 0) continue;
635       if (currentCluster->GetDetector() != sector) continue;
636       nclFound++;
637       if (idelta < 0){
638         ncl0++;
639       }
640       if (idelta > 0){
641         ncl1++;
642       }
643       riemanFit.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY,csigmaZ);
644       riemanFitW.AddPoint(currentCluster->GetX(), currentCluster->GetY(),currentCluster->GetZ(), csigmaY*TMath::Sqrt(1+TMath::Abs(idelta)),csigmaZ*TMath::Sqrt(1+TMath::Abs(idelta)));
645     }  // loop over neighbourhood for fitter filling 
646     if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
647     if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
648     if (nclFound < kDelta * kMinRatio) continue;    // if not enough clusters (7.5) found in neighbourhood goto next padrow
649     riemanFit.Update();
650     riemanFitW.Update();
651     Double_t chi2R=TMath::Sqrt(riemanFit.CalcChi2()/(2*nclFound-5));
652     Double_t chi2RW=TMath::Sqrt(riemanFitW.CalcChi2()/(2*nclFound-5));
653     Double_t paramYR[3], paramZR[3];
654     Double_t paramYRW[3], paramZRW[3];
655     //
656     paramYR[0]    = riemanFit.GetYat(xref);
657     paramZR[0]    = riemanFit.GetZat(xref);
658     paramYRW[0]   = riemanFitW.GetYat(xref);
659     paramZRW[0]   = riemanFitW.GetZat(xref);
660     //
661     paramYR[1]    = riemanFit.GetDYat(xref);
662     paramZR[1]    = riemanFit.GetDZat(xref);
663     paramYRW[1]   = riemanFitW.GetDYat(xref);
664     paramZRW[1]   = riemanFitW.GetDZat(xref);
665     //
666     Int_t reject=0;
667     if (chi2R>kChi2Cut) reject+=1;
668     if (chi2RW>kChi2Cut) reject+=2;
669     if (TMath::Abs(paramYR[0]-paramYRW[0])>kSigmaCut*csigmaY) reject+=4;
670     if (TMath::Abs(paramZR[0]-paramZRW[0])>kSigmaCut*csigmaZ) reject+=8;
671     if (TMath::Abs(paramYR[1]-paramYRW[1])>csigmaY) reject+=16;
672     if (TMath::Abs(paramZR[1]-paramZRW[1])>csigmaZ) reject+=32;
673     //
674     TTreeSRedirector *cstream = GetDebugStreamer();    
675     // get fit parameters from pol2 fit:     
676     Double_t tracky = paramYR[0];
677     Double_t trackz = paramZR[0];
678     Float_t  deltay = cluster0->GetY()-tracky;
679     Float_t  deltaz = cluster0->GetZ()-trackz;
680     Float_t  angley = paramYR[1];
681     Float_t  anglez = paramZR[1];
682     
683     Int_t padSize = 0;                          // short pads
684     if (cluster0->GetDetector() >= 36) {
685       padSize = 1;                              // medium pads 
686       if (cluster0->GetRow() > 63) padSize = 2; // long pads
687     }
688     Int_t ipad = TMath::Nint(cluster0->GetPad());
689     if (cstream){
690       Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
691       (*cstream)<<"Resol0"<<    
692         "run="<<fRun<<              //  run number
693         "event="<<fEvent<<          //  event number
694         "time="<<fTime<<            //  time stamp of event
695         "trigger="<<fTrigger<<      //  trigger
696         "mag="<<fMagF<<             //  magnetic field        
697         "padSize="<<padSize<<
698         //
699         "reject="<<reject<<
700         "cl.="<<cluster0<<          // cluster info
701         "xref="<<xref<<             // cluster reference X
702         //rieman fit
703         "yr="<<paramYR[0]<<         // track position y - rieman fit
704         "zr="<<paramZR[0]<<         // track position z - rieman fit 
705         "yrW="<<paramYRW[0]<<         // track position y - rieman fit - weight
706         "zrW="<<paramZRW[0]<<         // track position z - rieman fit - weight
707         "dyr="<<paramYR[1]<<         // track position y - rieman fit
708         "dzr="<<paramZR[1]<<         // track position z - rieman fit 
709         "dyrW="<<paramYRW[1]<<         // track position y - rieman fit - weight
710         "dzrW="<<paramZRW[1]<<         // track position z - rieman fit - weight
711         //
712         "angley="<<angley<<         // angle par fit
713         "anglez="<<anglez<<         // angle par fit
714         "zdr="<<zdrift<<            //
715         "dy="<<deltay<<
716         "dz="<<deltaz<<        
717         "sy="<<csigmaY<<
718         "sz="<<csigmaZ<<
719         "chi2R="<<chi2R<<
720         "chi2RW="<<chi2RW<<
721         "\n";
722     }    
723     if (reject) continue;
724
725     
726     // =========================================
727     // wirte collected information to histograms
728     // =========================================
729         
730     Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
731     fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
732     
733     
734     TH3F * his3 = 0;
735     his3 = (TH3F*)fRMSY->At(padSize);
736     if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
737     his3 = (TH3F*)fRMSZ->At(padSize);
738     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
739     
740     his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
741     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
742     his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
743     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
744     
745     
746     his3 = (TH3F*)fResolY->At(padSize);
747     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
748     his3 = (TH3F*)fResolZ->At(padSize);
749     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
750     his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
751     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
752     his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
753     if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );        
754     //=============================================================================================    
755     //
756     // Fill THN histograms
757     //
758     Double_t xvar[9];
759     xvar[1]=padSize;
760     xvar[2]=cluster0->GetZ();
761     xvar[3]=cluster0->GetMax();
762     
763     xvar[0]=deltay;
764     xvar[4]=cluster0->GetPad()-Int_t(cluster0->GetPad());      
765     xvar[5]=angley;
766     fHisDeltaY->Fill(xvar);
767     xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
768     fHisRMSY->Fill(xvar);
769     
770     xvar[0]=deltaz;
771     xvar[4]=cluster0->GetTimeBin()-Int_t(cluster0->GetTimeBin());
772     xvar[5]=anglez;
773     fHisDeltaZ->Fill(xvar);
774     xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
775     fHisRMSZ->Fill(xvar);
776     
777   }    // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
778 }  // FillResolutionHistoLocal(...)
779
780
781
782
783
784
785
786
787 void  AliTPCcalibTracks::SetStyle() const {
788    // 
789    // set style, can be called by all draw functions
790    // 
791    gROOT->SetStyle("Plain");
792    gStyle->SetFillColor(10);
793    gStyle->SetPadColor(10);
794    gStyle->SetCanvasColor(10);
795    gStyle->SetStatColor(10);
796    gStyle->SetPalette(1,0);
797    gStyle->SetNumberContours(60);
798 }
799
800
801
802 void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){ 
803    // 
804    // all functions are called, that produce pictures
805    // the histograms are written to the directory 'pathName'
806    // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
807    // 'stat' is also the number of minEntries for MakeResPlotsQTree
808    // 
809
810   if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
811   MakeResPlotsQTree(stat, pathName);
812 }
813    
814
815
816
817 void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
818    //
819    // Make tree with resolution parameters
820    // the result is written to 'resol.root' in directory 'pathname'
821    // file information are available in fileInfo
822    // available variables in the tree 'Resol':
823    //  Entries: number of entries for this resolution point
824    //  nbins:   number of bins that were accumulated
825    //  Dim:     direction, Dim==0: y-direction, Dim==1: z-direction
826    //  Pad:     padSize; short, medium and long
827    //  Length:  pad length, 0.75, 1, 1.5
828    //  QMean:   mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
829    //  Zc:      center of middle bin in drift direction
830    //  Zm:      mean dirftlength for accumulated Delta-Histograms
831    //  Zs:      width of driftlength bin
832    //  AngleC:  center of middle bin in Angle-Direction
833    //  AngleM:  mean angle for accumulated Delta-Histograms
834    //  AngleS:  width of Angle-bin
835    //  Resol:   sigma for gaus fit through Delta-Histograms
836    //  Sigma:   error of sigma for gaus fit through Delta Histograms
837    //  MeanR:   mean of the Delta-Histogram
838    //  SigmaR:  rms of the Delta-Histogram
839    //  RMSm:    mean of the gaus fit through RMS-Histogram
840    //  RMS:     sigma of the gaus fit through RMS-Histogram
841    //  RMSe0:   error of mean of gaus fit in RMS-Histogram
842    //  RMSe1:   error of sigma of gaus fit in RMS-Histogram
843    //  
844       
845   if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
846   if (GetDebugLevel() > -1) cout << "    relax, the calculation will take a while..." << endl;
847   
848    gSystem->MakeDirectory(pathName);
849    gSystem->ChangeDirectory(pathName);
850    TString kFileName = "resol.root";
851    TTreeSRedirector fTreeResol(kFileName.Data());
852    
853    TH3F *resArray[2][3][11];
854    TH3F *rmsArray[2][3][11];
855   
856    // load histograms from fArrayQDY and fArrayQDZ 
857    // into resArray and rmsArray
858    // that is all we need here
859    for (Int_t idim = 0; idim < 2; idim++){
860       for (Int_t ipad = 0; ipad < 3; ipad++){
861          for (Int_t iq = 0; iq <= 10; iq++){
862             rmsArray[idim][ipad][iq]=0;
863             resArray[idim][ipad][iq]=0;
864             Int_t bin = GetBin(iq,ipad); 
865             TH3F *hresl = 0;
866             if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
867             if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
868             if (!hresl) continue;
869             resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
870             resArray[idim][ipad][iq]->SetDirectory(0);
871             TH3F * hreslRMS = 0;
872             if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
873             if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
874             if (!hreslRMS) continue;
875             rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
876             rmsArray[idim][ipad][iq]->SetDirectory(0);
877          }
878       }
879    }
880    if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
881    
882    //--------------------------------------------------------------------------------------------
883    
884    char name[200];
885    Double_t qMean;
886    Double_t zMean, angleMean, zCenter, angleCenter;
887    Double_t zSigma, angleSigma;
888    TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
889    TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
890    TF1 *fitFunction = new TF1("fitFunction", "gaus");
891    Float_t entriesQ = 0;
892    Int_t loopCounter = 1;
893   
894    for (Int_t idim = 0; idim < 2; idim++){
895       // Loop y-z corrdinate
896       for (Int_t ipad = 0; ipad < 3; ipad++){
897          // loop pad type
898          for (Int_t iq = -1; iq < 10; iq++){
899             // LOOP Q
900            if (GetDebugLevel() > -1) 
901                cout << "Loop-counter, this is loop " << loopCounter << " of 66, (" 
902                   << (Int_t)((loopCounter)/66.*100) << "% done), " 
903                   << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << "  \r" << std::flush;
904             loopCounter++;
905             TH3F *hres = 0;
906             TH3F *hrms = 0;
907             qMean    = 0;
908             entriesQ = 0;
909             
910             // calculate qMean
911             if (iq == -1){
912                // integrated spectra
913                for (Int_t iql = 0; iql < 10; iql++){    
914                   Int_t bin = GetBin(iql,ipad); 
915                   TH3F *hresl = resArray[idim][ipad][iql];
916                   TH3F *hrmsl = rmsArray[idim][ipad][iql];
917                   if (!hresl) continue;
918                   if (!hrmsl) continue;     
919                   entriesQ += hresl->GetEntries();
920                   qMean += hresl->GetEntries() * GetQ(bin);      
921                   if (!hres) {
922                      hres = (TH3F*)hresl->Clone();
923                      hrms = (TH3F*)hrmsl->Clone();
924                   }
925                   else{
926                      hres->Add(hresl);
927                      hrms->Add(hrmsl);
928                   }
929                }
930                qMean /= entriesQ;
931                qMean *= -1.;  // integral mean charge
932             }
933             else {
934                // loop over neighboured Q-bins 
935                // accumulate entries from neighboured Q-bins
936                for (Int_t iql = iq - 1; iql <= iq + 1; iql++){              
937                   if (iql < 0) continue;
938                   Int_t bin = GetBin(iql,ipad);
939                   TH3F * hresl = resArray[idim][ipad][iql];
940                   TH3F * hrmsl = rmsArray[idim][ipad][iql];
941                   if (!hresl) continue;
942                   if (!hrmsl) continue;
943                   entriesQ += hresl->GetEntries(); 
944                   qMean += hresl->GetEntries() * GetQ(bin);      
945                   if (!hres) {
946                      hres = (TH3F*) hresl->Clone();
947                      hrms = (TH3F*) hrmsl->Clone();
948                   }
949                   else{
950                      hres->Add(hresl);
951                      hrms->Add(hrmsl);
952                   }
953                }
954                qMean/=entriesQ;
955             }
956             if (!hres) continue;
957             if (!hrms) continue;
958
959             TAxis *xAxisDriftLength = hres->GetXaxis();   // driftlength / z - axis
960             TAxis *yAxisAngle       = hres->GetYaxis();   // angle axis
961             TAxis *zAxisDelta       = hres->GetZaxis();   // delta axis
962             TAxis *zAxisRms         = hrms->GetZaxis();   // rms axis
963             
964             // loop over all angle bins
965             for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
966                angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
967                // loop over all driftlength bins
968                for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
969                   zCenter    = xAxisDriftLength->GetBinCenter(ibinxDL);
970                   zSigma     = xAxisDriftLength->GetBinWidth(ibinxDL);
971                   angleSigma = yAxisAngle->GetBinWidth(ibinyAngle); 
972                   zMean      = zCenter;      // changens, when more statistic is accumulated
973                   angleMean  = angleCenter;  // changens, when more statistic is accumulated
974                   
975                   // create 2 1D-Histograms, projectionRes and projectionRms
976                   // these histograms are delta histograms for given direction, padSize, chargeBin,
977                   // angleBin and driftLengthBin
978                   // later on they will be fitted with a gausian, its sigma is the resoltuion...
979                   snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
980                   // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
981                   projectionRes->SetNameTitle(name, name);
982                   snprintf(name,200,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
983                   // TH1D * projectionRms =  new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
984                   projectionRms->SetNameTitle(name, name);
985                   
986                   projectionRes->Reset();
987                   projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
988                   projectionRms->Reset();
989                   projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
990                   projectionRes->SetDirectory(0);
991                   projectionRms->SetDirectory(0);
992                   
993                   Double_t entries = 0;
994                   Int_t    nbins   = 0;   // counts, how many bins were accumulated
995                   zMean     = 0;
996                   angleMean = 0;
997                   entries   = 0;
998                   
999                   // fill projectionRes and projectionRms for given dim, ipad and iq, 
1000                   // as well as for given angleBin and driftlengthBin
1001                   // if this gives not enough statistic, include neighbourhood 
1002                   // (angle and driftlength) successifely
1003                   for (Int_t dbin = 0; dbin <= 8; dbin++){              // delta-bins around centered angleBin and driftlengthBin
1004                      for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) {   // delta-bins in angle direction
1005                         for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
1006                            if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue;   // add each bin only one time !
1007                            Int_t binx2 = ibinxDL + dbinx2;                       // position variable in x (driftlength) direction
1008                            Int_t biny2 = ibinyAngle + dbiny2;                    // position variable in y (angle)  direction
1009                            if (binx2 < 1 || biny2 < 1) continue;                 // don't go out of the histogram!
1010                            if (binx2 >= xAxisDriftLength->GetNbins()) continue;  // don't go out of the histogram!
1011                            if (biny2 >= yAxisAngle->GetNbins()) continue;        // don't go out of the histogram!
1012                            nbins++;                                              // count the number of accumulated bins
1013                            // Fill resolution histo
1014                            for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
1015                               // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3);     // unused variable
1016                               projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
1017                               entries   += hres->GetBinContent(binx2, biny2, ibin3);
1018                               zMean     += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
1019                               angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
1020                            }  // ibin3 loop
1021                            // fill RMS histo
1022                            for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
1023                               projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
1024                            }
1025                         }  //dbinx2 loop
1026                         if (entries > minEntries) break; // enough statistic accumulated
1027                      }  // dbiny2 loop
1028                      if (entries > minEntries) break;    // enough statistic accumulated
1029                   }  // dbin loop
1030                   if ( entries< minEntries) continue;  // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree  
1031                   zMean /= entries;
1032                   angleMean /= entries;
1033                   
1034                   if (entries > minEntries) {
1035                      //  when enough statistic is accumulated
1036                      //  fit Delta histograms with a gausian
1037                      //  of the gausian is the resolution (resol), its fit error is sigma
1038                      //  store also mean and RMS of the histogram
1039                      Float_t xmin     = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1040                      Float_t xmax     = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1041                      
1042 //                      projectionRes->Fit("gaus", "q0", "", xmin, xmax);
1043 //                      Float_t resol    = projectionRes->GetFunction("gaus")->GetParameter(2);
1044 //                      Float_t sigma    = projectionRes->GetFunction("gaus")->GetParError(2);
1045                      fitFunction->SetMaximum(xmax);
1046                      fitFunction->SetMinimum(xmin);
1047                      projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
1048                      Float_t resol    = fitFunction->GetParameter(2);
1049                      Float_t sigma    = fitFunction->GetParError(2);
1050                      
1051                      Float_t meanR    = projectionRes->GetMean();
1052                      Float_t sigmaR   = projectionRes->GetRMS();
1053                      // fit also RMS histograms with a gausian
1054                      // store mean and sigma of the gausian in rmsMean and rmsSigma
1055                      // store also the fit errors in errorRMS and errorSigma
1056                      xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
1057                      xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
1058                      
1059 //                      projectionRms->Fit("gaus","q0","",xmin,xmax);
1060 //                      Float_t rmsMean    = projectionRms->GetFunction("gaus")->GetParameter(1);
1061 //                      Float_t rmsSigma   = projectionRms->GetFunction("gaus")->GetParameter(2);
1062 //                      Float_t errorRMS   = projectionRms->GetFunction("gaus")->GetParError(1);
1063 //                      Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
1064                      projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
1065                      Float_t rmsMean    = fitFunction->GetParameter(1);
1066                      Float_t rmsSigma   = fitFunction->GetParameter(2);
1067                      Float_t errorRMS   = fitFunction->GetParError(1);
1068                      Float_t errorSigma = fitFunction->GetParError(2);
1069                     
1070                      Float_t length = 0.75;
1071                      if (ipad == 1) length = 1;
1072                      if (ipad == 2) length = 1.5;
1073                      
1074                      fTreeResol<<"Resol"<<
1075                         "Entries="<<entries<<      // number of entries for this resolution point
1076                         "nbins="<<nbins<<          // number of bins that were accumulated
1077                         "Dim="<<idim<<             // direction, Dim==0: y-direction, Dim==1: z-direction
1078                         "Pad="<<ipad<<             // padSize; short, medium and long
1079                         "Length="<<length<<        // pad length, 0.75, 1, 1.5
1080                         "QMean="<<qMean<<          // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1081                         "Zc="<<zCenter<<           // center of middle bin in drift direction
1082                         "Zm="<<zMean<<             // mean dirftlength for accumulated Delta-Histograms
1083                         "Zs="<<zSigma<<            // width of driftlength bin
1084                         "AngleC="<<angleCenter<<   // center of middle bin in Angle-Direction
1085                         "AngleM="<<angleMean<<     // mean angle for accumulated Delta-Histograms
1086                         "AngleS="<<angleSigma<<    // width of Angle-bin
1087                         "Resol="<<resol<<          // sigma for gaus fit through Delta-Histograms
1088                         "Sigma="<<sigma<<          // error of sigma for gaus fit through Delta Histograms
1089                         "MeanR="<<meanR<<          // mean of the Delta-Histogram
1090                         "SigmaR="<<sigmaR<<        // rms of the Delta-Histogram
1091                         "RMSm="<<rmsMean<<         // mean of the gaus fit through RMS-Histogram
1092                         "RMSs="<<rmsSigma<<        // sigma of the gaus fit through RMS-Histogram
1093                         "RMSe0="<<errorRMS<<       // error of mean of gaus fit in RMS-Histogram
1094                         "RMSe1="<<errorSigma<<     // error of sigma of gaus fit in RMS-Histogram
1095                         "\n";
1096                      if (GetDebugLevel() > 5) {
1097                         projectionRes->SetDirectory(fTreeResol.GetFile());
1098                         projectionRes->Write(projectionRes->GetName());
1099                         projectionRes->SetDirectory(0);
1100                         projectionRms->SetDirectory(fTreeResol.GetFile());
1101                         projectionRms->Write(projectionRms->GetName());
1102                         projectionRes->SetDirectory(0);
1103                      }
1104                   }  // if (projectionRes->GetSum() > minEntries)
1105                }  // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
1106             }  // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
1107             
1108          }  // iq-loop
1109       }  // ipad-loop
1110    }  // idim-loop
1111    delete projectionRes;
1112    delete projectionRms;
1113    
1114 //    TFile resolFile(fTreeResol.GetFile());
1115    TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
1116    fileInfo.Write("fileInfo");
1117 //    resolFile.Close();
1118 //    fTreeResol.GetFile()->Close();
1119    if (GetDebugLevel() > -1) cout << endl;
1120    if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
1121    gSystem->ChangeDirectory("..");
1122 }
1123
1124
1125
1126
1127
1128 Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
1129    // 
1130    // function to merge several AliTPCcalibTracks objects after PROOF calculation
1131    // The object's histograms are merged via their merge functions
1132    // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
1133    // 
1134    
1135   if (GetDebugLevel() > 0) cout << " *****  this is AliTPCcalibTracks::Merge(TCollection *collectionList)  *****"<< endl;  
1136    if (!collectionList) return 0;
1137    if (collectionList->IsEmpty()) return -1;
1138    
1139    if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl;     //    REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
1140    if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
1141    collectionList->Print();
1142    
1143    // create a list for each data member
1144    TList* deltaYList = new TList;
1145    TList* deltaZList = new TList;
1146    TList* arrayAmpRowList = new TList;
1147    TList* rejectedTracksList = new TList;
1148    TList* clusterCutHistoList = new TList;
1149    TList* arrayAmpList = new TList;
1150    TList* arrayQDYList = new TList;
1151    TList* arrayQDZList = new TList;
1152    TList* arrayQRMSYList = new TList;
1153    TList* arrayQRMSZList = new TList;
1154    TList* resolYList = new TList;
1155    TList* resolZList = new TList;
1156    TList* rMSYList = new TList;
1157    TList* rMSZList = new TList;
1158    
1159 //    TList* nRowsList = new TList;
1160 //    TList* nSectList = new TList;
1161 //    TList* fileNoList = new TList;
1162    
1163    TIterator *listIterator = collectionList->MakeIterator();
1164    AliTPCcalibTracks *calibTracks = 0;
1165    if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;    
1166    Int_t counter = 0;
1167    while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
1168       // loop over all entries in the collectionList and get dataMembers into lists
1169       
1170       arrayQDYList->Add(calibTracks->GetfArrayQDY());
1171       arrayQDZList->Add(calibTracks->GetfArrayQDZ());
1172       arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
1173       arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
1174       resolYList->Add(calibTracks->GetfResolY());
1175       resolZList->Add(calibTracks->GetfResolZ());
1176       rMSYList->Add(calibTracks->GetfRMSY());
1177       rMSZList->Add(calibTracks->GetfRMSZ());
1178       rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
1179       clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
1180       //
1181       if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
1182         fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());      
1183       //      fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
1184       counter++;
1185       if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
1186       AddHistos(calibTracks);
1187    }
1188    
1189    
1190    // merge data members
1191    if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl; 
1192    fClusterCutHisto->Merge(clusterCutHistoList);
1193    fRejectedTracksHisto->Merge(rejectedTracksList);
1194    
1195    TObjArray* objarray = 0;
1196    TH1* hist = 0;
1197    TList* histList = 0;
1198    TIterator *objListIterator = 0;
1199    
1200       
1201    if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
1202    // merge fArrayQDY
1203    for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
1204       objListIterator = arrayQDYList->MakeIterator();
1205       histList = new TList;
1206       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1207          // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
1208          hist = (TH3F*)(objarray->At(i));
1209          histList->Add(hist);
1210       }
1211       ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
1212       delete histList;
1213       delete objListIterator;
1214    }
1215
1216    if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
1217    // merge fArrayQDZ
1218    for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1219       objListIterator = arrayQDZList->MakeIterator();
1220       histList = new TList;
1221       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1222          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1223          if (!objarray) continue;
1224          hist = (TH3F*)(objarray->At(i));
1225          histList->Add(hist);
1226       }
1227       ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
1228       delete histList;
1229       delete objListIterator;
1230    }
1231
1232    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
1233    // merge fArrayQRMSY
1234    for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
1235       objListIterator = arrayQRMSYList->MakeIterator();
1236       histList = new TList;
1237       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1238          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1239          if (!objarray) continue;
1240          hist = (TH3F*)(objarray->At(i));
1241          histList->Add(hist);
1242       }
1243       ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
1244       delete histList;
1245       delete objListIterator;
1246    }   
1247
1248    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
1249    // merge fArrayQRMSZ
1250    for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
1251       objListIterator = arrayQRMSZList->MakeIterator();
1252       histList = new TList;
1253       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1254          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1255          if (!objarray) continue;
1256          hist = (TH3F*)(objarray->At(i));
1257          histList->Add(hist);
1258       }
1259       ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
1260       delete histList;
1261       delete objListIterator;
1262    } 
1263    
1264    
1265   
1266    
1267         
1268    
1269    if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
1270    // merge fResolY
1271    for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
1272       objListIterator = resolYList->MakeIterator();
1273       histList = new TList;
1274       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1275          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1276          if (!objarray) continue;
1277          hist = (TH3F*)(objarray->At(i));
1278          histList->Add(hist);
1279       }
1280       ((TH3F*)(fResolY->At(i)))->Merge(histList);
1281       delete histList;
1282       delete objListIterator;
1283    }
1284    
1285    // merge fResolZ
1286    for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1287       objListIterator = resolZList->MakeIterator();
1288       histList = new TList;
1289       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1290          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1291          if (!objarray) continue;
1292          hist = (TH3F*)(objarray->At(i));
1293          histList->Add(hist);
1294       }
1295       ((TH3F*)(fResolZ->At(i)))->Merge(histList);
1296       delete histList;
1297       delete objListIterator;
1298    }
1299
1300    // merge fRMSY
1301    for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
1302       objListIterator = rMSYList->MakeIterator();
1303       histList = new TList;
1304       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1305          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1306          if (!objarray) continue;
1307          hist = (TH3F*)(objarray->At(i));
1308          histList->Add(hist);
1309       }
1310       ((TH3F*)(fRMSY->At(i)))->Merge(histList);
1311       delete histList;
1312       delete objListIterator;
1313    }
1314          
1315    // merge fRMSZ
1316    for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
1317       objListIterator = rMSZList->MakeIterator();
1318       histList = new TList;
1319       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
1320          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
1321          if (!objarray) continue;
1322          hist = (TH3F*)(objarray->At(i));
1323          histList->Add(hist);
1324       }
1325       ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
1326       delete histList;
1327       delete objListIterator;
1328    }
1329    
1330    delete deltaYList;
1331    delete deltaZList;
1332    delete arrayAmpRowList;
1333    delete arrayAmpList;
1334    delete arrayQDYList;
1335    delete arrayQDZList;
1336    delete arrayQRMSYList;
1337    delete arrayQRMSZList;
1338    delete resolYList;
1339    delete resolZList;
1340    delete rMSYList;
1341    delete rMSZList;
1342    delete listIterator;
1343    
1344    if (GetDebugLevel() > 0) cout << "merging done!" << endl;
1345    
1346    return 1;
1347 }
1348
1349
1350
1351
1352 void AliTPCcalibTracks::MakeHistos(){
1353   //
1354   ////make THnSparse
1355   //
1356   //THnSparse  *fHisDeltaY;    // THnSparse - delta Y 
1357   //THnSparse  *fHisDeltaZ;    // THnSparse - delta Z 
1358   //THnSparse  *fHisRMSY;      // THnSparse - rms Y 
1359   //THnSparse  *fHisRMSZ;      // THnSparse - rms Z 
1360   //THnSparse  *fHisQmax;      // THnSparse - qmax 
1361   //THnSparse  *fHisQtot;      // THnSparse - qtot 
1362   // cluster  performance bins
1363   // 0 - variable of interest
1364   // 1 - pad type   - 0- short 1-medium 2-long pads
1365   // 2 - drift length - drift length -0-1
1366   // 3 - Qmax         - Qmax  - 2- 400
1367   // 4 - cog          - COG position - 0-1
1368   // 5 - tan(phi)     - local y angle
1369   // 6 - tan(theta)   - local z angle
1370   // 7 - sector       - sector number
1371   Double_t xminTrack[8], xmaxTrack[8];
1372   Int_t binsTrack[8];
1373   TString axisName[8];
1374   
1375   //
1376   binsTrack[0]=200;
1377   axisName[0]  ="var";
1378
1379   binsTrack[1] =3;
1380   xminTrack[1] =0; xmaxTrack[1]=3;
1381   axisName[1]  ="pad type";
1382   //
1383   binsTrack[2] =20;
1384   xminTrack[2] =-250; xmaxTrack[2]=250;
1385   axisName[2]  ="z";
1386   //
1387   binsTrack[3] =10;
1388   xminTrack[3] =1; xmaxTrack[3]=400;
1389   axisName[3]  ="Qmax";
1390   //
1391   binsTrack[4] =20;
1392   xminTrack[4] =0; xmaxTrack[4]=1;
1393   axisName[4]  ="cog";
1394   //
1395   binsTrack[5] =15;
1396   xminTrack[5] =-1.5; xmaxTrack[5]=1.5;
1397   axisName[5]  ="tan(angle)";
1398   //
1399   //
1400   xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1401   fHisDeltaY=new THnSparseF("#Delta_{y} (cm)","#Delta_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1402   xminTrack[0] =-1.5; xmaxTrack[0]=1.5;
1403   fHisDeltaZ=new THnSparseF("#Delta_{z} (cm)","#Delta_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1404   xminTrack[0] =0.; xmaxTrack[0]=0.5;
1405   fHisRMSY=new THnSparseF("#RMS_{y} (cm)","#RMS_{y} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1406   xminTrack[0] =0.; xmaxTrack[0]=0.5;
1407   fHisRMSZ=new THnSparseF("#RMS_{z} (cm)","#RMS_{z} (cm)", 6, binsTrack,xminTrack, xmaxTrack);
1408   xminTrack[0] =0.; xmaxTrack[0]=100;
1409   fHisQmax=new THnSparseF("Qmax (ADC)","Qmax (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
1410
1411   xminTrack[0] =0.; xmaxTrack[0]=250;
1412   fHisQtot=new THnSparseF("Qtot (ADC)","Qtot (ADC)", 6, binsTrack,xminTrack, xmaxTrack);
1413
1414
1415   for (Int_t ivar=0;ivar<6;ivar++){
1416     fHisDeltaY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1417     fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1418     fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1419     fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1420     fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1421     fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1422     fHisDeltaY->GetAxis(ivar)->SetTitle(axisName[ivar].Data());
1423     fHisDeltaZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1424     fHisRMSY->GetAxis(ivar)->SetName(axisName[ivar].Data());
1425     fHisRMSZ->GetAxis(ivar)->SetName(axisName[ivar].Data());
1426     fHisQmax->GetAxis(ivar)->SetName(axisName[ivar].Data());
1427     fHisQtot->GetAxis(ivar)->SetName(axisName[ivar].Data());
1428   }
1429
1430
1431   BinLogX(fHisDeltaY,3);
1432   BinLogX(fHisDeltaZ,3);
1433   BinLogX(fHisRMSY,3);
1434   BinLogX(fHisRMSZ,3);
1435   BinLogX(fHisQmax,3);
1436   BinLogX(fHisQtot,3);
1437
1438 }  
1439
1440 void    AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){
1441   //
1442   // Add histograms
1443   //
1444   if (calib->fHisDeltaY) fHisDeltaY->Add(calib->fHisDeltaY);
1445   if (calib->fHisDeltaZ) fHisDeltaZ->Add(calib->fHisDeltaZ);
1446   if (calib->fHisRMSY)   fHisRMSY->Add(calib->fHisRMSY);
1447   if (calib->fHisRMSZ)   fHisRMSZ->Add(calib->fHisRMSZ);
1448 }
1449
1450
1451
1452 void AliTPCcalibTracks::MakeSummaryTree(THnSparse *hisInput, TTreeSRedirector *pcstream, Int_t ptype){
1453   //
1454   // Dump summary info
1455   //
1456   //  0.OBJ: TAxis     var
1457   //  1.OBJ: TAxis     pad type
1458   //  2.OBJ: TAxis     z
1459   //  3.OBJ: TAxis     Qmax
1460   //  4.OBJ: TAxis     cog
1461   //  5.OBJ: TAxis     tan(angle)
1462   //
1463   if (ptype>3) return;
1464   Int_t idim[6]={0,1,2,3,4,5};
1465   TString hname[4]={"dy","dz","rmsy","rmsz"};
1466   //
1467   Int_t nbins5=hisInput->GetAxis(5)->GetNbins();
1468   Int_t first5=hisInput->GetAxis(5)->GetFirst();
1469   Int_t last5 =hisInput->GetAxis(5)->GetLast();
1470   //
1471   for (Int_t ibin5=first5-1; ibin5<=last5; ibin5+=1){   // axis 5 - local angle
1472     Bool_t bin5All=kTRUE;
1473     if (ibin5>=first5){
1474       hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5-1,first5),TMath::Min(ibin5+1,nbins5));
1475       if (ibin5==first5) hisInput->GetAxis(5)->SetRange(TMath::Max(ibin5,first5),TMath::Min(ibin5,nbins5));
1476       bin5All=kFALSE;
1477     }
1478     Double_t      x5= hisInput->GetAxis(5)->GetBinCenter(ibin5);
1479     THnSparse * his5= hisInput->Projection(5,idim);         //projected histogram according selection 5    
1480     Int_t nbins4=his5->GetAxis(4)->GetNbins();
1481     Int_t first4=his5->GetAxis(4)->GetFirst();
1482     Int_t last4 =his5->GetAxis(4)->GetLast();
1483     
1484     for (Int_t ibin4=first4-1; ibin4<=last4; ibin4+=1){   // axis 4 - cog
1485       Bool_t bin4All=kTRUE;
1486       if (ibin4>=first4){
1487         his5->GetAxis(4)->SetRange(TMath::Max(ibin4+1,first4),TMath::Min(ibin4-1,nbins4));
1488         if (ibin4==first4||ibin4==last4) his5->GetAxis(4)->SetRange(TMath::Max(ibin4,first4),TMath::Min(ibin4,nbins4));
1489         bin4All=kFALSE;
1490       }
1491       Double_t      x4= his5->GetAxis(4)->GetBinCenter(ibin4);
1492       THnSparse * his4= his5->Projection(4,idim);         //projected histogram according selection 4
1493       //
1494       Int_t nbins3=his4->GetAxis(3)->GetNbins();
1495       Int_t first3=his4->GetAxis(3)->GetFirst();
1496       Int_t last3 =his4->GetAxis(3)->GetLast();
1497       //
1498       for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){   // axis 3 - Qmax
1499         Bool_t bin3All=kTRUE;
1500         if (ibin3>=first3){
1501           his4->GetAxis(3)->SetRange(TMath::Max(ibin3,first3),TMath::Min(ibin3,nbins3));
1502           bin3All=kFALSE;
1503         }
1504         Double_t      x3= his4->GetAxis(3)->GetBinCenter(ibin3);
1505         THnSparse * his3= his4->Projection(3,idim);         //projected histogram according selection 3
1506         //
1507         Int_t nbins2    = his3->GetAxis(2)->GetNbins();
1508         Int_t first2    = his3->GetAxis(2)->GetFirst();
1509         Int_t last2     = his3->GetAxis(2)->GetLast();
1510         //
1511         for (Int_t ibin2=first2-1; ibin2<=last2; ibin2+=1){   // axis 2 - z     
1512           Bool_t bin2All=kTRUE;
1513           Double_t      x2= his3->GetAxis(2)->GetBinCenter(ibin2);
1514           if (ibin2>=first2){
1515             his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,first2),TMath::Min(ibin2+1,nbins2));
1516             if (ibin2==first2||ibin2==last2||TMath::Abs(x2)<20) his3->GetAxis(2)->SetRange(TMath::Max(ibin2,first2),TMath::Min(ibin2,nbins2));
1517             bin2All=kFALSE;
1518           }
1519           THnSparse * his2= his3->Projection(2,idim);         //projected histogram according selection 2
1520           //
1521           Int_t nbins1     = his2->GetAxis(1)->GetNbins();
1522           Int_t first1     = his2->GetAxis(1)->GetFirst();
1523           Int_t last1      = his2->GetAxis(1)->GetLast();
1524           for (Int_t ibin1=first1-1; ibin1<=last1; ibin1++){   //axis 1 - pad type 
1525             Bool_t bin1All=kTRUE;
1526             if (ibin1>=first1){
1527               his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));           
1528               bin1All=kFALSE;
1529             }
1530             Double_t       x1= TMath::Nint(his2->GetAxis(1)->GetBinCenter(ibin1)-0.5);
1531             TH1 * hisDelta = his2->Projection(0);
1532             Double_t entries = hisDelta->GetEntries();
1533             Double_t mean=0, rms=0;
1534             if (entries>10){
1535               mean    = hisDelta->GetMean();
1536               rms = hisDelta->GetRMS();   
1537               hisDelta->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
1538               mean    = hisDelta->GetMean();
1539               rms = hisDelta->GetRMS();
1540             }
1541             //
1542             (*pcstream)<<hname[ptype].Data()<<
1543               // flag - intgrated values for given bin
1544               "angleA="<<bin5All<<
1545               "cogA="<<bin4All<<
1546               "qmaxA="<<bin3All<<
1547               "zA="<<bin2All<<
1548               "ipadA="<<bin1All<<
1549               // center of bin value
1550               "angle="<<x5<<
1551               "cog="<<x4<<
1552               "qmax="<<x3<<
1553               "z="<<x2<<
1554               "ipad="<<x1<<
1555               // mean values
1556               //
1557               "entries="<<entries<<
1558               "mean="<<mean<<
1559               "rms="<<rms<<
1560               "ptype="<<ptype<<   //
1561               "\n";
1562             delete hisDelta;
1563             printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n",x5,x4,x3,x2,x1, entries,mean);          
1564           }//loop z
1565           delete his2;
1566         } // loop Q max
1567         delete his3;
1568       } // loop COG
1569       delete his4;
1570     }//loop local angle
1571     delete his5;
1572   }
1573 }