Adding the global event information in central place
[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 How to retrive it from file (created using calibration task):
40
41 gSystem->Load("libANALYSIS");
42 gSystem->Load("libTPCcalib");
43 TFile fcalib("CalibObjects.root");
44 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
45 AliTPCcalibTracks * calibTracks = ( AliTPCcalibTracks *)array->FindObject("calibTracks");
46
47
48 //USAGE of debug stream example
49  gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
50   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
51   AliXRDPROOFtoolkit tool;
52   TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
53   chainres->Lookup();
54 */
55
56
57                                                                //
58 ///////////////////////////////////////////////////////////////////////////////
59
60 //
61 // ROOT includes 
62 //
63 #include <iostream>
64 #include <fstream>
65 using namespace std;
66
67 #include <TROOT.h>
68 #include <TChain.h>
69 #include <TFile.h>
70 #include <TH3F.h>
71 #include <TProfile.h>
72
73 //
74 //#include <TPDGCode.h>
75 #include <TStyle.h>
76 #include "TLinearFitter.h"
77 //#include "TMatrixD.h"
78 #include "TTreeStream.h"
79 #include "TF1.h"
80 #include <TCanvas.h>
81 #include <TGraph2DErrors.h>
82 #include "TPostScript.h"
83 #include "TCint.h"
84
85 #include <TH2D.h>
86 #include <TF2.h>
87 #include <TSystem.h>
88 #include <TCollection.h>
89 #include <iostream>
90 #include <TLinearFitter.h>
91 #include <TString.h>
92
93 //
94 // AliROOT includes 
95 //
96 #include "AliMagF.h"
97 #include "AliTracker.h"
98 #include "AliESD.h"
99 #include "AliESDtrack.h"
100 #include "AliESDfriend.h"
101 #include "AliESDfriendTrack.h" 
102 #include "AliTPCseed.h"
103 #include "AliTPCclusterMI.h"
104 #include "AliTPCROC.h"
105
106 #include "AliTPCParamSR.h"
107 #include "AliTrackPointArray.h"
108 #include "AliTPCcalibTracks.h"
109 #include "AliTPCClusterParam.h"
110 #include "AliTPCcalibTracksCuts.h"
111 #include "AliTPCCalPadRegion.h"
112 #include "AliTPCCalPad.h"
113 #include "AliTPCCalROC.h"
114 #include "TText.h"
115 #include "TPaveText.h"
116 #include "TSystem.h"
117 #include "TStatToolkit.h"
118 #include "TCut.h"
119
120
121
122 ClassImp(AliTPCcalibTracks)
123
124
125 AliTPCcalibTracks::AliTPCcalibTracks():
126   AliTPCcalibBase(),
127   fClusterParam(0),
128   fROC(0),
129   fArrayAmpRow(0),
130   fArrayAmp(0), 
131   fArrayQDY(0), 
132   fArrayQDZ(0), 
133   fArrayQRMSY(0),
134   fArrayQRMSZ(0),
135   fArrayChargeVsDriftlength(0),
136   fcalPadRegionChargeVsDriftlength(0),
137   fDeltaY(0),
138   fDeltaZ(0),
139   fResolY(0),
140   fResolZ(0),
141   fRMSY(0),
142   fRMSZ(0),
143   fCuts(0),
144   fHclus(0),
145   fRejectedTracksHisto(0),
146   fHclusterPerPadrow(0),
147   fHclusterPerPadrowRaw(0),
148   fClusterCutHisto(0),
149   fCalPadClusterPerPad(0),
150   fCalPadClusterPerPadRaw(0)
151
152    // 
153    // AliTPCcalibTracks default constructor
154    //    
155   SetDebugLevel(1);
156   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl;  
157 }   
158
159
160
161 AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
162   AliTPCcalibBase(calibTracks),
163   fClusterParam(0),
164   fROC(0),
165   fArrayAmpRow(0),
166   fArrayAmp(0), 
167   fArrayQDY(0), 
168   fArrayQDZ(0), 
169   fArrayQRMSY(0),
170   fArrayQRMSZ(0),
171   fArrayChargeVsDriftlength(0),
172   fcalPadRegionChargeVsDriftlength(0),
173   fDeltaY(0),
174   fDeltaZ(0),
175   fResolY(0),
176   fResolZ(0),
177   fRMSY(0),
178   fRMSZ(0),
179   fCuts(0),
180   fHclus(0),
181   fRejectedTracksHisto(0),
182   fHclusterPerPadrow(0),
183   fHclusterPerPadrowRaw(0),
184   fClusterCutHisto(0),
185   fCalPadClusterPerPad(0),
186   fCalPadClusterPerPadRaw(0)
187 {
188    // 
189    // AliTPCcalibTracks copy constructor
190    // 
191   if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl;
192    
193    Bool_t dirStatus = TH1::AddDirectoryStatus();
194    TH1::AddDirectory(kFALSE);
195    
196    Int_t length = -1;
197    // backward compatibility: if the data member doesn't yet exist, it will not be merged
198    (calibTracks.fArrayAmpRow) ? length = calibTracks.fArrayAmpRow->GetEntriesFast() : length = -1;
199    fArrayAmpRow = new TObjArray(length);
200    fArrayAmp = new TObjArray(length);
201    for (Int_t i = 0; i < length; i++) {
202       fArrayAmpRow->AddAt( (TProfile*)calibTracks.fArrayAmpRow->At(i)->Clone(), i);
203       fArrayAmp->AddAt( ((TProfile*)calibTracks.fArrayAmp->At(i)->Clone()), i);
204    }
205    
206    (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
207    fArrayQDY= new TObjArray(length);
208    fArrayQDZ= new TObjArray(length);
209    fArrayQRMSY= new TObjArray(length);
210    fArrayQRMSZ= new TObjArray(length);
211    for (Int_t i = 0; i < length; i++) {
212       fArrayQDY->AddAt( ((TH1F*)calibTracks.fArrayQDY->At(i)->Clone()), i);
213       fArrayQDZ->AddAt( ((TH1F*)calibTracks.fArrayQDZ->At(i)->Clone()), i);
214       fArrayQRMSY->AddAt( ((TH1F*)calibTracks.fArrayQRMSY->At(i)->Clone()), i);
215       fArrayQRMSZ->AddAt( ((TH1F*)calibTracks.fArrayQRMSZ->At(i)->Clone()), i);
216    }
217    
218    (calibTracks.fResolY) ? length = calibTracks.fResolY->GetEntriesFast() : length = -1;
219    fResolY = new TObjArray(length);
220    fResolZ = new TObjArray(length);
221    fRMSY = new TObjArray(length);
222    fRMSZ = new TObjArray(length);
223    for (Int_t i = 0; i < length; i++) {
224       fResolY->AddAt( ((TH1F*)calibTracks.fResolY->At(i)->Clone()), i);
225       fResolZ->AddAt( ((TH1F*)calibTracks.fResolZ->At(i)->Clone()), i);
226       fRMSY->AddAt( ((TH1F*)calibTracks.fRMSY->At(i)->Clone()), i);
227       fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
228    } 
229    
230    (calibTracks.fArrayChargeVsDriftlength) ? length = calibTracks.fArrayChargeVsDriftlength->GetEntriesFast() : length = -1;
231    (calibTracks.fArrayChargeVsDriftlength) ? fArrayChargeVsDriftlength = new TObjArray(length) : fArrayChargeVsDriftlength = 0;
232    for (Int_t i = 0; i < length; i++) {
233       fArrayChargeVsDriftlength->AddAt( ((TProfile*)calibTracks.fArrayChargeVsDriftlength->At(i)->Clone()), i);
234    }
235    
236    fDeltaY =  (TH1F*)calibTracks.fDeltaY->Clone();
237    fDeltaZ =  (TH1F*)calibTracks.fDeltaZ->Clone();
238    fHclus = (TH1I*)calibTracks.fHclus->Clone();
239    fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
240    fRejectedTracksHisto    = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
241    fHclusterPerPadrow      = (TH1I*)calibTracks.fHclusterPerPadrow->Clone();
242    fHclusterPerPadrowRaw   = (TH1I*)calibTracks.fHclusterPerPadrowRaw->Clone();
243    fcalPadRegionChargeVsDriftlength = (AliTPCCalPadRegion*)calibTracks.fcalPadRegionChargeVsDriftlength->Clone();
244    fCalPadClusterPerPad    = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
245    fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
246
247    fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(), 
248       calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise());
249    SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle());
250    TH1::AddDirectory(dirStatus); // set status back to original status
251 //    cout << "+++++ end of copy constructor +++++" << endl;   // TO BE REMOVED
252 }
253
254
255 AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibTracks){
256   //
257   // assgnment operator
258   //
259   if (this != &calibTracks) {
260     new (this) AliTPCcalibTracks(calibTracks);
261   }
262   return *this;
263
264 }
265
266
267 AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam,  AliTPCcalibTracksCuts* cuts, Int_t logLevel) : 
268   AliTPCcalibBase(),
269   fClusterParam(0),
270   fROC(0),
271   fArrayAmpRow(0),
272   fArrayAmp(0), 
273   fArrayQDY(0), 
274   fArrayQDZ(0), 
275   fArrayQRMSY(0),
276   fArrayQRMSZ(0),
277   fArrayChargeVsDriftlength(0),
278   fcalPadRegionChargeVsDriftlength(0),
279   fDeltaY(0),
280   fDeltaZ(0),
281   fResolY(0),
282   fResolZ(0),
283   fRMSY(0),
284   fRMSZ(0),
285   fCuts(0),
286   fHclus(0),
287   fRejectedTracksHisto(0),
288   fHclusterPerPadrow(0),
289   fHclusterPerPadrowRaw(0),
290   fClusterCutHisto(0),
291   fCalPadClusterPerPad(0),
292   fCalPadClusterPerPadRaw(0)
293  {
294    // 
295    // AliTPCcalibTracks constructor
296    // specify 'name' and 'title' of your object
297    // specify 'clusterParam', (needed for TPC cluster error and shape parameterization)
298    // In the parameter 'cuts' the cuts are specified, that decide           
299    // weather a track will be accepted for calibration or not.              
300    //
301    // fDebugLevel - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStreamer, 6: waste your screen
302    // 
303    // All histograms are instatiated in this constructor.
304    // 
305    this->SetName(name);
306    this->SetTitle(title);
307
308    if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl;
309    G__SetCatchException(0);     
310    
311    fClusterParam = clusterParam;
312    if (fClusterParam){
313      fClusterParam->SetInstance(fClusterParam);
314    }
315    else {
316      Error("AliTPCcalibTracks","No cluster parametrization found! A valid clusterParam object is needed in the constructor. (To be found in 'TPCClusterParam.root'.)");
317    } 
318    fCuts = cuts;
319    SetDebugLevel(logLevel);
320    
321    TH1::AddDirectory(kFALSE);
322    
323    char chname[1000];
324    TProfile * prof1=0;
325    TH1F     * his1 =0;
326    fHclus = new TH1I("hclus","Number of clusters per track",160, 0, 160);     // valgrind 3
327    fRejectedTracksHisto    = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
328    fHclusterPerPadrow      = new TH1I("fHclusterPerPadrow", " clusters per padRow, used for the resolution tree", 160, 0, 160);
329    fHclusterPerPadrowRaw   = new TH1I("fHclusterPerPadrowRaw", " clusters per padRow, before cutting clusters", 160, 0, 160);
330    fCalPadClusterPerPad    = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
331    fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
332    fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
333    
334    // Amplitude  - sector - row histograms 
335    fArrayAmpRow = new TObjArray(72);
336    fArrayAmp    = new TObjArray(72);
337    fArrayChargeVsDriftlength = new TObjArray(72);
338    
339    for (Int_t i = 0; i < 36; i++){   
340       sprintf(chname,"Amp_row_Sector%d",i);
341       prof1 = new TProfile(chname,chname,63,0,64);          // valgrind 3   193,536 bytes in 354 blocks are still reachable 
342       prof1->SetXTitle("Pad row");
343       prof1->SetYTitle("Mean Max amplitude");
344       fArrayAmpRow->AddAt(prof1,i);
345       sprintf(chname,"Amp_row_Sector%d",i+36);
346       prof1 = new TProfile(chname,chname,96,0,97);       // valgrind 3   3,912 bytes in 6 blocks are possibly lost
347       prof1->SetXTitle("Pad row");
348       prof1->SetYTitle("Mean Max  amplitude");
349       fArrayAmpRow->AddAt(prof1,i+36);
350       
351       // amplitude
352       sprintf(chname,"Amp_Sector%d",i);
353       his1 = new TH1F(chname,chname,250,0,500);         // valgrind 
354       his1->SetXTitle("Max Amplitude (ADC)");
355       fArrayAmp->AddAt(his1,i);
356       sprintf(chname,"Amp_Sector%d",i+36);
357       his1 = new TH1F(chname,chname,200,0,600);         // valgrind 3   13,408,208 bytes in 229 blocks are still reachable
358       his1->SetXTitle("Max Amplitude (ADC)");
359       fArrayAmp->AddAt(his1,i+36);
360       
361       // driftlength
362       sprintf(chname, "driftlengt vs. charge, ROC %i", i);
363       prof1 = new TProfile(chname, chname, 500, 0, 250);
364       prof1->SetYTitle("Charge");
365       prof1->SetXTitle("Driftlength");
366       fArrayChargeVsDriftlength->AddAt(prof1,i);
367       sprintf(chname, "driftlengt vs. charge, ROC %i", i+36);
368       prof1 = new TProfile(chname, chname, 500, 0, 250);
369       prof1->SetYTitle("Charge");
370       prof1->SetXTitle("Driftlength");
371       fArrayChargeVsDriftlength->AddAt(prof1,i+36);
372    }
373    
374    TH1::AddDirectory(kFALSE);
375    
376    fDeltaY = new TH1F("DeltaY","DeltaY",100,-1,1);
377    fDeltaZ = new TH1F("DeltaZ","DeltaZ",100,-1,1);
378    
379    fResolY = new TObjArray(3);
380    fResolZ = new TObjArray(3);
381    fRMSY   = new TObjArray(3);
382    fRMSZ   = new TObjArray(3);
383    TH3F * his3D;
384    //
385    his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
386    fResolY->AddAt(his3D,0);     
387    his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
388    fResolY->AddAt(his3D,1);
389    his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
390    fResolY->AddAt(his3D,2);
391    //
392    his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
393    fResolZ->AddAt(his3D,0);
394    his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
395    fResolZ->AddAt(his3D,1);
396    his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
397    fResolZ->AddAt(his3D,2);
398    //
399    his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
400    fRMSY->AddAt(his3D,0);
401    his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
402    fRMSY->AddAt(his3D,1);
403    his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
404    fRMSY->AddAt(his3D,2);
405    //
406    his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
407    fRMSZ->AddAt(his3D,0);
408    his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
409    fRMSZ->AddAt(his3D,1);
410    his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
411    fRMSZ->AddAt(his3D,2);
412    //
413       
414    TH1::AddDirectory(kFALSE);
415    
416    fArrayQDY = new TObjArray(300);
417    fArrayQDZ = new TObjArray(300);
418    fArrayQRMSY = new TObjArray(300);
419    fArrayQRMSZ = new TObjArray(300);
420    for (Int_t iq = 0; iq <= 10; iq++){
421       for (Int_t ipad = 0; ipad < 3; ipad++){
422          Int_t   bin   = GetBin(iq, ipad);
423          Float_t qmean = GetQ(bin);
424          char name[200];
425          sprintf(name,"ResolY Pad%d Qmiddle%f",ipad, qmean);
426          his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1);
427          fArrayQDY->AddAt(his3D, bin);
428          sprintf(name,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
429          his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1);
430          fArrayQDZ->AddAt(his3D, bin);
431          sprintf(name,"RMSY Pad%d Qmiddle%f",ipad, qmean);
432          his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1);
433          fArrayQRMSY->AddAt(his3D, bin);
434          sprintf(name,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
435          his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1);
436          fArrayQRMSZ->AddAt(his3D, bin);
437       }
438    }
439    
440    fcalPadRegionChargeVsDriftlength = new AliTPCCalPadRegion("fcalPadRegionChargeVsDriftlength", "TProfiles with charge vs driftlength for each pad region");
441    TProfile *tempProf;
442    for (UInt_t padSize = 0; padSize < 3; padSize++) {
443       for (UInt_t isector = 0; isector < 36; isector++) {
444          if (padSize == 0) sprintf(chname, "driftlengt vs. charge, sector %i, short pads", isector);
445          if (padSize == 1) sprintf(chname, "driftlengt vs. charge, sector %i, medium  pads", isector);
446          if (padSize == 2) sprintf(chname, "driftlengt vs. charge, sector %i, long  pads", isector);
447          tempProf = new TProfile(chname, chname, 500, 0, 250);
448          tempProf->SetYTitle("Charge");
449          tempProf->SetXTitle("Driftlength");
450          fcalPadRegionChargeVsDriftlength->SetObject(tempProf, isector, padSize);
451       }
452    }
453    
454
455    if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl; 
456    cout << "end of main constructor" << endl; // TO BE REMOVED
457 }    
458
459
460 AliTPCcalibTracks::~AliTPCcalibTracks() {
461    // 
462    // AliTPCcalibTracks destructor
463    // 
464    
465   if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
466    Int_t length = 0;
467    if (fArrayAmpRow) length = fArrayAmpRow->GetEntriesFast();
468    for (Int_t i = 0; i < length; i++){
469       delete fArrayAmpRow->At(i);  
470       delete fArrayAmp->At(i);  
471    }
472    delete fArrayAmpRow;
473    delete fArrayAmp;
474    
475    delete fDeltaY;
476    delete fDeltaZ;
477    
478    if (fResolY) length = fResolY->GetEntriesFast();
479    for (Int_t i = 0; i < length; i++){
480       delete fResolY->At(i);
481       delete fResolZ->At(i);
482       delete fRMSY->At(i);
483       delete fRMSZ->At(i);
484    }
485    delete fResolY;
486    delete fResolZ;
487    delete fRMSY;
488    delete fRMSZ;
489    
490    if (fArrayQDY) length = fArrayQDY->GetEntriesFast();
491    for (Int_t i = 0; i < length; i++){
492       delete fArrayQDY->At(i);
493       delete fArrayQDZ->At(i);
494       delete fArrayQRMSY->At(i);
495       delete fArrayQRMSZ->At(i);
496    }
497    
498    if (fArrayChargeVsDriftlength) length = fArrayChargeVsDriftlength->GetEntriesFast();
499    for (Int_t i = 0; i < length; i++){
500       delete fArrayChargeVsDriftlength->At(i);
501    }
502    
503     
504    delete fArrayQDY;
505    delete fArrayQDZ;
506    delete fArrayQRMSY;
507    delete fArrayQRMSZ;
508    delete fArrayChargeVsDriftlength;
509    
510   delete fHclus;
511   delete fRejectedTracksHisto;
512   delete fClusterCutHisto;
513   delete fHclusterPerPadrow;
514   delete fHclusterPerPadrowRaw;
515   if (fCalPadClusterPerPad)    delete fCalPadClusterPerPad;
516   if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
517   if(fcalPadRegionChargeVsDriftlength) {
518      fcalPadRegionChargeVsDriftlength->Delete();
519      delete fcalPadRegionChargeVsDriftlength;
520   }
521 }
522    
523   
524
525 void AliTPCcalibTracks::Process(AliTPCseed *track){
526    // 
527    // To be called in the selector
528    // first AcceptTrack is evaluated, then calls all the following analyse functions: 
529    // FillResolutionHistoLocal(track)
530    // AlignUpDown(track, esd)
531    // 
532   if (GetDebugLevel() > 5) Info("Process","Starting to process the track...");
533    Int_t accpetStatus = AcceptTrack(track);
534    if (accpetStatus == 0) {
535       FillResolutionHistoLocal(track);
536       // AlignUpDown(track, esd);
537    }
538    else fRejectedTracksHisto->Fill(accpetStatus);
539 }
540
541
542
543 Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){
544   //
545   // calculate bins for given q and pad type 
546   // used in TObjArray
547   //
548   Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 );  
549   res *= 3;
550   res += pad;
551   return res;
552 }
553
554
555 Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){
556   //
557   // calculate bins for given iq and pad type 
558   // used in TObjArray
559   //
560   return iq * 3 + pad;;
561 }
562
563
564 Float_t AliTPCcalibTracks::GetQ(Int_t bin){
565    // 
566    // returns to bin belonging charge
567    // (bin / 3 + 3)^2
568    // 
569    Int_t bin0 = bin / 3;
570    bin0 += 3;
571    return bin0 * bin0;
572 }
573
574
575 Float_t AliTPCcalibTracks::GetPad(Int_t bin){
576    // 
577    // returns to bin belonging pad
578    // bin % 3
579    // 
580    return bin % 3; 
581 }
582
583
584
585 Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
586   //
587   // Function, that decides wheather a given track is accepted for 
588   // the analysis or not. 
589   // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts'
590   // Returns 0 if a track is accepted or an integer different from 0 
591   // to indicate the failed cut
592   //
593   const Int_t   kMinClusters  = fCuts->GetMinClusters();
594   const Float_t kMinRatio     = fCuts->GetMinRatio();
595   const Float_t kMax1pt       = fCuts->GetMax1pt();
596   const Float_t kEdgeYXCutNoise    = fCuts->GetEdgeYXCutNoise();
597   const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise();
598   
599   //
600   // edge induced noise tracks - NEXT RELEASE will be removed during tracking
601   if ( TMath::Abs(track->GetY() / track->GetX()) > kEdgeYXCutNoise )
602     if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return 1;
603   if (track->GetNumberOfClusters() < kMinClusters) return 2;
604   Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.);
605   if (ratio < kMinRatio) return 3;
606 //   Float_t mpt = track->Get1Pt();       // Get1Pt() doesn't exist any more
607   Float_t mpt = track->GetSigned1Pt();
608   if (TMath::Abs(mpt) > kMax1pt) return 4;
609   //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
610   //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
611   //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
612   
613   if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");  
614   return 0;
615 }
616
617
618 void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
619    //
620    // fill resolution histograms - localy - tracklet in the neighborhood
621    // write debug information to 'TPCSelectorDebug.root'
622    // 
623    // _ the main function, called during track analysis _
624    // 
625    // loop over all padrows along the track
626    // fit tracklets (length: 13 clusters) calculate mean chi^2 for this track-fit in Y and Z direction
627    // 
628    // loop again over all padrows along the track
629    // fit tracklet (clusters in given padrow +- kDelta padrows) 
630    // with polynom of 2nd order and two polynoms of 1st order
631    // take both polynoms of 1st order, calculate difference of their parameters
632    // add covariance matrixes and calculate chi2 of this difference
633    // if this chi2 is bigger than a given threshold, assume that the current cluster is
634    // a kink an goto next padrow
635    // if not:
636    // fill fArrayAmpRow, array with amplitudes vs. row for given sector
637    // fill fArrayAmp, array with amplitude histograms for give sector
638    // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fDeltaY, fDeltaZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
639    // 
640    // write debug information to 'TPCSelectorDebug.root'
641    // only for every kDeltaWriteDebugStream'th padrow to reduce data volume 
642    // and to avoid redundant data
643    // 
644
645   static TLinearFitter fFitterLinY1(2,"pol1");   //
646   static TLinearFitter fFitterLinZ1(2,"pol1");   // 
647   static TLinearFitter fFitterLinY2(2,"pol1");   // 
648   static TLinearFitter fFitterLinZ2(2,"pol1");   //
649   static TLinearFitter fFitterParY(3,"pol2");    // 
650   static TLinearFitter fFitterParZ(3,"pol2");    //
651
652   fFitterLinY1.StoreData(kFALSE);
653   fFitterLinZ1.StoreData(kFALSE);
654   fFitterLinY2.StoreData(kFALSE);
655   fFitterLinZ2.StoreData(kFALSE);
656   fFitterParY.StoreData(kFALSE);
657   fFitterParZ.StoreData(kFALSE);
658
659
660   if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
661    const Int_t   kDelta    = 10;          // delta rows to fit
662    const Float_t kMinRatio = 0.75;        // minimal ratio
663    //   const Float_t kCutChi2  = 6.;          // cut chi2 - left right  - kink removal
664    const Float_t kErrorFraction = 0.5;    // use only clusters with small interpolation error - for error param
665    const Int_t   kFirstLargePad = 127;    // medium pads -> long pads
666    const Float_t kLargePadSize  = 1.5;    // factor between medium and long pads' area
667    const Int_t   kDeltaWriteDebugStream  = 5;  // only for every kDeltaWriteDebugStream'th padrow debug information is calulated and written to debugstream
668    TVectorD paramY0(2);
669    TVectorD paramZ0(2);
670    TVectorD paramY1(2);
671    TVectorD paramZ1(2);
672    TVectorD paramY2(3);
673    TVectorD paramZ2(3);
674    TMatrixD matrixY0(2,2);
675    TMatrixD matrixZ0(2,2);
676    TMatrixD matrixY1(2,2);
677    TMatrixD matrixZ1(2,2);
678    
679    // estimate mean error
680    Int_t nTrackletsAll = 0;       // number of tracklets for given track
681    Float_t csigmaY     = 0;       // mean sigma for tracklet refit in Y direction
682    Float_t csigmaZ     = 0;       // mean sigma for tracklet refit in Z direction
683    Int_t nClusters     = 0;       // working variable, number of clusters per tracklet
684    Int_t sectorG       = -1;      // working variable, sector of tracklet, has to stay constant for one tracklet
685    
686    fHclus->Fill(track->GetNumberOfClusters());      // for statistics overview
687    // ---------------------------------------------------------------------
688    for (Int_t irow = 0; irow < 159; irow++){
689       // loop over all rows along the track
690       // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
691       // calculate mean chi^2 for this track-fit in Y and Z direction
692       AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
693       if (!cluster0) continue;  // no cluster found
694       Int_t sector = cluster0->GetDetector();
695       fHclusterPerPadrowRaw->Fill(irow);
696       
697       Int_t ipad = TMath::Nint(cluster0->GetPad());
698       Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
699       fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
700       
701       if (sector != sectorG){
702          // track leaves sector before it crossed enough rows to fit / initialization
703          nClusters = 0;
704          fFitterParY.ClearPoints();
705          fFitterParZ.ClearPoints();
706          sectorG = sector;
707       }
708       else {
709          nClusters++;
710          Double_t x = cluster0->GetX();
711          fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
712          fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
713          //
714          if ( nClusters >= kDelta + 3 ){  
715          // if more than 13 (kDelta+3) clusters were added to the fitters
716          // fit the tracklet, increase trackletCounter
717          fFitterParY.Eval();
718          fFitterParZ.Eval();
719          nTrackletsAll++;
720          csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
721          csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
722          nClusters = -1;
723          fFitterParY.ClearPoints();
724          fFitterParZ.ClearPoints();
725          }
726       }
727    }      // for (Int_t irow = 0; irow < 159; irow++)
728    // mean chi^2 for all tracklet fits in Y and in Z direction: 
729    csigmaY = TMath::Sqrt(csigmaY / nTrackletsAll);
730    csigmaZ = TMath::Sqrt(csigmaZ / nTrackletsAll);
731    // ---------------------------------------------------------------------
732  
733    for (Int_t irow = 0; irow < 159; irow++){
734       // loop again over all rows along the track
735       // do analysis
736       // 
737       Int_t nclFound = 0;  // number of clusters in the neighborhood
738       Int_t ncl0 = 0;      // number of clusters in rows < rowOfCenterCluster
739       Int_t ncl1 = 0;      // number of clusters in rows > rowOfCenterCluster
740       AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
741       if (!cluster0) continue;
742       Int_t sector = cluster0->GetDetector();
743       Float_t xref = cluster0->GetX();
744          
745       // Make Fit
746       fFitterParY.ClearPoints();
747       fFitterParZ.ClearPoints();
748       fFitterLinY1.ClearPoints();
749       fFitterLinZ1.ClearPoints();
750       fFitterLinY2.ClearPoints();
751       fFitterLinZ2.ClearPoints();
752       
753       // fit tracklet (clusters in given padrow +- kDelta padrows) 
754       // with polynom of 2nd order and two polynoms of 1st order
755       // take both polynoms of 1st order, calculate difference of their parameters
756       // add covariance matrixes and calculate chi2 of this difference
757       // if this chi2 is bigger than a given threshold, assume that the current cluster is
758       // a kink an goto next padrow
759       
760       for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
761          // loop over irow +- kDelta rows (neighboured rows)
762          // 
763          // 
764          if (idelta == 0) continue;                                // don't use center cluster
765          if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
766          AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
767          if (!currentCluster) continue;
768          if (currentCluster->GetType() < 0) continue;
769          if (currentCluster->GetDetector() != sector) continue;
770          Double_t x = currentCluster->GetX() - xref;  // x = differece: current cluster - cluster @ irow
771          nclFound++;
772          if (idelta < 0){
773          ncl0++;
774          fFitterLinY1.AddPoint(&x, currentCluster->GetY(), csigmaY);
775          fFitterLinZ1.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
776          }
777          if (idelta > 0){
778          ncl1++;
779          fFitterLinY2.AddPoint(&x, currentCluster->GetY(), csigmaY);
780          fFitterLinZ2.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
781          }
782          fFitterParY.AddPoint(&x, currentCluster->GetY(), csigmaY);  
783          fFitterParZ.AddPoint(&x, currentCluster->GetZ(), csigmaZ);  
784       }  // loop over neighbourhood for fitter filling 
785
786
787       
788       if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
789       if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
790       if (nclFound < kDelta * kMinRatio) continue;    // if not enough clusters (7.5) found in neighbourhood goto next padrow
791       fFitterParY.Eval();
792       fFitterParZ.Eval();
793       Double_t chi2 = (fFitterParY.GetChisquare() + fFitterParZ.GetChisquare()) / (2. * nclFound - 6.);
794       //if (chi2 > kCutChi2) fRejectedTracksHisto->Fill(9);
795       //if (chi2 > kCutChi2) fClusterCutHisto->Fill(2, irow);
796       //if (chi2 > kCutChi2) continue;   // if chi^2 is too big goto next padrow
797       TTreeSRedirector *cstream = GetDebugStreamer();
798       if (cstream){
799         (*cstream)<<"Cut9"<<
800           "chi2="<<chi2<<
801           "\n";
802       }
803       // REMOVE KINK
804       // only when there are enough clusters (4) in each direction
805       if (ncl0 > 4){
806          fFitterLinY1.Eval();
807          fFitterLinZ1.Eval();
808       }
809       if (ncl1 > 4){
810          fFitterLinY2.Eval();
811          fFitterLinZ2.Eval();
812       }
813       
814       if (ncl0 > 4 && ncl1 > 4){
815          fFitterLinY1.GetCovarianceMatrix(matrixY0);
816          fFitterLinY2.GetCovarianceMatrix(matrixY1);
817          fFitterLinZ1.GetCovarianceMatrix(matrixZ0);
818          fFitterLinZ2.GetCovarianceMatrix(matrixZ1);
819          fFitterLinY2.GetParameters(paramY1);
820          fFitterLinZ2.GetParameters(paramZ1);
821          fFitterLinY1.GetParameters(paramY0);
822          fFitterLinZ1.GetParameters(paramZ0);
823          paramY0 -= paramY1;
824          paramZ0 -= paramZ1;
825          matrixY0 += matrixY1;
826          matrixZ0 += matrixZ1;
827          Double_t chi2 = 0;
828          
829          TMatrixD difY(2, 1, paramY0.GetMatrixArray());
830          TMatrixD difYT(1, 2, paramY0.GetMatrixArray());
831          matrixY0.Invert();
832          TMatrixD mulY(matrixY0, TMatrixD::kMult, difY);
833          TMatrixD chi2Y(difYT, TMatrixD::kMult, mulY);
834          chi2 += chi2Y(0, 0);
835          
836          TMatrixD difZ(2, 1, paramZ0.GetMatrixArray());
837          TMatrixD difZT(1, 2, paramZ0.GetMatrixArray());
838          matrixZ0.Invert();
839          TMatrixD mulZ(matrixZ0, TMatrixD::kMult, difZ);
840          TMatrixD chi2Z(difZT, TMatrixD::kMult, mulZ);
841          chi2 += chi2Z(0, 0);      
842          
843          // REMOVE KINK - TO be fixed - proper chi2 calculation for curved track to be implemented
844          //if (chi2 * 0.25 > kCutChi2) fRejectedTracksHisto->Fill(8);
845          //if (chi2 * 0.25 > kCutChi2) fClusterCutHisto->Fill(3, irow);
846          //if (chi2 * 0.25 > kCutChi2) continue;   // if chi2 is too big goto next padrow
847          // fit tracklet with polynom of 2nd order and two polynoms of 1st order
848          // take both polynoms of 1st order, calculate difference of their parameters
849          // add covariance matrixes and calculate chi2 of this difference
850          // if this chi2 is bigger than a given threshold, assume that the current cluster is
851          // a kink an goto next padrow
852
853          TTreeSRedirector *cstream = GetDebugStreamer();
854          if (cstream){
855            (*cstream)<<"Cut8"<<
856              "chi2="<<chi2<<
857              "\n";
858          }       
859       }
860       
861       // current padrow has no kink
862       
863       // get fit parameters from pol2 fit: 
864       Double_t paramY[4], paramZ[4];
865       paramY[0] = fFitterParY.GetParameter(0);
866       paramY[1] = fFitterParY.GetParameter(1);
867       paramY[2] = fFitterParY.GetParameter(2);
868       paramZ[0] = fFitterParZ.GetParameter(0);
869       paramZ[1] = fFitterParZ.GetParameter(1);
870       paramZ[2] = fFitterParZ.GetParameter(2);    
871       
872       Double_t tracky = paramY[0];
873       Double_t trackz = paramZ[0];
874       Float_t  deltay = tracky - cluster0->GetY();
875       Float_t  deltaz = trackz - cluster0->GetZ();
876       Float_t  angley = paramY[1] - paramY[0] / xref;
877       Float_t  anglez = paramZ[1];
878       
879       Float_t max = cluster0->GetMax();
880       UInt_t isegment = cluster0->GetDetector() % 36;
881       Int_t padSize = 0;                          // short pads
882       if (cluster0->GetDetector() >= 36) {
883          padSize = 1;                              // medium pads 
884          if (cluster0->GetRow() > 63) padSize = 2; // long pads
885       }
886
887       // =========================================
888       // wirte collected information to histograms
889       // =========================================
890       
891       TProfile *profAmpRow =  (TProfile*)fArrayAmpRow->At(sector);
892       if ( irow >= kFirstLargePad) max /= kLargePadSize;
893       profAmpRow->Fill( (Double_t)cluster0->GetRow(), max );
894       TH1F *hisAmp =  (TH1F*)fArrayAmp->At(sector);
895       hisAmp->Fill(max);
896       
897       // remove the following two lines one day:
898       TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector);
899       profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max );
900       
901       TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize));
902       profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max );
903       
904       fHclusterPerPadrow->Fill(irow);   // fill histogram showing clusters per padrow
905       Int_t ipad = TMath::Nint(cluster0->GetPad());
906       Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
907       fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
908    
909          
910       TH3F * his3 = 0;
911       his3 = (TH3F*)fRMSY->At(padSize);
912       if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
913       his3 = (TH3F*)fRMSZ->At(padSize);
914       if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
915       
916       his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
917       if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
918       his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
919       if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
920    
921       
922       // Fill resolution histograms
923       Bool_t useForResol = kTRUE;
924       if (fFitterParY.GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE;
925    
926       if (cstream){
927         Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
928         Float_t sy = cluster0->GetSigmaY2();
929         Float_t sz = cluster0->GetSigmaZ2();
930         (*cstream)<<"Resol0"<<
931           "run="<<fRun<<              //  run number
932           "event="<<fEvent<<          //  event number
933           "time="<<fTime<<            //  time stamp of event
934           "trigger="<<fTrigger<<      //  trigger
935           "mag="<<fMagF<<             //  magnetic field              
936           "padSize="<<padSize<<
937           "angley="<<angley<<
938           "anglez="<<anglez<<
939           "zdr="<<zdrift<<
940           "dy="<<deltay<<
941           "dz="<<deltaz<<
942           "sy="<<sy<<
943           "sz="<<sz<<
944           "\n";
945       }
946
947       if (useForResol){
948          fDeltaY->Fill(deltay);
949          fDeltaZ->Fill(deltaz);
950          his3 = (TH3F*)fResolY->At(padSize);
951          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
952          his3 = (TH3F*)fResolZ->At(padSize);
953          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
954          his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
955          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
956          his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
957          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
958       }
959       
960       //=============================================================================================
961       
962       if (useForResol && nclFound > 2 * kMinRatio * kDelta 
963           && irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){
964         if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
965          FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta);
966       }  // if (useForResol && nclFound > 2 * kMinRatio * kDelta)
967    
968    }    // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
969 }  // FillResolutionHistoLocal(...)
970
971
972
973 void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t  angley, Float_t  anglez, Int_t nclFound, Int_t kDelta) {
974    // 
975    //  - debug part of FillResolutionHistoLocal - 
976    // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data
977    // called only for GetStreamLevel() > 4
978    // fill resolution trees
979    //
980       
981    Int_t sector = cluster0->GetDetector();
982    Float_t xref = cluster0->GetX();
983    Int_t padSize = 0;                          // short pads
984    if (cluster0->GetDetector() >= 36) {
985       padSize = 1;                              // medium pads 
986       if (cluster0->GetRow() > 63) padSize = 2; // long pads
987    }
988       
989    static TLinearFitter fitY0(3, "pol2");
990    static TLinearFitter fitZ0(3, "pol2");
991    static TLinearFitter fitY2(5, "hyp4");
992    static TLinearFitter fitZ2(5, "hyp4");
993    static TLinearFitter fitY2Q(5, "hyp4");
994    static TLinearFitter fitZ2Q(5, "hyp4");
995    static TLinearFitter fitY2S(5, "hyp4");
996    static TLinearFitter fitZ2S(5, "hyp4");
997    fitY0.ClearPoints();
998    fitZ0.ClearPoints();
999    fitY2.ClearPoints();
1000    fitZ2.ClearPoints();
1001    fitY2Q.ClearPoints();
1002    fitZ2Q.ClearPoints();
1003    fitY2S.ClearPoints();
1004    fitZ2S.ClearPoints();
1005    
1006    for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
1007      // loop over irow +- kDelta rows (neighboured rows)
1008      // 
1009      // 
1010      if (idelta == 0) continue;
1011      if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
1012      AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta);
1013      if (!cluster) continue;
1014      if (cluster->GetType() < 0) continue;
1015      if (cluster->GetDetector() != sector) continue;
1016      Double_t x = cluster->GetX() - xref;
1017      Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) );
1018      Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) );
1019      //
1020      Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1021      Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1022      Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), 
1023                                                            TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1024      Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), 
1025                                                            TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1026      Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1027                                                          TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1028                                                          TMath::Sqrt(cluster0->GetSigmaY2()), 0 );
1029      Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1030                                                         TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1031                                                         TMath::Sqrt(cluster0->GetSigmaZ2()),0 );
1032      sigmaYS  = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.);
1033      sigmaZS  = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.);
1034      //
1035      if (kDelta != 0){
1036        fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
1037        fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
1038      }
1039      Double_t xxx[4];
1040      xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0;
1041      xxx[1] = x;
1042      xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0;
1043      xxx[3] = x * x;    
1044      fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
1045      fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
1046      fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
1047      fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
1048      fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
1049      fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
1050      //
1051    }  // neigbouhood-loop
1052       //
1053    fitY0.Eval();
1054    fitZ0.Eval();
1055    fitY2.Eval();
1056    fitZ2.Eval();
1057    fitY2Q.Eval();
1058    fitZ2Q.Eval();
1059    fitY2S.Eval();
1060    fitZ2S.Eval();
1061    Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.);
1062    Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.);
1063    Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.);
1064    Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.);
1065    Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.);
1066    Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.);
1067    Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.);
1068    Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.);
1069    //
1070    static  TVectorD    parY0(3);
1071    static  TMatrixD    matY0(3, 3);
1072    static  TVectorD    parZ0(3);
1073    static  TMatrixD    matZ0(3, 3);
1074    fitY0.GetParameters(parY0);
1075    fitY0.GetCovarianceMatrix(matY0);
1076    fitZ0.GetParameters(parZ0);
1077    fitZ0.GetCovarianceMatrix(matZ0);
1078    //
1079    static  TVectorD    parY2(5);
1080    static  TMatrixD    matY2(5,5);
1081    static  TVectorD    parZ2(5);
1082    static  TMatrixD    matZ2(5,5);
1083    fitY2.GetParameters(parY2);
1084    fitY2.GetCovarianceMatrix(matY2);
1085    fitZ2.GetParameters(parZ2);
1086    fitZ2.GetCovarianceMatrix(matZ2);
1087    //
1088    static  TVectorD    parY2Q(5);
1089    static  TMatrixD    matY2Q(5,5);
1090    static  TVectorD    parZ2Q(5);
1091    static  TMatrixD    matZ2Q(5,5);
1092    fitY2Q.GetParameters(parY2Q);
1093    fitY2Q.GetCovarianceMatrix(matY2Q);
1094    fitZ2Q.GetParameters(parZ2Q);
1095    fitZ2Q.GetCovarianceMatrix(matZ2Q);
1096    static  TVectorD    parY2S(5);
1097    static  TMatrixD    matY2S(5,5);
1098    static  TVectorD    parZ2S(5);
1099    static  TMatrixD    matZ2S(5,5);
1100    fitY2S.GetParameters(parY2S);
1101    fitY2S.GetCovarianceMatrix(matY2S);
1102    fitZ2S.GetParameters(parZ2S);
1103    fitZ2S.GetCovarianceMatrix(matZ2S);
1104    Float_t sigmaY0   = TMath::Sqrt(matY0(0,0));
1105    Float_t sigmaZ0   = TMath::Sqrt(matZ0(0,0));
1106    Float_t sigmaDY0  = TMath::Sqrt(matY0(1,1));
1107    Float_t sigmaDZ0  = TMath::Sqrt(matZ0(1,1));
1108    Float_t sigmaY2   = TMath::Sqrt(matY2(1,1));
1109    Float_t sigmaZ2   = TMath::Sqrt(matZ2(1,1));
1110    Float_t sigmaDY2  = TMath::Sqrt(matY2(3,3));
1111    Float_t sigmaDZ2  = TMath::Sqrt(matZ2(3,3));
1112    Float_t sigmaY2Q  = TMath::Sqrt(matY2Q(1,1));
1113    Float_t sigmaZ2Q  = TMath::Sqrt(matZ2Q(1,1));
1114    Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
1115    Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
1116    Float_t sigmaY2S  = TMath::Sqrt(matY2S(1,1));
1117    Float_t sigmaZ2S  = TMath::Sqrt(matZ2S(1,1));
1118    Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
1119    Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
1120    
1121    // Error parameters
1122    Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
1123    Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
1124    //
1125    Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1126                                                   TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1127    Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1128                                                   TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1129    Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1130                                                         TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1131    Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1132                                                         TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1133    
1134    
1135    // RMS parameters
1136    Float_t meanRMSY = 0;
1137    Float_t meanRMSZ = 0;
1138    Int_t   nclRMS = 0;
1139    for (Int_t idelta = -2; idelta <= 2; idelta++){
1140      // loop over neighbourhood
1141      if (idelta+irow < 0 || idelta+irow > 159) continue;
1142      //         if (idelta+irow>159) continue;
1143      AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
1144      if (!cluster) continue;
1145      meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
1146      meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
1147      nclRMS++;
1148    }
1149    meanRMSY /= nclRMS; 
1150    meanRMSZ /= nclRMS; 
1151    
1152    Float_t rmsY      = TMath::Sqrt(cluster0->GetSigmaY2());  
1153    Float_t rmsZ      = TMath::Sqrt(cluster0->GetSigmaZ2());
1154    Float_t rmsYT     = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1155                                               TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1156    Float_t rmsZT     = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1157                                               TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1158    Float_t rmsYT0    = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1159                                               TMath::Abs(angley));
1160    Float_t rmsZT0    = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1161                                                  TMath::Abs(anglez));
1162    Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1163                                                   TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1164    Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1165                                                   TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1166    Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1167                                                       TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1168                                                       rmsY,meanRMSY);
1169    Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1170                                                       TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1171                                                       rmsZ,meanRMSZ);
1172    //
1173    // cluster debug
1174    TTreeSRedirector *cstream = GetDebugStreamer();
1175    if (cstream){
1176      (*cstream)<<"ResolCl"<<    // valgrind 3   40,000 bytes in 1 blocks are possibly 
1177        "run="<<fRun<<              //  run number
1178        "event="<<fEvent<<          //  event number
1179        "time="<<fTime<<            //  time stamp of event
1180        "trigger="<<fTrigger<<      //  trigger
1181        "mag="<<fMagF<<             //  magnetic field         
1182        "Sector="<<sector<<
1183        "Cl.="<<cluster0<<
1184        "CSigmaY0="<<csigmaY0<<   // cluster errorY
1185        "CSigmaYQ="<<csigmaYQ<<   // cluster errorY - q scaled
1186        "CSigmaYS="<<csigmaYS<<   // cluster errorY - q scaled
1187        "CSigmaZ0="<<csigmaZ0<<   // 
1188        "CSigmaZQ="<<csigmaZQ<<
1189        "CSigmaZS="<<csigmaZS<<
1190        "shapeYF="<<rmsYFactor<<
1191        "shapeZF="<<rmsZFactor<<
1192        "rmsY="<<rmsY<<
1193        "rmsZ="<<rmsZ<<
1194        "rmsYM="<<meanRMSY<<
1195        "rmsZM="<<meanRMSZ<<
1196        "rmsYT="<<rmsYT<<
1197        "rmsZT="<<rmsZT<<
1198        "rmsYT0="<<rmsYT0<<
1199        "rmsZT0="<<rmsZT0<<
1200        "rmsYS="<<rmsYSigma<<  
1201        "rmsZS="<<rmsZSigma<<
1202        "padSize="<<padSize<<
1203        "Ncl="<<nclFound<<       
1204        "PY0.="<<&parY0<<
1205        "PZ0.="<<&parZ0<<
1206        "SigmaY0="<<sigmaY0<< 
1207        "SigmaZ0="<<sigmaZ0<< 
1208        "angley="<<angley<<
1209        "anglez="<<anglez<<
1210        "\n";
1211      
1212      //       tracklet dubug
1213      (*cstream)<<"ResolTr"<<    
1214        "run="<<fRun<<              //  run number
1215        "event="<<fEvent<<          //  event number
1216        "time="<<fTime<<            //  time stamp of event
1217        "trigger="<<fTrigger<<      //  trigger
1218        "mag="<<fMagF<<             //  magnetic field         
1219        "padSize="<<padSize<<
1220        "IPad="<<padSize<<
1221        "Sector="<<sector<<
1222        "Ncl="<<nclFound<<       
1223        "chi2Y0="<<chi2Y0<<
1224        "chi2Z0="<<chi2Z0<<
1225        "chi2Y2="<<chi2Y2<<
1226        "chi2Z2="<<chi2Z2<<
1227        "chi2Y2Q="<<chi2Y2Q<<
1228        "chi2Z2Q="<<chi2Z2Q<<
1229        "chi2Y2S="<<chi2Y2S<<
1230        "chi2Z2S="<<chi2Z2S<<
1231        "PY0.="<<&parY0<<
1232        "PZ0.="<<&parZ0<<
1233        "PY2.="<<&parY2<<
1234        "PZ2.="<<&parZ2<<
1235        "PY2Q.="<<&parY2Q<<
1236        "PZ2Q.="<<&parZ2Q<<
1237        "PY2S.="<<&parY2S<<
1238        "PZ2S.="<<&parZ2S<<
1239        "SigmaY0="<<sigmaY0<< 
1240        "SigmaZ0="<<sigmaZ0<< 
1241        "SigmaDY0="<<sigmaDY0<< 
1242        "SigmaDZ0="<<sigmaDZ0<< 
1243        "SigmaY2="<<sigmaY2<< 
1244        "SigmaZ2="<<sigmaZ2<< 
1245        "SigmaDY2="<<sigmaDY2<< 
1246        "SigmaDZ2="<<sigmaDZ2<< 
1247        "SigmaY2Q="<<sigmaY2Q<< 
1248        "SigmaZ2Q="<<sigmaZ2Q<< 
1249        "SigmaDY2Q="<<sigmaDY2Q<< 
1250        "SigmaDZ2Q="<<sigmaDZ2Q<< 
1251        "SigmaY2S="<<sigmaY2S<< 
1252        "SigmaZ2S="<<sigmaZ2S<< 
1253        "SigmaDY2S="<<sigmaDY2S<< 
1254        "SigmaDZ2S="<<sigmaDZ2S<< 
1255        "angley="<<angley<<
1256        "anglez="<<anglez<<
1257        "\n";
1258    }  
1259 }
1260
1261
1262
1263
1264
1265 TH2D * AliTPCcalibTracks::MakeDiff(TH2D * hfit, TF2 * func){
1266    // 
1267    // creates a new histogram which contains the difference between
1268    // the histogram hfit and the function func
1269    // 
1270   TH2D * result = (TH2D*)hfit->Clone();      // valgrind 3   40,139 bytes in 11 blocks are still reachable
1271   result->SetTitle(Form("%s fit residuals",result->GetTitle()));
1272   result->SetName(Form("%s fit residuals",result->GetName()));
1273   TAxis *xaxis = hfit->GetXaxis();
1274   TAxis *yaxis = hfit->GetYaxis();
1275   Double_t x[2];
1276   for (Int_t biny = 0; biny <= yaxis->GetNbins(); biny++) {
1277     x[1]  = yaxis->GetBinCenter(biny);
1278     for (Int_t binx = 0; binx <= xaxis->GetNbins(); binx++) {
1279       x[0]  = xaxis->GetBinCenter(binx);
1280       Int_t bin = hfit->GetBin(binx, biny);
1281       Double_t val = hfit->GetBinContent(bin);
1282 //      result->SetBinContent( bin, (val - func->Eval(x[0], x[1])) / func->Eval(x[0], x[1]) );
1283       result->SetBinContent( bin, (val / func->Eval(x[0], x[1])) - 1 );
1284     }    
1285   }
1286   return result;
1287 }
1288
1289
1290 void  AliTPCcalibTracks::SetStyle() const {
1291    // 
1292    // set style, can be called by all draw functions
1293    // 
1294    gROOT->SetStyle("Plain");
1295    gStyle->SetFillColor(10);
1296    gStyle->SetPadColor(10);
1297    gStyle->SetCanvasColor(10);
1298    gStyle->SetStatColor(10);
1299    gStyle->SetPalette(1,0);
1300    gStyle->SetNumberContours(60);
1301 }
1302
1303
1304 void AliTPCcalibTracks::Draw(Option_t* opt){
1305    // 
1306    // draw-function of AliTPCcalibTracks
1307    // will draws some exemplaric pictures
1308    // 
1309
1310   if (GetDebugLevel() > 6) Info("Draw", "Drawing an exemplaric picture.");
1311    SetStyle();
1312    Double_t min = 0;
1313    Double_t max = 0;
1314    TCanvas *c1 = new TCanvas();
1315    c1->Divide(0, 3);
1316    TVirtualPad *upperThird = c1->GetPad(1);
1317    TVirtualPad *middleThird = c1->GetPad(2);
1318    TVirtualPad *lowerThird = c1->GetPad(3);
1319    upperThird->Divide(2,0);
1320    TVirtualPad *upleft  = upperThird->GetPad(1);
1321    TVirtualPad *upright = upperThird->GetPad(2);
1322    middleThird->Divide(2,0);
1323    TVirtualPad *middleLeft  = middleThird->GetPad(1);
1324    TVirtualPad *middleRight = middleThird->GetPad(2);
1325    lowerThird->Divide(2,0);
1326    TVirtualPad *downLeft  = lowerThird->GetPad(1);
1327    TVirtualPad *downRight = lowerThird->GetPad(2);
1328    
1329    
1330    upleft->cd(0);
1331    min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1332    max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1333    fDeltaY->SetAxisRange(min, max);
1334    fDeltaY->Fit("gaus","q","",min, max);        // valgrind 3  7 block possibly lost   2,400 bytes in 1 blocks are still reachable
1335    c1->Update();
1336    
1337    upright->cd(0);
1338    max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1339    min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1340    fDeltaZ->SetAxisRange(min, max);
1341    fDeltaZ->Fit("gaus","q","",min, max);
1342    c1->Update();
1343    
1344    middleLeft->cd();
1345    fHclus->Draw(opt);
1346    
1347    middleRight->cd();
1348    fRejectedTracksHisto->Draw(opt);
1349    TPaveText *pt = new TPaveText(0.6,0.6, 0.8,0.8, "NDC");
1350    TText *t1 = pt->AddText("1: kEdgeThetaCutNoise");
1351    TText *t2 = pt->AddText("2: kMinClusters");
1352    TText *t3 = pt->AddText("3: kMinRatio");
1353    TText *t4 = pt->AddText("4: kMax1pt");
1354    t1 = t1; t2 = t2; t3 = t3; t4 = t4;    // avoid compiler warnings
1355    pt->SetToolTipText("Legend for failed cuts");
1356    pt->Draw();
1357    
1358    downLeft->cd();
1359    fHclusterPerPadrowRaw->Draw(opt);
1360    
1361    downRight->cd();
1362    fHclusterPerPadrow->Draw(opt);
1363 }
1364
1365
1366 void AliTPCcalibTracks::MakeReport(Int_t stat, char* pathName){ 
1367    // 
1368    // all functions are called, that produce pictures
1369    // the histograms are written to the directory 'pathName'
1370    // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
1371    // 'stat' is also the number of minEntries for MakeResPlotsQTree
1372    // 
1373
1374   if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
1375    MakeAmpPlots(stat, pathName);
1376    MakeDeltaPlots(pathName);
1377    FitResolutionNew(pathName);
1378    FitRMSNew(pathName);
1379    MakeChargeVsDriftLengthPlots(pathName);
1380 //    MakeResPlotsQ(1, 1); 
1381    MakeResPlotsQTree(stat, pathName);
1382 }
1383    
1384
1385 void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){ 
1386    // 
1387    // creates several plots:
1388    // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps
1389    // fArrayAmp.ps: one histogram per sector, the histogram shows the charge per cluster
1390    // fArrayAmpRow.ps: one histogram per sector, mean max. amplitude vs. pad row with landau fit
1391    // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1392    // Empty histograms (sectors without data) are not written to file
1393    // the ps-files are written to the directory 'pathName', that is created if it does not exist
1394    // 'stat': only histograms with more than 'stat' entries are written to file.
1395    // 
1396    
1397    SetStyle();
1398    gSystem->MakeDirectory(pathName);
1399    gSystem->ChangeDirectory(pathName);
1400    
1401    TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1402    TPostScript *ps; 
1403    // histograms with accumulated amplitude for all IROCs and OROCs
1404    TH1F *allAmpHisIROC = ((TH1F*)(fArrayAmp->At(0))->Clone());
1405    allAmpHisIROC->SetName("Amp all IROCs");
1406    allAmpHisIROC->SetTitle("Amp all IROCs");
1407    TH1F *allAmpHisOROC = ((TH1F*)(fArrayAmp->At(36))->Clone());
1408    allAmpHisOROC->SetName("Amp all OROCs");
1409    allAmpHisOROC->SetTitle("Amp all OROCs");
1410    
1411    ps = new TPostScript("fArrayAmp.ps", 112);
1412    if (GetDebugLevel() > -1) cout << "creating fArrayAmp.ps..." << endl;
1413    for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++){
1414       if ( ((TH1F*)fArrayAmp->At(i))->GetEntries() < stat  ) continue;
1415       ps->NewPage();
1416       ((TH1F*)fArrayAmp->At(i))->Draw();
1417       c1->Update();              // valgrind 3
1418       if (i > 0 && i < 36) {
1419          allAmpHisIROC->Add(((TH1F*)fArrayAmp->At(i)));
1420          allAmpHisOROC->Add(((TH1F*)fArrayAmp->At(i+36)));
1421       }
1422    }
1423    ps->NewPage();
1424    allAmpHisIROC->Draw();
1425    c1->Update();              // valgrind
1426    ps->NewPage();
1427    allAmpHisOROC->Draw();
1428    c1->Update();
1429    ps->Close();
1430    delete ps;
1431    
1432    TH1F *his = 0;
1433    Double_t min = 0;
1434    Double_t max = 0;
1435    ps = new TPostScript("fArrayAmpRow.ps", 112);
1436    if (GetDebugLevel() > -1) cout << "creating fArrayAmpRow.ps..." << endl;
1437    for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++){
1438       his = (TH1F*)fArrayAmpRow->At(i);
1439       if (his->GetEntries() < stat) continue;
1440       ps->NewPage();
1441       min = TMath::Max( his->GetBinCenter(his->GetMaximumBin() )-100., 0.);
1442       max = his->GetBinCenter(5*his->GetMaximumBin()) + 100;
1443       his->SetAxisRange(min, max);
1444       his->Fit("pol3", "q", "", min, max);
1445       // his->Draw("error");    // don't use this line when you don't want to have empty pages in the ps-file
1446       c1->Update();
1447    }
1448    ps->Close();
1449    delete ps;
1450    delete c1;
1451    gSystem->ChangeDirectory("..");
1452 }
1453
1454
1455 void AliTPCcalibTracks::MakeDeltaPlots(char* pathName){
1456    // 
1457    // creates several plots:
1458    // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1459    // the ps-files are written to the directory 'pathName', that is created if it does not exist
1460    // 
1461    
1462    SetStyle();
1463    gSystem->MakeDirectory(pathName);
1464    gSystem->ChangeDirectory(pathName);
1465    
1466    TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1467    TPostScript *ps; 
1468    Double_t min = 0;
1469    Double_t max = 0;
1470    
1471    ps = new TPostScript("DeltaYZ.ps", 112);
1472    if (GetDebugLevel() > -1) cout << "creating DeltaYZ.ps..." << endl;
1473    min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1474    max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1475    fDeltaY->SetAxisRange(min, max);
1476    ps->NewPage();
1477    fDeltaY->Fit("gaus","q","",min, max);        // valgrind 3  7 block possibly lost   2,400 bytes in 1 blocks are still reachable
1478    c1->Update();
1479    ps->NewPage();
1480    max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1481    min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1482    fDeltaZ->SetAxisRange(min, max);
1483    fDeltaZ->Fit("gaus","q","",min, max);
1484    c1->Update();
1485    ps->Close();
1486    delete ps;
1487    delete c1;
1488    gSystem->ChangeDirectory("..");
1489 }
1490
1491
1492 void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(char* pathName){
1493    // 
1494    // creates charge vs. driftlength plots, one TProfile for each ROC
1495    // is not correct like this, should be one TProfile for each sector and padsize
1496    // 
1497    
1498    SetStyle();
1499    gSystem->MakeDirectory(pathName);
1500    gSystem->ChangeDirectory(pathName);
1501    
1502    TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1503    TPostScript *ps; 
1504    ps = new TPostScript("chargeVsDriftlengthOld.ps", 112);
1505    if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlength.ps..." << endl;
1506    TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone());
1507    TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone());
1508    chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC");
1509    chargeVsDriftlengthAllIROCs->SetTitle("charge vs. driftlength, all IROCs");
1510    chargeVsDriftlengthAllOROCs->SetName("allAmpHisOROC");
1511    chargeVsDriftlengthAllOROCs->SetTitle("charge vs. driftlength, all OROCs");
1512    
1513    for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) {
1514       ((TProfile*)fArrayChargeVsDriftlength->At(i))->Draw();
1515       c1->Update();
1516       if (i > 0 && i < 36) { 
1517          chargeVsDriftlengthAllIROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i)));
1518          chargeVsDriftlengthAllOROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i+36)));
1519       }
1520       ps->NewPage();
1521    }
1522    chargeVsDriftlengthAllIROCs->Draw();
1523    c1->Update();              // valgrind
1524    ps->NewPage();
1525    chargeVsDriftlengthAllOROCs->Draw();
1526    c1->Update();
1527    ps->Close();  
1528    delete ps;
1529    delete c1;
1530    gSystem->ChangeDirectory("..");
1531 }   
1532
1533
1534 void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(char* pathName){
1535    // 
1536    // creates charge vs. driftlength plots, one TProfile for each ROC
1537    // under development....
1538    // 
1539    
1540    SetStyle();
1541    gSystem->MakeDirectory(pathName);
1542    gSystem->ChangeDirectory(pathName);
1543    
1544    TCanvas* c1 = new TCanvas("c1", "c1", 700,(Int_t)(TMath::Sqrt(2)*700));     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1545 //    TCanvas c1("c1", "c1", 500,(sqrt(2)*500))
1546    c1->Divide(0,3);
1547    TPostScript *ps; 
1548    ps = new TPostScript("chargeVsDriftlength.ps", 111);
1549    if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl;
1550    
1551    TProfile *chargeVsDriftlengthAllShortPads  = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone());
1552    TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone());
1553    TProfile *chargeVsDriftlengthAllLongPads   = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)->Clone());
1554    chargeVsDriftlengthAllShortPads->SetName("allAmpHisShortPads");
1555    chargeVsDriftlengthAllShortPads->SetTitle("charge vs. driftlength, all sectors, short pads");
1556    chargeVsDriftlengthAllMediumPads->SetName("allAmpHisMediumPads");
1557    chargeVsDriftlengthAllMediumPads->SetTitle("charge vs. driftlength, all sectors, medium pads");
1558    chargeVsDriftlengthAllLongPads->SetName("allAmpHisLongPads");
1559    chargeVsDriftlengthAllLongPads->SetTitle("charge vs. driftlength, all sectors, long pads");
1560    
1561    for (Int_t i = 0; i < 36; i++) {
1562       c1->cd(1)->SetGridx();
1563       c1->cd(1)->SetGridy();
1564       ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,0))->Draw();
1565       c1->cd(2)->SetGridx();
1566       c1->cd(2)->SetGridy();
1567       ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,1))->Draw();
1568       c1->cd(3)->SetGridx();
1569       c1->cd(3)->SetGridy();
1570       ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,2))->Draw();
1571       c1->Update();
1572       chargeVsDriftlengthAllShortPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0));
1573       chargeVsDriftlengthAllMediumPads->Add((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1));
1574       chargeVsDriftlengthAllLongPads->Add(  (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2));
1575       ps->NewPage();
1576    }
1577    c1->cd(1)->SetGridx();
1578    c1->cd(1)->SetGridy();
1579    chargeVsDriftlengthAllShortPads->Draw();
1580    c1->cd(2)->SetGridx();
1581    c1->cd(2)->SetGridy();
1582    chargeVsDriftlengthAllMediumPads->Draw();
1583    c1->cd(3)->SetGridx();
1584    c1->cd(3)->SetGridy();
1585    chargeVsDriftlengthAllLongPads->Draw();
1586    c1->Update();              // valgrind
1587 //    ps->NewPage();
1588    ps->Close();  
1589    delete ps;
1590    delete c1;
1591    gSystem->ChangeDirectory("..");
1592 }   
1593
1594
1595
1596 void AliTPCcalibTracks::FitResolutionNew(char* pathName){
1597    // 
1598    // calculates different resulution fits in Y and Z direction
1599    // the histograms are written to 'ResolutionYZ.ps'
1600    // writes calculated resolution to 'resol.txt'
1601    // all files are stored in the directory pathName
1602    // 
1603    
1604    SetStyle();
1605    gSystem->MakeDirectory(pathName);
1606    gSystem->ChangeDirectory(pathName);
1607    
1608    TCanvas c;
1609    c.Divide(2,1); 
1610    if (GetDebugLevel() > -1) cout << "creating ResolutionYZ.ps..." << endl;
1611    TPostScript *ps = new TPostScript("ResolutionYZ.ps", 112); 
1612    TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1613    fres->SetParameter(0,0.02);
1614    fres->SetParameter(1,0.0054);
1615    fres->SetParameter(2,0.13);  
1616    
1617    TH1::AddDirectory(kTRUE);  // TH3F::FitSlicesZ() writes histograms into the current directory
1618    
1619    // create histogramw for Y-resolution
1620    TH3F * hisResY0 = (TH3F*)fResolY->At(0);
1621    hisResY0->FitSlicesZ();
1622    TH2D * hisResY02 = (TH2D*)gDirectory->Get("Resol Y0_2");
1623    TH3F * hisResY1 = (TH3F*)fResolY->At(1); 
1624    hisResY1->FitSlicesZ();
1625    TH2D * hisResY12 = (TH2D*)gDirectory->Get("Resol Y1_2");
1626    TH3F * hisResY2 = (TH3F*)fResolY->At(2);
1627    hisResY2->FitSlicesZ();
1628    TH2D * hisResY22 = (TH2D*)gDirectory->Get("Resol Y2_2");
1629     //
1630    ps->NewPage();
1631    c.cd(1);
1632    hisResY02->Fit(fres, "q");      // valgrind    132,072 bytes in 6 blocks are indirectly lost
1633    hisResY02->Draw("surf1"); 
1634    c.cd(2);
1635    MakeDiff(hisResY02,fres)->Draw("surf1");
1636    c.Update();
1637    //   c.SaveAs("ResolutionYPad0.eps");
1638    ps->NewPage();
1639    c.cd(1);
1640    hisResY12->Fit(fres, "q");
1641    hisResY12->Draw("surf1");
1642    c.cd(2);
1643    MakeDiff(hisResY12,fres)->Draw("surf1");
1644    c.Update();
1645    //   c.SaveAs("ResolutionYPad1.eps");
1646    ps->NewPage();
1647    c.cd(1);
1648    hisResY22->Fit(fres, "q");
1649    hisResY22->Draw("surf1"); 
1650    c.cd(2);
1651    MakeDiff(hisResY22,fres)->Draw("surf1");
1652    c.Update();
1653 //    c.SaveAs("ResolutionYPad2.eps");
1654    
1655    // create histogramw for Z-resolution
1656    TH3F * hisResZ0 = (TH3F*)fResolZ->At(0);
1657    hisResZ0->FitSlicesZ();
1658    TH2D * hisResZ02 = (TH2D*)gDirectory->Get("Resol Z0_2");
1659    TH3F * hisResZ1 = (TH3F*)fResolZ->At(1);
1660    hisResZ1->FitSlicesZ();
1661    TH2D * hisResZ12 = (TH2D*)gDirectory->Get("Resol Z1_2");
1662    TH3F * hisResZ2 = (TH3F*)fResolZ->At(2);
1663    hisResZ2->FitSlicesZ();
1664    TH2D * hisResZ22 = (TH2D*)gDirectory->Get("Resol Z2_2");
1665    
1666    ps->NewPage();
1667    c.cd(1);
1668    hisResZ02->Fit(fres, "q");
1669    hisResZ02->Draw("surf1");
1670    c.cd(2);
1671    MakeDiff(hisResZ02,fres)->Draw("surf1");
1672    c.Update();
1673 //    c.SaveAs("ResolutionZPad0.eps");
1674    ps->NewPage();
1675    c.cd(1);
1676    hisResZ12->Fit(fres, "q");
1677    hisResZ12->Draw("surf1"); 
1678    c.cd(2);
1679    MakeDiff(hisResZ12,fres)->Draw("surf1");
1680    c.Update();
1681 //    c.SaveAs("ResolutionZPad1.eps");
1682    ps->NewPage();
1683    c.cd(1);
1684    hisResZ22->Fit(fres, "q");
1685    hisResZ22->Draw("surf1");  
1686    c.cd(2);
1687    MakeDiff(hisResZ22,fres)->Draw("surf1");
1688    c.Update();
1689 //    c.SaveAs("ResolutionZPad2.eps");
1690    ps->Close();
1691    delete ps;
1692    
1693    // write calculated resoltuions to 'resol.txt'
1694    ofstream fresol("resol.txt");
1695    fresol<<"Pad 0.75 cm"<<"\n";
1696    hisResY02->Fit(fres, "q");                     // valgrind
1697    fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1698    hisResZ02->Fit(fres, "q");
1699    fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1700    //
1701    fresol<<"Pad 1.00 cm"<<1<<"\n";
1702    hisResY12->Fit(fres, "q");                     // valgrind
1703    fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1704    hisResZ12->Fit(fres, "q");
1705    fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1706    //
1707    fresol<<"Pad 1.50 cm"<<0<<"\n";
1708    hisResY22->Fit(fres, "q");
1709    fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1710    hisResZ22->Fit(fres, "q");
1711    fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1712    
1713    TH1::AddDirectory(kFALSE);
1714    gSystem->ChangeDirectory("..");
1715    delete fres;
1716 }
1717
1718
1719 void AliTPCcalibTracks::FitRMSNew(char* pathName){
1720    // 
1721    // calculates different resulution-rms fits in Y and Z direction
1722    // the histograms are written to 'RMS_YZ.ps'
1723    // writes calculated resolution rms to 'rms.txt'
1724    // all files are stored in the directory pathName
1725    // 
1726    
1727    SetStyle();
1728    gSystem->MakeDirectory(pathName);
1729    gSystem->ChangeDirectory(pathName);
1730    
1731    TCanvas c;        // valgrind 3   42,120 bytes in 405 blocks are still reachable   23,816 bytes in 229 blocks are still reachable
1732    c.Divide(2,1); 
1733    if (GetDebugLevel() > -1) cout << "creating RMS_YZ.ps..." << endl;
1734    TPostScript *ps = new TPostScript("RMS_YZ.ps", 112); 
1735    TF2 *frms = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1736    frms->SetParameter(0,0.02);
1737    frms->SetParameter(1,0.0054);
1738    frms->SetParameter(2,0.13);  
1739    
1740    TH1::AddDirectory(kTRUE);  // TH3F::FitSlicesZ() writes histograms into the current directory
1741    
1742    // create histogramw for Y-RMS   
1743    TH3F * hisResY0 = (TH3F*)fRMSY->At(0);
1744    hisResY0->FitSlicesZ();
1745    TH2D * hisResY02 = (TH2D*)gDirectory->Get("RMS Y0_1");
1746    TH3F * hisResY1 = (TH3F*)fRMSY->At(1);
1747    hisResY1->FitSlicesZ();
1748    TH2D * hisResY12 = (TH2D*)gDirectory->Get("RMS Y1_1");
1749    TH3F * hisResY2 = (TH3F*)fRMSY->At(2);
1750    hisResY2->FitSlicesZ();
1751    TH2D * hisResY22 = (TH2D*)gDirectory->Get("RMS Y2_1");
1752    //
1753    ps->NewPage();
1754    c.cd(1);
1755    hisResY02->Fit(frms, "qn0"); 
1756    hisResY02->Draw("surf1"); 
1757    c.cd(2);
1758    MakeDiff(hisResY02,frms)->Draw("surf1");
1759    c.Update();
1760 //    c.SaveAs("RMSYPad0.eps");
1761    ps->NewPage();
1762    c.cd(1);
1763    hisResY12->Fit(frms, "qn0");               // valgrind   several blocks possibly lost
1764    hisResY12->Draw("surf1");
1765    c.cd(2);
1766    MakeDiff(hisResY12,frms)->Draw("surf1");
1767    c.Update();
1768 //    c.SaveAs("RMSYPad1.eps");
1769    ps->NewPage();
1770    c.cd(1);
1771    hisResY22->Fit(frms, "qn0");
1772    hisResY22->Draw("surf1"); 
1773    c.cd(2);
1774    MakeDiff(hisResY22,frms)->Draw("surf1");
1775    c.Update();
1776 //    c.SaveAs("RMSYPad2.eps");
1777    
1778    // create histogramw for Z-RMS   
1779    TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0);
1780    hisResZ0->FitSlicesZ();
1781    TH2D * hisResZ02 = (TH2D*)gDirectory->Get("RMS Z0_1");
1782    TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1); 
1783    hisResZ1->FitSlicesZ();
1784    TH2D * hisResZ12 = (TH2D*)gDirectory->Get("RMS Z1_1");
1785    TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2); 
1786    hisResZ2->FitSlicesZ();
1787    TH2D * hisResZ22 = (TH2D*)gDirectory->Get("RMS Z2_1");
1788    //
1789    ps->NewPage();
1790    c.cd(1);
1791    hisResZ02->Fit(frms, "qn0");         // valgrind
1792    hisResZ02->Draw("surf1");
1793    c.cd(2);
1794    MakeDiff(hisResZ02,frms)->Draw("surf1");
1795    c.Update();
1796 //    c.SaveAs("RMSZPad0.eps");
1797    ps->NewPage();
1798    c.cd(1);
1799    hisResZ12->Fit(frms, "qn0");
1800    hisResZ12->Draw("surf1"); 
1801    c.cd(2);
1802    MakeDiff(hisResZ12,frms)->Draw("surf1");
1803    c.Update();
1804 //    c.SaveAs("RMSZPad1.eps");
1805    ps->NewPage();
1806    c.cd(1);
1807    hisResZ22->Fit(frms, "qn0");         // valgrind  1 block possibly lost
1808    hisResZ22->Draw("surf1");  
1809    c.cd(2);
1810    MakeDiff(hisResZ22,frms)->Draw("surf1");
1811    c.Update();
1812 //    c.SaveAs("RMSZPad2.eps");
1813    
1814    // write calculated resoltuion rms to 'rms.txt'
1815    ofstream filerms("rms.txt");
1816    filerms<<"Pad 0.75 cm"<<"\n";
1817    hisResY02->Fit(frms, "qn0");         // valgrind   23 blocks indirectly lost
1818    filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1819    hisResZ02->Fit(frms, "qn0");         // valgrind   23 blocks indirectly lost
1820    filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1821    //
1822    filerms<<"Pad 1.00 cm"<<1<<"\n";
1823    hisResY12->Fit(frms, "qn0");         // valgrind      3,256 bytes in 22 blocks are indirectly lost 
1824    filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1825    hisResZ12->Fit(frms, "qn0");         // valgrind    66,036 bytes in 3 blocks are still reachable
1826    filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1827    //
1828    filerms<<"Pad 1.50 cm"<<0<<"\n";
1829    hisResY22->Fit(frms, "qn0");      // valgrind   40,139 bytes in 11 blocks are still reachable   330,180 bytes in 15 blocks are possibly lost
1830    filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1831    hisResZ22->Fit(frms, "qn0");
1832    filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1833
1834    TH1::AddDirectory(kFALSE);
1835    gSystem->ChangeDirectory("..");
1836    ps->Close();
1837    delete ps;
1838    delete frms;
1839 }
1840
1841
1842 void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){
1843    //
1844    // Make tree with resolution parameters
1845    // the result is written to 'resol.root' in directory 'pathname'
1846    // file information are available in fileInfo
1847    // available variables in the tree 'Resol':
1848    //  Entries: number of entries for this resolution point
1849    //  nbins:   number of bins that were accumulated
1850    //  Dim:     direction, Dim==0: y-direction, Dim==1: z-direction
1851    //  Pad:     padSize; short, medium and long
1852    //  Length:  pad length, 0.75, 1, 1.5
1853    //  QMean:   mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1854    //  Zc:      center of middle bin in drift direction
1855    //  Zm:      mean dirftlength for accumulated Delta-Histograms
1856    //  Zs:      width of driftlength bin
1857    //  AngleC:  center of middle bin in Angle-Direction
1858    //  AngleM:  mean angle for accumulated Delta-Histograms
1859    //  AngleS:  width of Angle-bin
1860    //  Resol:   sigma for gaus fit through Delta-Histograms
1861    //  Sigma:   error of sigma for gaus fit through Delta Histograms
1862    //  MeanR:   mean of the Delta-Histogram
1863    //  SigmaR:  rms of the Delta-Histogram
1864    //  RMSm:    mean of the gaus fit through RMS-Histogram
1865    //  RMS:     sigma of the gaus fit through RMS-Histogram
1866    //  RMSe0:   error of mean of gaus fit in RMS-Histogram
1867    //  RMSe1:   error of sigma of gaus fit in RMS-Histogram
1868    //  
1869       
1870   if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
1871   if (GetDebugLevel() > -1) cout << "    relax, the calculation will take a while..." << endl;
1872   
1873    gSystem->MakeDirectory(pathName);
1874    gSystem->ChangeDirectory(pathName);
1875    TString kFileName = "resol.root";
1876    TTreeSRedirector fTreeResol(kFileName.Data());
1877    
1878    TH3F *resArray[2][3][11];
1879    TH3F *rmsArray[2][3][11];
1880   
1881    // load histograms from fArrayQDY and fArrayQDZ 
1882    // into resArray and rmsArray
1883    // that is all we need here
1884    for (Int_t idim = 0; idim < 2; idim++){
1885       for (Int_t ipad = 0; ipad < 3; ipad++){
1886          for (Int_t iq = 0; iq <= 10; iq++){
1887             rmsArray[idim][ipad][iq]=0;
1888             resArray[idim][ipad][iq]=0;
1889             Int_t bin = GetBin(iq,ipad); 
1890             TH3F *hresl = 0;
1891             if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
1892             if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
1893             if (!hresl) continue;
1894             resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
1895             resArray[idim][ipad][iq]->SetDirectory(0);
1896             TH3F * hreslRMS = 0;
1897             if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
1898             if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
1899             if (!hreslRMS) continue;
1900             rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
1901             rmsArray[idim][ipad][iq]->SetDirectory(0);
1902          }
1903       }
1904    }
1905    if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
1906    
1907    //--------------------------------------------------------------------------------------------
1908    
1909    char name[200];
1910    Double_t qMean;
1911    Double_t zMean, angleMean, zCenter, angleCenter;
1912    Double_t zSigma, angleSigma;
1913    TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
1914    TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
1915    TF1 *fitFunction = new TF1("fitFunction", "gaus");
1916    Float_t entriesQ = 0;
1917    Int_t loopCounter = 1;
1918   
1919    for (Int_t idim = 0; idim < 2; idim++){
1920       // Loop y-z corrdinate
1921       for (Int_t ipad = 0; ipad < 3; ipad++){
1922          // loop pad type
1923          for (Int_t iq = -1; iq < 10; iq++){
1924             // LOOP Q
1925            if (GetDebugLevel() > -1) 
1926                cout << "Loop-counter, this is loop " << loopCounter << " of 66, (" 
1927                   << (Int_t)((loopCounter)/66.*100) << "% done), " 
1928                   << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << "  \r" << std::flush;
1929             loopCounter++;
1930             TH3F *hres = 0;
1931             TH3F *hrms = 0;
1932             qMean    = 0;
1933             entriesQ = 0;
1934             
1935             // calculate qMean
1936             if (iq == -1){
1937                // integrated spectra
1938                for (Int_t iql = 0; iql < 10; iql++){    
1939                   Int_t bin = GetBin(iql,ipad); 
1940                   TH3F *hresl = resArray[idim][ipad][iql];
1941                   TH3F *hrmsl = rmsArray[idim][ipad][iql];
1942                   if (!hresl) continue;
1943                   if (!hrmsl) continue;     
1944                   entriesQ += hresl->GetEntries();
1945                   qMean += hresl->GetEntries() * GetQ(bin);      
1946                   if (!hres) {
1947                      hres = (TH3F*)hresl->Clone();
1948                      hrms = (TH3F*)hrmsl->Clone();
1949                   }
1950                   else{
1951                      hres->Add(hresl);
1952                      hrms->Add(hrmsl);
1953                   }
1954                }
1955                qMean /= entriesQ;
1956                qMean *= -1.;  // integral mean charge
1957             }
1958             else {
1959                // loop over neighboured Q-bins 
1960                // accumulate entries from neighboured Q-bins
1961                for (Int_t iql = iq - 1; iql <= iq + 1; iql++){              
1962                   if (iql < 0) continue;
1963                   Int_t bin = GetBin(iql,ipad);
1964                   TH3F * hresl = resArray[idim][ipad][iql];
1965                   TH3F * hrmsl = rmsArray[idim][ipad][iql];
1966                   if (!hresl) continue;
1967                   if (!hrmsl) continue;
1968                   entriesQ += hresl->GetEntries(); 
1969                   qMean += hresl->GetEntries() * GetQ(bin);      
1970                   if (!hres) {
1971                      hres = (TH3F*) hresl->Clone();
1972                      hrms = (TH3F*) hrmsl->Clone();
1973                   }
1974                   else{
1975                      hres->Add(hresl);
1976                      hrms->Add(hrmsl);
1977                   }
1978                }
1979                qMean/=entriesQ;
1980             }
1981       
1982             TAxis *xAxisDriftLength = hres->GetXaxis();   // driftlength / z - axis
1983             TAxis *yAxisAngle       = hres->GetYaxis();   // angle axis
1984             TAxis *zAxisDelta       = hres->GetZaxis();   // delta axis
1985             TAxis *zAxisRms         = hrms->GetZaxis();   // rms axis
1986             
1987             // loop over all angle bins
1988             for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
1989                angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
1990                // loop over all driftlength bins
1991                for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
1992                   zCenter    = xAxisDriftLength->GetBinCenter(ibinxDL);
1993                   zSigma     = xAxisDriftLength->GetBinWidth(ibinxDL);
1994                   angleSigma = yAxisAngle->GetBinWidth(ibinyAngle); 
1995                   zMean      = zCenter;      // changens, when more statistic is accumulated
1996                   angleMean  = angleCenter;  // changens, when more statistic is accumulated
1997                   
1998                   // create 2 1D-Histograms, projectionRes and projectionRms
1999                   // these histograms are delta histograms for given direction, padSize, chargeBin,
2000                   // angleBin and driftLengthBin
2001                   // later on they will be fitted with a gausian, its sigma is the resoltuion...
2002                   sprintf(name,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
2003                   // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
2004                   projectionRes->SetNameTitle(name, name);
2005                   sprintf(name,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
2006                   // TH1D * projectionRms =  new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
2007                   projectionRms->SetNameTitle(name, name);
2008                   
2009                   projectionRes->Reset();
2010                   projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
2011                   projectionRms->Reset();
2012                   projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
2013                   projectionRes->SetDirectory(0);
2014                   projectionRms->SetDirectory(0);
2015                   
2016                   Double_t entries = 0;
2017                   Int_t    nbins   = 0;   // counts, how many bins were accumulated
2018                   zMean     = 0;
2019                   angleMean = 0;
2020                   entries   = 0;
2021                   
2022                   // fill projectionRes and projectionRms for given dim, ipad and iq, 
2023                   // as well as for given angleBin and driftlengthBin
2024                   // if this gives not enough statistic, include neighbourhood 
2025                   // (angle and driftlength) successifely
2026                   for (Int_t dbin = 0; dbin <= 8; dbin++){              // delta-bins around centered angleBin and driftlengthBin
2027                      for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) {   // delta-bins in angle direction
2028                         for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
2029                            if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue;   // add each bin only one time !
2030                            Int_t binx2 = ibinxDL + dbinx2;                       // position variable in x (driftlength) direction
2031                            Int_t biny2 = ibinyAngle + dbiny2;                    // position variable in y (angle)  direction
2032                            if (binx2 < 1 || biny2 < 1) continue;                 // don't go out of the histogram!
2033                            if (binx2 >= xAxisDriftLength->GetNbins()) continue;  // don't go out of the histogram!
2034                            if (biny2 >= yAxisAngle->GetNbins()) continue;        // don't go out of the histogram!
2035                            nbins++;                                              // count the number of accumulated bins
2036                            // Fill resolution histo
2037                            for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
2038                               // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3);     // unused variable
2039                               projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
2040                               entries   += hres->GetBinContent(binx2, biny2, ibin3);
2041                               zMean     += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
2042                               angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
2043                            }  // ibin3 loop
2044                            // fill RMS histo
2045                            for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
2046                               projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
2047                            }
2048                         }  //dbinx2 loop
2049                         if (entries > minEntries) break; // enough statistic accumulated
2050                      }  // dbiny2 loop
2051                      if (entries > minEntries) break;    // enough statistic accumulated
2052                   }  // dbin loop
2053                   if ( entries< minEntries) continue;  // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree  
2054                   zMean /= entries;
2055                   angleMean /= entries;
2056                   
2057                   if (entries > minEntries) {
2058                      //  when enough statistic is accumulated
2059                      //  fit Delta histograms with a gausian
2060                      //  of the gausian is the resolution (resol), its fit error is sigma
2061                      //  store also mean and RMS of the histogram
2062                      Float_t xmin     = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2063                      Float_t xmax     = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2064                      
2065 //                      projectionRes->Fit("gaus", "q0", "", xmin, xmax);
2066 //                      Float_t resol    = projectionRes->GetFunction("gaus")->GetParameter(2);
2067 //                      Float_t sigma    = projectionRes->GetFunction("gaus")->GetParError(2);
2068                      fitFunction->SetMaximum(xmax);
2069                      fitFunction->SetMinimum(xmin);
2070                      projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
2071                      Float_t resol    = fitFunction->GetParameter(2);
2072                      Float_t sigma    = fitFunction->GetParError(2);
2073                      
2074                      Float_t meanR    = projectionRes->GetMean();
2075                      Float_t sigmaR   = projectionRes->GetRMS();
2076                      // fit also RMS histograms with a gausian
2077                      // store mean and sigma of the gausian in rmsMean and rmsSigma
2078                      // store also the fit errors in errorRMS and errorSigma
2079                      xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2080                      xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2081                      
2082 //                      projectionRms->Fit("gaus","q0","",xmin,xmax);
2083 //                      Float_t rmsMean    = projectionRms->GetFunction("gaus")->GetParameter(1);
2084 //                      Float_t rmsSigma   = projectionRms->GetFunction("gaus")->GetParameter(2);
2085 //                      Float_t errorRMS   = projectionRms->GetFunction("gaus")->GetParError(1);
2086 //                      Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
2087                      projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
2088                      Float_t rmsMean    = fitFunction->GetParameter(1);
2089                      Float_t rmsSigma   = fitFunction->GetParameter(2);
2090                      Float_t errorRMS   = fitFunction->GetParError(1);
2091                      Float_t errorSigma = fitFunction->GetParError(2);
2092                     
2093                      Float_t length = 0.75;
2094                      if (ipad == 1) length = 1;
2095                      if (ipad == 2) length = 1.5;
2096                      
2097                      fTreeResol<<"Resol"<<
2098                         "Entries="<<entries<<      // number of entries for this resolution point
2099                         "nbins="<<nbins<<          // number of bins that were accumulated
2100                         "Dim="<<idim<<             // direction, Dim==0: y-direction, Dim==1: z-direction
2101                         "Pad="<<ipad<<             // padSize; short, medium and long
2102                         "Length="<<length<<        // pad length, 0.75, 1, 1.5
2103                         "QMean="<<qMean<<          // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
2104                         "Zc="<<zCenter<<           // center of middle bin in drift direction
2105                         "Zm="<<zMean<<             // mean dirftlength for accumulated Delta-Histograms
2106                         "Zs="<<zSigma<<            // width of driftlength bin
2107                         "AngleC="<<angleCenter<<   // center of middle bin in Angle-Direction
2108                         "AngleM="<<angleMean<<     // mean angle for accumulated Delta-Histograms
2109                         "AngleS="<<angleSigma<<    // width of Angle-bin
2110                         "Resol="<<resol<<          // sigma for gaus fit through Delta-Histograms
2111                         "Sigma="<<sigma<<          // error of sigma for gaus fit through Delta Histograms
2112                         "MeanR="<<meanR<<          // mean of the Delta-Histogram
2113                         "SigmaR="<<sigmaR<<        // rms of the Delta-Histogram
2114                         "RMSm="<<rmsMean<<         // mean of the gaus fit through RMS-Histogram
2115                         "RMSs="<<rmsSigma<<        // sigma of the gaus fit through RMS-Histogram
2116                         "RMSe0="<<errorRMS<<       // error of mean of gaus fit in RMS-Histogram
2117                         "RMSe1="<<errorSigma<<     // error of sigma of gaus fit in RMS-Histogram
2118                         "\n";
2119                      if (GetDebugLevel() > 5) {
2120                         projectionRes->SetDirectory(fTreeResol.GetFile());
2121                         projectionRes->Write(projectionRes->GetName());
2122                         projectionRes->SetDirectory(0);
2123                         projectionRms->SetDirectory(fTreeResol.GetFile());
2124                         projectionRms->Write(projectionRms->GetName());
2125                         projectionRes->SetDirectory(0);
2126                      }
2127                   }  // if (projectionRes->GetSum() > minEntries)
2128                }  // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
2129             }  // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
2130             
2131          }  // iq-loop
2132       }  // ipad-loop
2133    }  // idim-loop
2134    delete projectionRes;
2135    delete projectionRms;
2136    
2137 //    TFile resolFile(fTreeResol.GetFile());
2138    TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
2139    fileInfo.Write("fileInfo");
2140 //    resolFile.Close();
2141 //    fTreeResol.GetFile()->Close();
2142    if (GetDebugLevel() > -1) cout << endl;
2143    if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
2144    gSystem->ChangeDirectory("..");
2145 }
2146
2147
2148
2149
2150
2151 Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
2152    // 
2153    // function to merge several AliTPCcalibTracks objects after PROOF calculation
2154    // The object's histograms are merged via their merge functions
2155    // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
2156    // 
2157    
2158   if (GetDebugLevel() > 0) cout << " *****  this is AliTPCcalibTracks::Merge(TCollection *collectionList)  *****"<< endl;  
2159    if (!collectionList) return 0;
2160    if (collectionList->IsEmpty()) return -1;
2161    
2162    if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl;     //    REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
2163    if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
2164    collectionList->Print();
2165    
2166    // create a list for each data member
2167    TList* deltaYList = new TList;
2168    TList* deltaZList = new TList;
2169    TList* arrayAmpRowList = new TList;
2170    TList* rejectedTracksList = new TList;
2171    TList* hclusList = new TList;
2172    TList* clusterCutHistoList = new TList;
2173    TList* arrayAmpList = new TList;
2174    TList* arrayQDYList = new TList;
2175    TList* arrayQDZList = new TList;
2176    TList* arrayQRMSYList = new TList;
2177    TList* arrayQRMSZList = new TList;
2178    TList* arrayChargeVsDriftlengthList = new TList;
2179    TList* calPadRegionChargeVsDriftlengthList = new TList;
2180    TList* hclusterPerPadrowList = new TList;
2181    TList* hclusterPerPadrowRawList = new TList;
2182    TList* resolYList = new TList;
2183    TList* resolZList = new TList;
2184    TList* rMSYList = new TList;
2185    TList* rMSZList = new TList;
2186    
2187 //    TList* nRowsList = new TList;
2188 //    TList* nSectList = new TList;
2189 //    TList* fileNoList = new TList;
2190    
2191    TIterator *listIterator = collectionList->MakeIterator();
2192    AliTPCcalibTracks *calibTracks = 0;
2193    if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;    
2194    Int_t counter = 0;
2195    while ( (calibTracks = (AliTPCcalibTracks*)listIterator->Next()) ){
2196       // loop over all entries in the collectionList and get dataMembers into lists
2197       if (!calibTracks) continue;
2198       
2199       deltaYList->Add( calibTracks->GetfDeltaY() );
2200       deltaZList->Add( calibTracks->GetfDeltaZ() );
2201       arrayAmpRowList->Add(calibTracks->GetfArrayAmpRow());
2202       arrayAmpList->Add(calibTracks->GetfArrayAmp());
2203       arrayQDYList->Add(calibTracks->GetfArrayQDY());
2204       arrayQDZList->Add(calibTracks->GetfArrayQDZ());
2205       arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
2206       arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
2207       resolYList->Add(calibTracks->GetfResolY());
2208       resolZList->Add(calibTracks->GetfResolZ());
2209       rMSYList->Add(calibTracks->GetfRMSY());
2210       rMSZList->Add(calibTracks->GetfRMSZ());
2211       arrayChargeVsDriftlengthList->Add(calibTracks->GetfArrayChargeVsDriftlength());
2212       calPadRegionChargeVsDriftlengthList->Add(calibTracks->GetCalPadRegionchargeVsDriftlength());
2213       hclusList->Add(calibTracks->GetfHclus());
2214       rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
2215       clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
2216       hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow());
2217       hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw());
2218       fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
2219       fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
2220       counter++;
2221       if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
2222    }
2223    
2224    
2225    // merge data members
2226    if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl; 
2227    fDeltaY->Merge(deltaYList);
2228    fDeltaZ->Merge(deltaZList);
2229    fHclus->Merge(hclusList);
2230    fClusterCutHisto->Merge(clusterCutHistoList);
2231    fRejectedTracksHisto->Merge(rejectedTracksList);
2232    fHclusterPerPadrow->Merge(hclusterPerPadrowList);
2233    fHclusterPerPadrowRaw->Merge(hclusterPerPadrowRawList);
2234    
2235    TObjArray* objarray = 0;
2236    TH1* hist = 0;
2237    TList* histList = 0;
2238    TIterator *objListIterator = 0;
2239    
2240    if (GetDebugLevel() > 0) cout << "merging fArrayAmpRows..." << endl;
2241    // merge fArrayAmpRows
2242    for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) {  // loop over data member, i<72
2243       objListIterator = arrayAmpRowList->MakeIterator();
2244       histList = new TList;
2245       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2246          // loop over arrayAmpRowList, get TObjArray, get object at position i, cast it into TProfile
2247          if (!objarray) continue;
2248          hist = (TProfile*)(objarray->At(i));
2249          histList->Add(hist);
2250       }
2251       ((TProfile*)(fArrayAmpRow->At(i)))->Merge(histList);
2252       delete histList;
2253       delete objListIterator;
2254    }
2255    
2256    if (GetDebugLevel() > 0) cout << "merging fArrayAmps..." << endl;
2257    // merge fArrayAmps
2258    for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) {  // loop over data member, i<72
2259       TIterator *objListIterator = arrayAmpList->MakeIterator();
2260       histList = new TList;
2261       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2262          // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F
2263          if (!objarray) continue;
2264          hist = (TH1F*)(objarray->At(i));
2265          histList->Add(hist);
2266       }
2267       ((TH1F*)(fArrayAmp->At(i)))->Merge(histList);
2268       delete histList;
2269       delete objListIterator;
2270    }
2271    
2272    if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
2273    // merge fArrayQDY
2274    for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
2275       objListIterator = arrayQDYList->MakeIterator();
2276       histList = new TList;
2277       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2278          // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
2279          if (!objarray) continue;
2280          hist = (TH3F*)(objarray->At(i));
2281          histList->Add(hist);
2282       }
2283       ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
2284       delete histList;
2285       delete objListIterator;
2286    }
2287
2288    if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
2289    // merge fArrayQDZ
2290    for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2291       objListIterator = arrayQDZList->MakeIterator();
2292       histList = new TList;
2293       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2294          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2295          if (!objarray) continue;
2296          hist = (TH3F*)(objarray->At(i));
2297          histList->Add(hist);
2298       }
2299       ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
2300       delete histList;
2301       delete objListIterator;
2302    }
2303
2304    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
2305    // merge fArrayQRMSY
2306    for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
2307       objListIterator = arrayQRMSYList->MakeIterator();
2308       histList = new TList;
2309       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2310          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2311          if (!objarray) continue;
2312          hist = (TH3F*)(objarray->At(i));
2313          histList->Add(hist);
2314       }
2315       ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
2316       delete histList;
2317       delete objListIterator;
2318    }   
2319
2320    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
2321    // merge fArrayQRMSZ
2322    for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2323       objListIterator = arrayQRMSZList->MakeIterator();
2324       histList = new TList;
2325       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2326          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2327          if (!objarray) continue;
2328          hist = (TH3F*)(objarray->At(i));
2329          histList->Add(hist);
2330       }
2331       ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
2332       delete histList;
2333       delete objListIterator;
2334    } 
2335    
2336    if (GetDebugLevel() > 0) cout << "merging fArrayChargeVsDriftlength..." << endl;
2337    // merge fArrayChargeVsDriftlength
2338    for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { // loop over data member, i < 300
2339       objListIterator = arrayChargeVsDriftlengthList->MakeIterator();
2340       histList = new TList;
2341       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2342          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TProfile
2343          if (!objarray) continue;
2344          hist = (TProfile*)(objarray->At(i));
2345          histList->Add(hist);
2346       }
2347       ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Merge(histList);
2348       delete histList;
2349       delete objListIterator;
2350    }    
2351    
2352    if (GetDebugLevel() > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl;
2353    // merge fcalPadRegionChargeVsDriftlength
2354    AliTPCCalPadRegion *cpr = 0x0;
2355    
2356    /*
2357    TIterator *regionIterator = fcalPadRegionChargeVsDriftlength->MakeIterator();
2358    while (hist = (TProfile*)regionIterator->Next()) {
2359       // loop over all calPadRegion's in destination calibTracks object
2360          objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2361          while (( cpr =  (AliTPCCalPadRegion*)objListIterator->Next() )) { 
2362
2363       
2364       hist->Merge(...);
2365    }
2366    */
2367    
2368    for (UInt_t isec = 0; isec < 36; isec++) {
2369       for (UInt_t padSize = 0; padSize < 3; padSize++){
2370          objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2371          histList = new TList;
2372          while (( cpr =  (AliTPCCalPadRegion*)objListIterator->Next() )) { 
2373             // loop over calPadRegionChargeVsDriftlengthList, get AliTPCCalPadRegion, get object 
2374             if (!cpr) continue;
2375             hist = (TProfile*)cpr->GetObject(isec, padSize);
2376             histList->Add(hist);
2377          }
2378          ((TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isec, padSize)))->Merge(histList);
2379          delete histList;
2380          delete objListIterator;
2381       }
2382    }
2383   
2384    
2385         
2386    
2387    if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
2388    // merge fResolY
2389    for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
2390       objListIterator = resolYList->MakeIterator();
2391       histList = new TList;
2392       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2393          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2394          if (!objarray) continue;
2395          hist = (TH3F*)(objarray->At(i));
2396          histList->Add(hist);
2397       }
2398       ((TH3F*)(fResolY->At(i)))->Merge(histList);
2399       delete histList;
2400       delete objListIterator;
2401    }
2402    
2403    // merge fResolZ
2404    for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2405       objListIterator = resolZList->MakeIterator();
2406       histList = new TList;
2407       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2408          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2409          if (!objarray) continue;
2410          hist = (TH3F*)(objarray->At(i));
2411          histList->Add(hist);
2412       }
2413       ((TH3F*)(fResolZ->At(i)))->Merge(histList);
2414       delete histList;
2415       delete objListIterator;
2416    }
2417
2418    // merge fRMSY
2419    for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
2420       objListIterator = rMSYList->MakeIterator();
2421       histList = new TList;
2422       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2423          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2424          if (!objarray) continue;
2425          hist = (TH3F*)(objarray->At(i));
2426          histList->Add(hist);
2427       }
2428       ((TH3F*)(fRMSY->At(i)))->Merge(histList);
2429       delete histList;
2430       delete objListIterator;
2431    }
2432          
2433    // merge fRMSZ
2434    for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2435       objListIterator = rMSZList->MakeIterator();
2436       histList = new TList;
2437       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2438          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2439          if (!objarray) continue;
2440          hist = (TH3F*)(objarray->At(i));
2441          histList->Add(hist);
2442       }
2443       ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
2444       delete histList;
2445       delete objListIterator;
2446    }
2447    
2448    delete deltaYList;
2449    delete deltaZList;
2450    delete arrayAmpRowList;
2451    delete arrayAmpList;
2452    delete arrayQDYList;
2453    delete arrayQDZList;
2454    delete arrayQRMSYList;
2455    delete arrayQRMSZList;
2456    delete resolYList;
2457    delete resolZList;
2458    delete rMSYList;
2459    delete rMSZList;
2460    delete listIterator;
2461    
2462    if (GetDebugLevel() > 0) cout << "merging done!" << endl;
2463    
2464    return 1;
2465 }
2466
2467
2468 AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks){
2469    // 
2470    // function to test AliTPCcalibTrack::Merge:
2471    // in the file 'f' is a AliTPCcalibTrack object with name "calibTracks"
2472    // this object is appended 'nCalTracks' times to a TList
2473    // A new AliTPCcalibTrack object is created which merges the list
2474    // this object is returned
2475    // 
2476    /*
2477    // .L AliTPCcalibTracks.cxx+g
2478    .L libTPCcalib.so
2479    TFile f("Output.root");
2480    AliTPCcalibTracks* calTracks = (AliTPCcalibTracks*)f.Get("calibTracks");
2481    //f.Close();
2482    TFile clusterParamFile("/u/lbozyk/calibration/workdir/calibTracks/TPCClusterParam.root");
2483    AliTPCClusterParam *clusterParam  =  (AliTPCClusterParam *) clusterParamFile.Get("Param"); 
2484    clusterParamFile.Close();
2485
2486    AliTPCcalibTracks::TestMerge(calTracks, clusterParam);
2487    */
2488    TList *list = new TList();
2489    if (ct == 0 || clusterParam == 0) return 0;
2490    cout << "making list with " << nCalTracks << " AliTPCcalibTrack objects" << endl;
2491    for (Int_t i = 0; i < nCalTracks; i++) {
2492       if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl;
2493       list->Add(new AliTPCcalibTracks(*ct));
2494    }
2495    
2496    // only for check at the end
2497    AliTPCcalibTracks* cal1 = new AliTPCcalibTracks(*ct);
2498    Double_t cal1Entries = ((TH1F*)cal1->GetfArrayAmpRow()->At(5))->GetEntries();
2499 //    Double_t cal1Entries = 5; //((TH1F*)ct->GetfArrayAmpRow()->At(5))->GetEntries();
2500
2501    cout  << "The list contains " << list->GetEntries() << " entries. " << endl;
2502   
2503    
2504    AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018);
2505    AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts, 5);
2506    cal->Merge(list);
2507    
2508    cout << "cal->GetfArrayAmpRow()->At(5)->Print():" << endl;
2509    cal->GetfArrayAmpRow()->At(5)->Print();
2510    Double_t calEntries = ((TH1F*)cal->GetfArrayAmpRow()->At(5))->GetEntries();
2511    
2512    cout << "cal1->GetfArrayAmpRow()->At(5))->GetEntries() = " << cal1Entries << endl;
2513    cout << " cal->GetfArrayAmpRow()->At(5))->GetEntries() = " <<  calEntries << endl;
2514    printf("That means there were %f / %f = %f AliTPCcalibTracks-Objects merged. \n", 
2515       calEntries, cal1Entries, ((Double_t)calEntries/cal1Entries));
2516    
2517    return cal;
2518
2519 }
2520
2521
2522 void  AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints, Int_t verbose){
2523   //
2524   // Make position corrections
2525   // for the moment Only using debug streamer 
2526   // chainres  - debug tree
2527   // param     - parameters to be updated
2528   // maxPoints - maximal number of points using for fit
2529   // verbose   - print info flag
2530   //
2531   // Current implementation - only using debug streamers
2532   // 
2533   
2534   /*    
2535     //Defaults
2536     Int_t maxPoints=100000;
2537   */
2538   /*
2539     Usage: 
2540     //0. Load libraries
2541     gSystem->Load("libANALYSIS");
2542     gSystem->Load("libSTAT");
2543     gSystem->Load("libTPCcalib");
2544     
2545
2546     //1. Load Parameters to be modified:
2547     //e.g:
2548     
2549     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
2550     AliCDBManager::Instance()->SetRun(0) 
2551     AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
2552
2553     //2. Load chain from debug streamers
2554     //
2555     //e.g
2556     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2557     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2558     AliXRDPROOFtoolkit tool;  
2559     TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
2560     chainres->Lookup();
2561     //3. Do fits and store results
2562     // 
2563     AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ;
2564     TFile f("paramout.root","recreate");
2565     param->Write("clusterParam");
2566     f.Close();
2567     //4. Verification
2568     TFile f2("paramout.root");
2569     AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
2570     param2->SetInstance(param2);
2571     chainres->Draw("fitZ0:AliTPCClusterParam::SPosCorrection(1,0,Cl.fPad,Cl.fTimeBin,Cl.fZ,Cl.fSigmaY2,Cl.fSigmaZ2,Cl.fMax)","Cl.fDetector<36","",10000); // should be line 
2572     
2573    */
2574
2575
2576   TStatToolkit toolkit;
2577   Double_t chi2;
2578   TVectorD fitParamY0;
2579   TVectorD fitParamY1;
2580   TVectorD fitParamZ0;
2581   TVectorD fitParamZ1;
2582   TMatrixD covMatrix;
2583   Int_t npoints;
2584   
2585   chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
2586   chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
2587   chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
2588   chainres->SetAlias("st","(sin(dt)-dt)");
2589   //
2590   chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
2591   //
2592   //
2593   //
2594   TCut cutA("1");
2595   TString fstringY="";  
2596   //
2597   fstringY+="(dp)++";            //1
2598   fstringY+="(dp)*di++";         //2
2599   fstringY+="(sp)++";            //3
2600   fstringY+="(sp)*di++";         //4
2601   TString fstringZ="";  
2602   fstringZ+="(dt)++";            //1
2603   fstringZ+="(dt)*di++";         //2
2604   fstringZ+="(st)++";            //3
2605   fstringZ+="(st)*di++";         //4
2606   //
2607   // Z corrections
2608   //
2609   TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
2610   printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2611   param->fPosZcor[0] = (TVectorD*) fitParamZ0.Clone();
2612   //
2613   TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
2614   //
2615   printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2616   param->fPosZcor[1] = (TVectorD*) fitParamZ1.Clone();
2617   param->fPosZcor[2] = (TVectorD*) fitParamZ1.Clone();
2618   //
2619   // Y corrections
2620   //   
2621   TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
2622   printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2623   param->fPosYcor[0] = (TVectorD*) fitParamY0.Clone();
2624   
2625
2626   TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
2627   //
2628   printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2629   param->fPosYcor[1] = (TVectorD*) fitParamY1.Clone();
2630   param->fPosYcor[2] = (TVectorD*) fitParamY1.Clone();
2631   //
2632   //
2633   //
2634   chainres->SetAlias("fitZ0",strZ0->Data());
2635   chainres->SetAlias("fitZ1",strZ1->Data());
2636   chainres->SetAlias("fitY0",strY0->Data());
2637   chainres->SetAlias("fitY1",strY1->Data());
2638   //  chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000);   
2639 }
2640
2641
2642