Adding charge and position corrections as a function of relative position to the...
[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           "padSize="<<padSize<<
932           "angley="<<angley<<
933           "anglez="<<anglez<<
934           "zdr="<<zdrift<<
935           "dy="<<deltay<<
936           "dz="<<deltaz<<
937           "sy="<<sy<<
938           "sz="<<sz<<
939           "\n";
940       }
941
942       if (useForResol){
943          fDeltaY->Fill(deltay);
944          fDeltaZ->Fill(deltaz);
945          his3 = (TH3F*)fResolY->At(padSize);
946          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
947          his3 = (TH3F*)fResolZ->At(padSize);
948          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
949          his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
950          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
951          his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
952          if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
953       }
954       
955       //=============================================================================================
956       
957       if (useForResol && nclFound > 2 * kMinRatio * kDelta 
958           && irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){
959         if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
960          FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta);
961       }  // if (useForResol && nclFound > 2 * kMinRatio * kDelta)
962    
963    }    // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
964 }  // FillResolutionHistoLocal(...)
965
966
967
968 void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t  angley, Float_t  anglez, Int_t nclFound, Int_t kDelta) {
969    // 
970    //  - debug part of FillResolutionHistoLocal - 
971    // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data
972    // called only for GetStreamLevel() > 4
973    // fill resolution trees
974    //
975       
976    Int_t sector = cluster0->GetDetector();
977    Float_t xref = cluster0->GetX();
978    Int_t padSize = 0;                          // short pads
979    if (cluster0->GetDetector() >= 36) {
980       padSize = 1;                              // medium pads 
981       if (cluster0->GetRow() > 63) padSize = 2; // long pads
982    }
983       
984    static TLinearFitter fitY0(3, "pol2");
985    static TLinearFitter fitZ0(3, "pol2");
986    static TLinearFitter fitY2(5, "hyp4");
987    static TLinearFitter fitZ2(5, "hyp4");
988    static TLinearFitter fitY2Q(5, "hyp4");
989    static TLinearFitter fitZ2Q(5, "hyp4");
990    static TLinearFitter fitY2S(5, "hyp4");
991    static TLinearFitter fitZ2S(5, "hyp4");
992    fitY0.ClearPoints();
993    fitZ0.ClearPoints();
994    fitY2.ClearPoints();
995    fitZ2.ClearPoints();
996    fitY2Q.ClearPoints();
997    fitZ2Q.ClearPoints();
998    fitY2S.ClearPoints();
999    fitZ2S.ClearPoints();
1000    
1001    for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
1002      // loop over irow +- kDelta rows (neighboured rows)
1003      // 
1004      // 
1005      if (idelta == 0) continue;
1006      if (idelta + irow < 0 || idelta + irow > 159) continue;   // don't go out of ROC
1007      AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta);
1008      if (!cluster) continue;
1009      if (cluster->GetType() < 0) continue;
1010      if (cluster->GetDetector() != sector) continue;
1011      Double_t x = cluster->GetX() - xref;
1012      Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) );
1013      Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) );
1014      //
1015      Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1016      Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1017      Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), 
1018                                                            TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1019      Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), 
1020                                                            TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1021      Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1022                                                          TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1023                                                          TMath::Sqrt(cluster0->GetSigmaY2()), 0 );
1024      Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1025                                                         TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1026                                                         TMath::Sqrt(cluster0->GetSigmaZ2()),0 );
1027      sigmaYS  = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.);
1028      sigmaZS  = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.);
1029      //
1030      if (kDelta != 0){
1031        fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
1032        fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
1033      }
1034      Double_t xxx[4];
1035      xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0;
1036      xxx[1] = x;
1037      xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0;
1038      xxx[3] = x * x;    
1039      fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
1040      fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
1041      fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
1042      fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
1043      fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
1044      fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
1045      //
1046    }  // neigbouhood-loop
1047       //
1048    fitY0.Eval();
1049    fitZ0.Eval();
1050    fitY2.Eval();
1051    fitZ2.Eval();
1052    fitY2Q.Eval();
1053    fitZ2Q.Eval();
1054    fitY2S.Eval();
1055    fitZ2S.Eval();
1056    Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.);
1057    Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.);
1058    Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.);
1059    Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.);
1060    Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.);
1061    Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.);
1062    Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.);
1063    Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.);
1064    //
1065    static  TVectorD    parY0(3);
1066    static  TMatrixD    matY0(3, 3);
1067    static  TVectorD    parZ0(3);
1068    static  TMatrixD    matZ0(3, 3);
1069    fitY0.GetParameters(parY0);
1070    fitY0.GetCovarianceMatrix(matY0);
1071    fitZ0.GetParameters(parZ0);
1072    fitZ0.GetCovarianceMatrix(matZ0);
1073    //
1074    static  TVectorD    parY2(5);
1075    static  TMatrixD    matY2(5,5);
1076    static  TVectorD    parZ2(5);
1077    static  TMatrixD    matZ2(5,5);
1078    fitY2.GetParameters(parY2);
1079    fitY2.GetCovarianceMatrix(matY2);
1080    fitZ2.GetParameters(parZ2);
1081    fitZ2.GetCovarianceMatrix(matZ2);
1082    //
1083    static  TVectorD    parY2Q(5);
1084    static  TMatrixD    matY2Q(5,5);
1085    static  TVectorD    parZ2Q(5);
1086    static  TMatrixD    matZ2Q(5,5);
1087    fitY2Q.GetParameters(parY2Q);
1088    fitY2Q.GetCovarianceMatrix(matY2Q);
1089    fitZ2Q.GetParameters(parZ2Q);
1090    fitZ2Q.GetCovarianceMatrix(matZ2Q);
1091    static  TVectorD    parY2S(5);
1092    static  TMatrixD    matY2S(5,5);
1093    static  TVectorD    parZ2S(5);
1094    static  TMatrixD    matZ2S(5,5);
1095    fitY2S.GetParameters(parY2S);
1096    fitY2S.GetCovarianceMatrix(matY2S);
1097    fitZ2S.GetParameters(parZ2S);
1098    fitZ2S.GetCovarianceMatrix(matZ2S);
1099    Float_t sigmaY0   = TMath::Sqrt(matY0(0,0));
1100    Float_t sigmaZ0   = TMath::Sqrt(matZ0(0,0));
1101    Float_t sigmaDY0  = TMath::Sqrt(matY0(1,1));
1102    Float_t sigmaDZ0  = TMath::Sqrt(matZ0(1,1));
1103    Float_t sigmaY2   = TMath::Sqrt(matY2(1,1));
1104    Float_t sigmaZ2   = TMath::Sqrt(matZ2(1,1));
1105    Float_t sigmaDY2  = TMath::Sqrt(matY2(3,3));
1106    Float_t sigmaDZ2  = TMath::Sqrt(matZ2(3,3));
1107    Float_t sigmaY2Q  = TMath::Sqrt(matY2Q(1,1));
1108    Float_t sigmaZ2Q  = TMath::Sqrt(matZ2Q(1,1));
1109    Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
1110    Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
1111    Float_t sigmaY2S  = TMath::Sqrt(matY2S(1,1));
1112    Float_t sigmaZ2S  = TMath::Sqrt(matZ2S(1,1));
1113    Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
1114    Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
1115    
1116    // Error parameters
1117    Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
1118    Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
1119    //
1120    Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1121                                                   TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1122    Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1123                                                   TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1124    Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1125                                                         TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1126    Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1127                                                         TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1128    
1129    
1130    // RMS parameters
1131    Float_t meanRMSY = 0;
1132    Float_t meanRMSZ = 0;
1133    Int_t   nclRMS = 0;
1134    for (Int_t idelta = -2; idelta <= 2; idelta++){
1135      // loop over neighbourhood
1136      if (idelta+irow < 0 || idelta+irow > 159) continue;
1137      //         if (idelta+irow>159) continue;
1138      AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
1139      if (!cluster) continue;
1140      meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
1141      meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
1142      nclRMS++;
1143    }
1144    meanRMSY /= nclRMS; 
1145    meanRMSZ /= nclRMS; 
1146    
1147    Float_t rmsY      = TMath::Sqrt(cluster0->GetSigmaY2());  
1148    Float_t rmsZ      = TMath::Sqrt(cluster0->GetSigmaZ2());
1149    Float_t rmsYT     = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1150                                               TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1151    Float_t rmsZT     = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1152                                               TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1153    Float_t rmsYT0    = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1154                                               TMath::Abs(angley));
1155    Float_t rmsZT0    = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1156                                                  TMath::Abs(anglez));
1157    Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1158                                                   TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1159    Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1160                                                   TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1161    Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1162                                                       TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1163                                                       rmsY,meanRMSY);
1164    Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1165                                                       TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1166                                                       rmsZ,meanRMSZ);
1167    //
1168    // cluster debug
1169    TTreeSRedirector *cstream = GetDebugStreamer();
1170    if (cstream){
1171      (*cstream)<<"ResolCl"<<    // valgrind 3   40,000 bytes in 1 blocks are possibly 
1172        "Sector="<<sector<<
1173        "Cl.="<<cluster0<<
1174        "CSigmaY0="<<csigmaY0<<   // cluster errorY
1175        "CSigmaYQ="<<csigmaYQ<<   // cluster errorY - q scaled
1176        "CSigmaYS="<<csigmaYS<<   // cluster errorY - q scaled
1177        "CSigmaZ0="<<csigmaZ0<<   // 
1178        "CSigmaZQ="<<csigmaZQ<<
1179        "CSigmaZS="<<csigmaZS<<
1180        "shapeYF="<<rmsYFactor<<
1181        "shapeZF="<<rmsZFactor<<
1182        "rmsY="<<rmsY<<
1183        "rmsZ="<<rmsZ<<
1184        "rmsYM="<<meanRMSY<<
1185        "rmsZM="<<meanRMSZ<<
1186        "rmsYT="<<rmsYT<<
1187        "rmsZT="<<rmsZT<<
1188        "rmsYT0="<<rmsYT0<<
1189        "rmsZT0="<<rmsZT0<<
1190        "rmsYS="<<rmsYSigma<<  
1191        "rmsZS="<<rmsZSigma<<
1192        "padSize="<<padSize<<
1193        "Ncl="<<nclFound<<       
1194        "PY0.="<<&parY0<<
1195        "PZ0.="<<&parZ0<<
1196        "SigmaY0="<<sigmaY0<< 
1197        "SigmaZ0="<<sigmaZ0<< 
1198        "angley="<<angley<<
1199        "anglez="<<anglez<<
1200        "\n";
1201      
1202      //       tracklet dubug
1203      (*cstream)<<"ResolTr"<<    
1204        "padSize="<<padSize<<
1205        "IPad="<<padSize<<
1206        "Sector="<<sector<<
1207        "Ncl="<<nclFound<<       
1208        "chi2Y0="<<chi2Y0<<
1209        "chi2Z0="<<chi2Z0<<
1210        "chi2Y2="<<chi2Y2<<
1211        "chi2Z2="<<chi2Z2<<
1212        "chi2Y2Q="<<chi2Y2Q<<
1213        "chi2Z2Q="<<chi2Z2Q<<
1214        "chi2Y2S="<<chi2Y2S<<
1215        "chi2Z2S="<<chi2Z2S<<
1216        "PY0.="<<&parY0<<
1217        "PZ0.="<<&parZ0<<
1218        "PY2.="<<&parY2<<
1219        "PZ2.="<<&parZ2<<
1220        "PY2Q.="<<&parY2Q<<
1221        "PZ2Q.="<<&parZ2Q<<
1222        "PY2S.="<<&parY2S<<
1223        "PZ2S.="<<&parZ2S<<
1224        "SigmaY0="<<sigmaY0<< 
1225        "SigmaZ0="<<sigmaZ0<< 
1226        "SigmaDY0="<<sigmaDY0<< 
1227        "SigmaDZ0="<<sigmaDZ0<< 
1228        "SigmaY2="<<sigmaY2<< 
1229        "SigmaZ2="<<sigmaZ2<< 
1230        "SigmaDY2="<<sigmaDY2<< 
1231        "SigmaDZ2="<<sigmaDZ2<< 
1232        "SigmaY2Q="<<sigmaY2Q<< 
1233        "SigmaZ2Q="<<sigmaZ2Q<< 
1234        "SigmaDY2Q="<<sigmaDY2Q<< 
1235        "SigmaDZ2Q="<<sigmaDZ2Q<< 
1236        "SigmaY2S="<<sigmaY2S<< 
1237        "SigmaZ2S="<<sigmaZ2S<< 
1238        "SigmaDY2S="<<sigmaDY2S<< 
1239        "SigmaDZ2S="<<sigmaDZ2S<< 
1240        "angley="<<angley<<
1241        "anglez="<<anglez<<
1242        "\n";
1243    }  
1244 }
1245
1246
1247
1248
1249
1250 TH2D * AliTPCcalibTracks::MakeDiff(TH2D * hfit, TF2 * func){
1251    // 
1252    // creates a new histogram which contains the difference between
1253    // the histogram hfit and the function func
1254    // 
1255   TH2D * result = (TH2D*)hfit->Clone();      // valgrind 3   40,139 bytes in 11 blocks are still reachable
1256   result->SetTitle(Form("%s fit residuals",result->GetTitle()));
1257   result->SetName(Form("%s fit residuals",result->GetName()));
1258   TAxis *xaxis = hfit->GetXaxis();
1259   TAxis *yaxis = hfit->GetYaxis();
1260   Double_t x[2];
1261   for (Int_t biny = 0; biny <= yaxis->GetNbins(); biny++) {
1262     x[1]  = yaxis->GetBinCenter(biny);
1263     for (Int_t binx = 0; binx <= xaxis->GetNbins(); binx++) {
1264       x[0]  = xaxis->GetBinCenter(binx);
1265       Int_t bin = hfit->GetBin(binx, biny);
1266       Double_t val = hfit->GetBinContent(bin);
1267 //      result->SetBinContent( bin, (val - func->Eval(x[0], x[1])) / func->Eval(x[0], x[1]) );
1268       result->SetBinContent( bin, (val / func->Eval(x[0], x[1])) - 1 );
1269     }    
1270   }
1271   return result;
1272 }
1273
1274
1275 void  AliTPCcalibTracks::SetStyle() const {
1276    // 
1277    // set style, can be called by all draw functions
1278    // 
1279    gROOT->SetStyle("Plain");
1280    gStyle->SetFillColor(10);
1281    gStyle->SetPadColor(10);
1282    gStyle->SetCanvasColor(10);
1283    gStyle->SetStatColor(10);
1284    gStyle->SetPalette(1,0);
1285    gStyle->SetNumberContours(60);
1286 }
1287
1288
1289 void AliTPCcalibTracks::Draw(Option_t* opt){
1290    // 
1291    // draw-function of AliTPCcalibTracks
1292    // will draws some exemplaric pictures
1293    // 
1294
1295   if (GetDebugLevel() > 6) Info("Draw", "Drawing an exemplaric picture.");
1296    SetStyle();
1297    Double_t min = 0;
1298    Double_t max = 0;
1299    TCanvas *c1 = new TCanvas();
1300    c1->Divide(0, 3);
1301    TVirtualPad *upperThird = c1->GetPad(1);
1302    TVirtualPad *middleThird = c1->GetPad(2);
1303    TVirtualPad *lowerThird = c1->GetPad(3);
1304    upperThird->Divide(2,0);
1305    TVirtualPad *upleft  = upperThird->GetPad(1);
1306    TVirtualPad *upright = upperThird->GetPad(2);
1307    middleThird->Divide(2,0);
1308    TVirtualPad *middleLeft  = middleThird->GetPad(1);
1309    TVirtualPad *middleRight = middleThird->GetPad(2);
1310    lowerThird->Divide(2,0);
1311    TVirtualPad *downLeft  = lowerThird->GetPad(1);
1312    TVirtualPad *downRight = lowerThird->GetPad(2);
1313    
1314    
1315    upleft->cd(0);
1316    min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1317    max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1318    fDeltaY->SetAxisRange(min, max);
1319    fDeltaY->Fit("gaus","q","",min, max);        // valgrind 3  7 block possibly lost   2,400 bytes in 1 blocks are still reachable
1320    c1->Update();
1321    
1322    upright->cd(0);
1323    max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1324    min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1325    fDeltaZ->SetAxisRange(min, max);
1326    fDeltaZ->Fit("gaus","q","",min, max);
1327    c1->Update();
1328    
1329    middleLeft->cd();
1330    fHclus->Draw(opt);
1331    
1332    middleRight->cd();
1333    fRejectedTracksHisto->Draw(opt);
1334    TPaveText *pt = new TPaveText(0.6,0.6, 0.8,0.8, "NDC");
1335    TText *t1 = pt->AddText("1: kEdgeThetaCutNoise");
1336    TText *t2 = pt->AddText("2: kMinClusters");
1337    TText *t3 = pt->AddText("3: kMinRatio");
1338    TText *t4 = pt->AddText("4: kMax1pt");
1339    t1 = t1; t2 = t2; t3 = t3; t4 = t4;    // avoid compiler warnings
1340    pt->SetToolTipText("Legend for failed cuts");
1341    pt->Draw();
1342    
1343    downLeft->cd();
1344    fHclusterPerPadrowRaw->Draw(opt);
1345    
1346    downRight->cd();
1347    fHclusterPerPadrow->Draw(opt);
1348 }
1349
1350
1351 void AliTPCcalibTracks::MakeReport(Int_t stat, char* pathName){ 
1352    // 
1353    // all functions are called, that produce pictures
1354    // the histograms are written to the directory 'pathName'
1355    // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
1356    // 'stat' is also the number of minEntries for MakeResPlotsQTree
1357    // 
1358
1359   if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
1360    MakeAmpPlots(stat, pathName);
1361    MakeDeltaPlots(pathName);
1362    FitResolutionNew(pathName);
1363    FitRMSNew(pathName);
1364    MakeChargeVsDriftLengthPlots(pathName);
1365 //    MakeResPlotsQ(1, 1); 
1366    MakeResPlotsQTree(stat, pathName);
1367 }
1368    
1369
1370 void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){ 
1371    // 
1372    // creates several plots:
1373    // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps
1374    // fArrayAmp.ps: one histogram per sector, the histogram shows the charge per cluster
1375    // fArrayAmpRow.ps: one histogram per sector, mean max. amplitude vs. pad row with landau fit
1376    // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1377    // Empty histograms (sectors without data) are not written to file
1378    // the ps-files are written to the directory 'pathName', that is created if it does not exist
1379    // 'stat': only histograms with more than 'stat' entries are written to file.
1380    // 
1381    
1382    SetStyle();
1383    gSystem->MakeDirectory(pathName);
1384    gSystem->ChangeDirectory(pathName);
1385    
1386    TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1387    TPostScript *ps; 
1388    // histograms with accumulated amplitude for all IROCs and OROCs
1389    TH1F *allAmpHisIROC = ((TH1F*)(fArrayAmp->At(0))->Clone());
1390    allAmpHisIROC->SetName("Amp all IROCs");
1391    allAmpHisIROC->SetTitle("Amp all IROCs");
1392    TH1F *allAmpHisOROC = ((TH1F*)(fArrayAmp->At(36))->Clone());
1393    allAmpHisOROC->SetName("Amp all OROCs");
1394    allAmpHisOROC->SetTitle("Amp all OROCs");
1395    
1396    ps = new TPostScript("fArrayAmp.ps", 112);
1397    if (GetDebugLevel() > -1) cout << "creating fArrayAmp.ps..." << endl;
1398    for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++){
1399       if ( ((TH1F*)fArrayAmp->At(i))->GetEntries() < stat  ) continue;
1400       ps->NewPage();
1401       ((TH1F*)fArrayAmp->At(i))->Draw();
1402       c1->Update();              // valgrind 3
1403       if (i > 0 && i < 36) {
1404          allAmpHisIROC->Add(((TH1F*)fArrayAmp->At(i)));
1405          allAmpHisOROC->Add(((TH1F*)fArrayAmp->At(i+36)));
1406       }
1407    }
1408    ps->NewPage();
1409    allAmpHisIROC->Draw();
1410    c1->Update();              // valgrind
1411    ps->NewPage();
1412    allAmpHisOROC->Draw();
1413    c1->Update();
1414    ps->Close();
1415    delete ps;
1416    
1417    TH1F *his = 0;
1418    Double_t min = 0;
1419    Double_t max = 0;
1420    ps = new TPostScript("fArrayAmpRow.ps", 112);
1421    if (GetDebugLevel() > -1) cout << "creating fArrayAmpRow.ps..." << endl;
1422    for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++){
1423       his = (TH1F*)fArrayAmpRow->At(i);
1424       if (his->GetEntries() < stat) continue;
1425       ps->NewPage();
1426       min = TMath::Max( his->GetBinCenter(his->GetMaximumBin() )-100., 0.);
1427       max = his->GetBinCenter(5*his->GetMaximumBin()) + 100;
1428       his->SetAxisRange(min, max);
1429       his->Fit("pol3", "q", "", min, max);
1430       // his->Draw("error");    // don't use this line when you don't want to have empty pages in the ps-file
1431       c1->Update();
1432    }
1433    ps->Close();
1434    delete ps;
1435    delete c1;
1436    gSystem->ChangeDirectory("..");
1437 }
1438
1439
1440 void AliTPCcalibTracks::MakeDeltaPlots(char* pathName){
1441    // 
1442    // creates several plots:
1443    // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1444    // the ps-files are written to the directory 'pathName', that is created if it does not exist
1445    // 
1446    
1447    SetStyle();
1448    gSystem->MakeDirectory(pathName);
1449    gSystem->ChangeDirectory(pathName);
1450    
1451    TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1452    TPostScript *ps; 
1453    Double_t min = 0;
1454    Double_t max = 0;
1455    
1456    ps = new TPostScript("DeltaYZ.ps", 112);
1457    if (GetDebugLevel() > -1) cout << "creating DeltaYZ.ps..." << endl;
1458    min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1459    max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1460    fDeltaY->SetAxisRange(min, max);
1461    ps->NewPage();
1462    fDeltaY->Fit("gaus","q","",min, max);        // valgrind 3  7 block possibly lost   2,400 bytes in 1 blocks are still reachable
1463    c1->Update();
1464    ps->NewPage();
1465    max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1466    min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1467    fDeltaZ->SetAxisRange(min, max);
1468    fDeltaZ->Fit("gaus","q","",min, max);
1469    c1->Update();
1470    ps->Close();
1471    delete ps;
1472    delete c1;
1473    gSystem->ChangeDirectory("..");
1474 }
1475
1476
1477 void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(char* pathName){
1478    // 
1479    // creates charge vs. driftlength plots, one TProfile for each ROC
1480    // is not correct like this, should be one TProfile for each sector and padsize
1481    // 
1482    
1483    SetStyle();
1484    gSystem->MakeDirectory(pathName);
1485    gSystem->ChangeDirectory(pathName);
1486    
1487    TCanvas* c1 = new TCanvas();     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1488    TPostScript *ps; 
1489    ps = new TPostScript("chargeVsDriftlengthOld.ps", 112);
1490    if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlength.ps..." << endl;
1491    TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone());
1492    TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone());
1493    chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC");
1494    chargeVsDriftlengthAllIROCs->SetTitle("charge vs. driftlength, all IROCs");
1495    chargeVsDriftlengthAllOROCs->SetName("allAmpHisOROC");
1496    chargeVsDriftlengthAllOROCs->SetTitle("charge vs. driftlength, all OROCs");
1497    
1498    for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) {
1499       ((TProfile*)fArrayChargeVsDriftlength->At(i))->Draw();
1500       c1->Update();
1501       if (i > 0 && i < 36) { 
1502          chargeVsDriftlengthAllIROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i)));
1503          chargeVsDriftlengthAllOROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i+36)));
1504       }
1505       ps->NewPage();
1506    }
1507    chargeVsDriftlengthAllIROCs->Draw();
1508    c1->Update();              // valgrind
1509    ps->NewPage();
1510    chargeVsDriftlengthAllOROCs->Draw();
1511    c1->Update();
1512    ps->Close();  
1513    delete ps;
1514    delete c1;
1515    gSystem->ChangeDirectory("..");
1516 }   
1517
1518
1519 void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(char* pathName){
1520    // 
1521    // creates charge vs. driftlength plots, one TProfile for each ROC
1522    // under development....
1523    // 
1524    
1525    SetStyle();
1526    gSystem->MakeDirectory(pathName);
1527    gSystem->ChangeDirectory(pathName);
1528    
1529    TCanvas* c1 = new TCanvas("c1", "c1", 700,(Int_t)(TMath::Sqrt(2)*700));     // valgrind 3 ???  634 bytes in 28 blocks are still reachable
1530 //    TCanvas c1("c1", "c1", 500,(sqrt(2)*500))
1531    c1->Divide(0,3);
1532    TPostScript *ps; 
1533    ps = new TPostScript("chargeVsDriftlength.ps", 111);
1534    if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl;
1535    
1536    TProfile *chargeVsDriftlengthAllShortPads  = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone());
1537    TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone());
1538    TProfile *chargeVsDriftlengthAllLongPads   = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)->Clone());
1539    chargeVsDriftlengthAllShortPads->SetName("allAmpHisShortPads");
1540    chargeVsDriftlengthAllShortPads->SetTitle("charge vs. driftlength, all sectors, short pads");
1541    chargeVsDriftlengthAllMediumPads->SetName("allAmpHisMediumPads");
1542    chargeVsDriftlengthAllMediumPads->SetTitle("charge vs. driftlength, all sectors, medium pads");
1543    chargeVsDriftlengthAllLongPads->SetName("allAmpHisLongPads");
1544    chargeVsDriftlengthAllLongPads->SetTitle("charge vs. driftlength, all sectors, long pads");
1545    
1546    for (Int_t i = 0; i < 36; i++) {
1547       c1->cd(1)->SetGridx();
1548       c1->cd(1)->SetGridy();
1549       ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,0))->Draw();
1550       c1->cd(2)->SetGridx();
1551       c1->cd(2)->SetGridy();
1552       ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,1))->Draw();
1553       c1->cd(3)->SetGridx();
1554       c1->cd(3)->SetGridy();
1555       ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,2))->Draw();
1556       c1->Update();
1557       chargeVsDriftlengthAllShortPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0));
1558       chargeVsDriftlengthAllMediumPads->Add((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1));
1559       chargeVsDriftlengthAllLongPads->Add(  (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2));
1560       ps->NewPage();
1561    }
1562    c1->cd(1)->SetGridx();
1563    c1->cd(1)->SetGridy();
1564    chargeVsDriftlengthAllShortPads->Draw();
1565    c1->cd(2)->SetGridx();
1566    c1->cd(2)->SetGridy();
1567    chargeVsDriftlengthAllMediumPads->Draw();
1568    c1->cd(3)->SetGridx();
1569    c1->cd(3)->SetGridy();
1570    chargeVsDriftlengthAllLongPads->Draw();
1571    c1->Update();              // valgrind
1572 //    ps->NewPage();
1573    ps->Close();  
1574    delete ps;
1575    delete c1;
1576    gSystem->ChangeDirectory("..");
1577 }   
1578
1579
1580
1581 void AliTPCcalibTracks::FitResolutionNew(char* pathName){
1582    // 
1583    // calculates different resulution fits in Y and Z direction
1584    // the histograms are written to 'ResolutionYZ.ps'
1585    // writes calculated resolution to 'resol.txt'
1586    // all files are stored in the directory pathName
1587    // 
1588    
1589    SetStyle();
1590    gSystem->MakeDirectory(pathName);
1591    gSystem->ChangeDirectory(pathName);
1592    
1593    TCanvas c;
1594    c.Divide(2,1); 
1595    if (GetDebugLevel() > -1) cout << "creating ResolutionYZ.ps..." << endl;
1596    TPostScript *ps = new TPostScript("ResolutionYZ.ps", 112); 
1597    TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1598    fres->SetParameter(0,0.02);
1599    fres->SetParameter(1,0.0054);
1600    fres->SetParameter(2,0.13);  
1601    
1602    TH1::AddDirectory(kTRUE);  // TH3F::FitSlicesZ() writes histograms into the current directory
1603    
1604    // create histogramw for Y-resolution
1605    TH3F * hisResY0 = (TH3F*)fResolY->At(0);
1606    hisResY0->FitSlicesZ();
1607    TH2D * hisResY02 = (TH2D*)gDirectory->Get("Resol Y0_2");
1608    TH3F * hisResY1 = (TH3F*)fResolY->At(1); 
1609    hisResY1->FitSlicesZ();
1610    TH2D * hisResY12 = (TH2D*)gDirectory->Get("Resol Y1_2");
1611    TH3F * hisResY2 = (TH3F*)fResolY->At(2);
1612    hisResY2->FitSlicesZ();
1613    TH2D * hisResY22 = (TH2D*)gDirectory->Get("Resol Y2_2");
1614     //
1615    ps->NewPage();
1616    c.cd(1);
1617    hisResY02->Fit(fres, "q");      // valgrind    132,072 bytes in 6 blocks are indirectly lost
1618    hisResY02->Draw("surf1"); 
1619    c.cd(2);
1620    MakeDiff(hisResY02,fres)->Draw("surf1");
1621    c.Update();
1622    //   c.SaveAs("ResolutionYPad0.eps");
1623    ps->NewPage();
1624    c.cd(1);
1625    hisResY12->Fit(fres, "q");
1626    hisResY12->Draw("surf1");
1627    c.cd(2);
1628    MakeDiff(hisResY12,fres)->Draw("surf1");
1629    c.Update();
1630    //   c.SaveAs("ResolutionYPad1.eps");
1631    ps->NewPage();
1632    c.cd(1);
1633    hisResY22->Fit(fres, "q");
1634    hisResY22->Draw("surf1"); 
1635    c.cd(2);
1636    MakeDiff(hisResY22,fres)->Draw("surf1");
1637    c.Update();
1638 //    c.SaveAs("ResolutionYPad2.eps");
1639    
1640    // create histogramw for Z-resolution
1641    TH3F * hisResZ0 = (TH3F*)fResolZ->At(0);
1642    hisResZ0->FitSlicesZ();
1643    TH2D * hisResZ02 = (TH2D*)gDirectory->Get("Resol Z0_2");
1644    TH3F * hisResZ1 = (TH3F*)fResolZ->At(1);
1645    hisResZ1->FitSlicesZ();
1646    TH2D * hisResZ12 = (TH2D*)gDirectory->Get("Resol Z1_2");
1647    TH3F * hisResZ2 = (TH3F*)fResolZ->At(2);
1648    hisResZ2->FitSlicesZ();
1649    TH2D * hisResZ22 = (TH2D*)gDirectory->Get("Resol Z2_2");
1650    
1651    ps->NewPage();
1652    c.cd(1);
1653    hisResZ02->Fit(fres, "q");
1654    hisResZ02->Draw("surf1");
1655    c.cd(2);
1656    MakeDiff(hisResZ02,fres)->Draw("surf1");
1657    c.Update();
1658 //    c.SaveAs("ResolutionZPad0.eps");
1659    ps->NewPage();
1660    c.cd(1);
1661    hisResZ12->Fit(fres, "q");
1662    hisResZ12->Draw("surf1"); 
1663    c.cd(2);
1664    MakeDiff(hisResZ12,fres)->Draw("surf1");
1665    c.Update();
1666 //    c.SaveAs("ResolutionZPad1.eps");
1667    ps->NewPage();
1668    c.cd(1);
1669    hisResZ22->Fit(fres, "q");
1670    hisResZ22->Draw("surf1");  
1671    c.cd(2);
1672    MakeDiff(hisResZ22,fres)->Draw("surf1");
1673    c.Update();
1674 //    c.SaveAs("ResolutionZPad2.eps");
1675    ps->Close();
1676    delete ps;
1677    
1678    // write calculated resoltuions to 'resol.txt'
1679    ofstream fresol("resol.txt");
1680    fresol<<"Pad 0.75 cm"<<"\n";
1681    hisResY02->Fit(fres, "q");                     // valgrind
1682    fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1683    hisResZ02->Fit(fres, "q");
1684    fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1685    //
1686    fresol<<"Pad 1.00 cm"<<1<<"\n";
1687    hisResY12->Fit(fres, "q");                     // valgrind
1688    fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1689    hisResZ12->Fit(fres, "q");
1690    fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1691    //
1692    fresol<<"Pad 1.50 cm"<<0<<"\n";
1693    hisResY22->Fit(fres, "q");
1694    fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1695    hisResZ22->Fit(fres, "q");
1696    fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1697    
1698    TH1::AddDirectory(kFALSE);
1699    gSystem->ChangeDirectory("..");
1700    delete fres;
1701 }
1702
1703
1704 void AliTPCcalibTracks::FitRMSNew(char* pathName){
1705    // 
1706    // calculates different resulution-rms fits in Y and Z direction
1707    // the histograms are written to 'RMS_YZ.ps'
1708    // writes calculated resolution rms to 'rms.txt'
1709    // all files are stored in the directory pathName
1710    // 
1711    
1712    SetStyle();
1713    gSystem->MakeDirectory(pathName);
1714    gSystem->ChangeDirectory(pathName);
1715    
1716    TCanvas c;        // valgrind 3   42,120 bytes in 405 blocks are still reachable   23,816 bytes in 229 blocks are still reachable
1717    c.Divide(2,1); 
1718    if (GetDebugLevel() > -1) cout << "creating RMS_YZ.ps..." << endl;
1719    TPostScript *ps = new TPostScript("RMS_YZ.ps", 112); 
1720    TF2 *frms = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1721    frms->SetParameter(0,0.02);
1722    frms->SetParameter(1,0.0054);
1723    frms->SetParameter(2,0.13);  
1724    
1725    TH1::AddDirectory(kTRUE);  // TH3F::FitSlicesZ() writes histograms into the current directory
1726    
1727    // create histogramw for Y-RMS   
1728    TH3F * hisResY0 = (TH3F*)fRMSY->At(0);
1729    hisResY0->FitSlicesZ();
1730    TH2D * hisResY02 = (TH2D*)gDirectory->Get("RMS Y0_1");
1731    TH3F * hisResY1 = (TH3F*)fRMSY->At(1);
1732    hisResY1->FitSlicesZ();
1733    TH2D * hisResY12 = (TH2D*)gDirectory->Get("RMS Y1_1");
1734    TH3F * hisResY2 = (TH3F*)fRMSY->At(2);
1735    hisResY2->FitSlicesZ();
1736    TH2D * hisResY22 = (TH2D*)gDirectory->Get("RMS Y2_1");
1737    //
1738    ps->NewPage();
1739    c.cd(1);
1740    hisResY02->Fit(frms, "qn0"); 
1741    hisResY02->Draw("surf1"); 
1742    c.cd(2);
1743    MakeDiff(hisResY02,frms)->Draw("surf1");
1744    c.Update();
1745 //    c.SaveAs("RMSYPad0.eps");
1746    ps->NewPage();
1747    c.cd(1);
1748    hisResY12->Fit(frms, "qn0");               // valgrind   several blocks possibly lost
1749    hisResY12->Draw("surf1");
1750    c.cd(2);
1751    MakeDiff(hisResY12,frms)->Draw("surf1");
1752    c.Update();
1753 //    c.SaveAs("RMSYPad1.eps");
1754    ps->NewPage();
1755    c.cd(1);
1756    hisResY22->Fit(frms, "qn0");
1757    hisResY22->Draw("surf1"); 
1758    c.cd(2);
1759    MakeDiff(hisResY22,frms)->Draw("surf1");
1760    c.Update();
1761 //    c.SaveAs("RMSYPad2.eps");
1762    
1763    // create histogramw for Z-RMS   
1764    TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0);
1765    hisResZ0->FitSlicesZ();
1766    TH2D * hisResZ02 = (TH2D*)gDirectory->Get("RMS Z0_1");
1767    TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1); 
1768    hisResZ1->FitSlicesZ();
1769    TH2D * hisResZ12 = (TH2D*)gDirectory->Get("RMS Z1_1");
1770    TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2); 
1771    hisResZ2->FitSlicesZ();
1772    TH2D * hisResZ22 = (TH2D*)gDirectory->Get("RMS Z2_1");
1773    //
1774    ps->NewPage();
1775    c.cd(1);
1776    hisResZ02->Fit(frms, "qn0");         // valgrind
1777    hisResZ02->Draw("surf1");
1778    c.cd(2);
1779    MakeDiff(hisResZ02,frms)->Draw("surf1");
1780    c.Update();
1781 //    c.SaveAs("RMSZPad0.eps");
1782    ps->NewPage();
1783    c.cd(1);
1784    hisResZ12->Fit(frms, "qn0");
1785    hisResZ12->Draw("surf1"); 
1786    c.cd(2);
1787    MakeDiff(hisResZ12,frms)->Draw("surf1");
1788    c.Update();
1789 //    c.SaveAs("RMSZPad1.eps");
1790    ps->NewPage();
1791    c.cd(1);
1792    hisResZ22->Fit(frms, "qn0");         // valgrind  1 block possibly lost
1793    hisResZ22->Draw("surf1");  
1794    c.cd(2);
1795    MakeDiff(hisResZ22,frms)->Draw("surf1");
1796    c.Update();
1797 //    c.SaveAs("RMSZPad2.eps");
1798    
1799    // write calculated resoltuion rms to 'rms.txt'
1800    ofstream filerms("rms.txt");
1801    filerms<<"Pad 0.75 cm"<<"\n";
1802    hisResY02->Fit(frms, "qn0");         // valgrind   23 blocks indirectly lost
1803    filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1804    hisResZ02->Fit(frms, "qn0");         // valgrind   23 blocks indirectly lost
1805    filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1806    //
1807    filerms<<"Pad 1.00 cm"<<1<<"\n";
1808    hisResY12->Fit(frms, "qn0");         // valgrind      3,256 bytes in 22 blocks are indirectly lost 
1809    filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1810    hisResZ12->Fit(frms, "qn0");         // valgrind    66,036 bytes in 3 blocks are still reachable
1811    filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1812    //
1813    filerms<<"Pad 1.50 cm"<<0<<"\n";
1814    hisResY22->Fit(frms, "qn0");      // valgrind   40,139 bytes in 11 blocks are still reachable   330,180 bytes in 15 blocks are possibly lost
1815    filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1816    hisResZ22->Fit(frms, "qn0");
1817    filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1818
1819    TH1::AddDirectory(kFALSE);
1820    gSystem->ChangeDirectory("..");
1821    ps->Close();
1822    delete ps;
1823    delete frms;
1824 }
1825
1826
1827 void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){
1828    //
1829    // Make tree with resolution parameters
1830    // the result is written to 'resol.root' in directory 'pathname'
1831    // file information are available in fileInfo
1832    // available variables in the tree 'Resol':
1833    //  Entries: number of entries for this resolution point
1834    //  nbins:   number of bins that were accumulated
1835    //  Dim:     direction, Dim==0: y-direction, Dim==1: z-direction
1836    //  Pad:     padSize; short, medium and long
1837    //  Length:  pad length, 0.75, 1, 1.5
1838    //  QMean:   mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1839    //  Zc:      center of middle bin in drift direction
1840    //  Zm:      mean dirftlength for accumulated Delta-Histograms
1841    //  Zs:      width of driftlength bin
1842    //  AngleC:  center of middle bin in Angle-Direction
1843    //  AngleM:  mean angle for accumulated Delta-Histograms
1844    //  AngleS:  width of Angle-bin
1845    //  Resol:   sigma for gaus fit through Delta-Histograms
1846    //  Sigma:   error of sigma for gaus fit through Delta Histograms
1847    //  MeanR:   mean of the Delta-Histogram
1848    //  SigmaR:  rms of the Delta-Histogram
1849    //  RMSm:    mean of the gaus fit through RMS-Histogram
1850    //  RMS:     sigma of the gaus fit through RMS-Histogram
1851    //  RMSe0:   error of mean of gaus fit in RMS-Histogram
1852    //  RMSe1:   error of sigma of gaus fit in RMS-Histogram
1853    //  
1854       
1855   if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
1856   if (GetDebugLevel() > -1) cout << "    relax, the calculation will take a while..." << endl;
1857   
1858    gSystem->MakeDirectory(pathName);
1859    gSystem->ChangeDirectory(pathName);
1860    TString kFileName = "resol.root";
1861    TTreeSRedirector fTreeResol(kFileName.Data());
1862    
1863    TH3F *resArray[2][3][11];
1864    TH3F *rmsArray[2][3][11];
1865   
1866    // load histograms from fArrayQDY and fArrayQDZ 
1867    // into resArray and rmsArray
1868    // that is all we need here
1869    for (Int_t idim = 0; idim < 2; idim++){
1870       for (Int_t ipad = 0; ipad < 3; ipad++){
1871          for (Int_t iq = 0; iq <= 10; iq++){
1872             rmsArray[idim][ipad][iq]=0;
1873             resArray[idim][ipad][iq]=0;
1874             Int_t bin = GetBin(iq,ipad); 
1875             TH3F *hresl = 0;
1876             if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
1877             if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
1878             if (!hresl) continue;
1879             resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
1880             resArray[idim][ipad][iq]->SetDirectory(0);
1881             TH3F * hreslRMS = 0;
1882             if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
1883             if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
1884             if (!hreslRMS) continue;
1885             rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
1886             rmsArray[idim][ipad][iq]->SetDirectory(0);
1887          }
1888       }
1889    }
1890    if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
1891    
1892    //--------------------------------------------------------------------------------------------
1893    
1894    char name[200];
1895    Double_t qMean;
1896    Double_t zMean, angleMean, zCenter, angleCenter;
1897    Double_t zSigma, angleSigma;
1898    TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
1899    TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
1900    TF1 *fitFunction = new TF1("fitFunction", "gaus");
1901    Float_t entriesQ = 0;
1902    Int_t loopCounter = 1;
1903   
1904    for (Int_t idim = 0; idim < 2; idim++){
1905       // Loop y-z corrdinate
1906       for (Int_t ipad = 0; ipad < 3; ipad++){
1907          // loop pad type
1908          for (Int_t iq = -1; iq < 10; iq++){
1909             // LOOP Q
1910            if (GetDebugLevel() > -1) 
1911                cout << "Loop-counter, this is loop " << loopCounter << " of 66, (" 
1912                   << (Int_t)((loopCounter)/66.*100) << "% done), " 
1913                   << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << "  \r" << std::flush;
1914             loopCounter++;
1915             TH3F *hres = 0;
1916             TH3F *hrms = 0;
1917             qMean    = 0;
1918             entriesQ = 0;
1919             
1920             // calculate qMean
1921             if (iq == -1){
1922                // integrated spectra
1923                for (Int_t iql = 0; iql < 10; iql++){    
1924                   Int_t bin = GetBin(iql,ipad); 
1925                   TH3F *hresl = resArray[idim][ipad][iql];
1926                   TH3F *hrmsl = rmsArray[idim][ipad][iql];
1927                   if (!hresl) continue;
1928                   if (!hrmsl) continue;     
1929                   entriesQ += hresl->GetEntries();
1930                   qMean += hresl->GetEntries() * GetQ(bin);      
1931                   if (!hres) {
1932                      hres = (TH3F*)hresl->Clone();
1933                      hrms = (TH3F*)hrmsl->Clone();
1934                   }
1935                   else{
1936                      hres->Add(hresl);
1937                      hrms->Add(hrmsl);
1938                   }
1939                }
1940                qMean /= entriesQ;
1941                qMean *= -1.;  // integral mean charge
1942             }
1943             else {
1944                // loop over neighboured Q-bins 
1945                // accumulate entries from neighboured Q-bins
1946                for (Int_t iql = iq - 1; iql <= iq + 1; iql++){              
1947                   if (iql < 0) continue;
1948                   Int_t bin = GetBin(iql,ipad);
1949                   TH3F * hresl = resArray[idim][ipad][iql];
1950                   TH3F * hrmsl = rmsArray[idim][ipad][iql];
1951                   if (!hresl) continue;
1952                   if (!hrmsl) continue;
1953                   entriesQ += hresl->GetEntries(); 
1954                   qMean += hresl->GetEntries() * GetQ(bin);      
1955                   if (!hres) {
1956                      hres = (TH3F*) hresl->Clone();
1957                      hrms = (TH3F*) hrmsl->Clone();
1958                   }
1959                   else{
1960                      hres->Add(hresl);
1961                      hrms->Add(hrmsl);
1962                   }
1963                }
1964                qMean/=entriesQ;
1965             }
1966       
1967             TAxis *xAxisDriftLength = hres->GetXaxis();   // driftlength / z - axis
1968             TAxis *yAxisAngle       = hres->GetYaxis();   // angle axis
1969             TAxis *zAxisDelta       = hres->GetZaxis();   // delta axis
1970             TAxis *zAxisRms         = hrms->GetZaxis();   // rms axis
1971             
1972             // loop over all angle bins
1973             for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
1974                angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
1975                // loop over all driftlength bins
1976                for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
1977                   zCenter    = xAxisDriftLength->GetBinCenter(ibinxDL);
1978                   zSigma     = xAxisDriftLength->GetBinWidth(ibinxDL);
1979                   angleSigma = yAxisAngle->GetBinWidth(ibinyAngle); 
1980                   zMean      = zCenter;      // changens, when more statistic is accumulated
1981                   angleMean  = angleCenter;  // changens, when more statistic is accumulated
1982                   
1983                   // create 2 1D-Histograms, projectionRes and projectionRms
1984                   // these histograms are delta histograms for given direction, padSize, chargeBin,
1985                   // angleBin and driftLengthBin
1986                   // later on they will be fitted with a gausian, its sigma is the resoltuion...
1987                   sprintf(name,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
1988                   // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1989                   projectionRes->SetNameTitle(name, name);
1990                   sprintf(name,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
1991                   // TH1D * projectionRms =  new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1992                   projectionRms->SetNameTitle(name, name);
1993                   
1994                   projectionRes->Reset();
1995                   projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1996                   projectionRms->Reset();
1997                   projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1998                   projectionRes->SetDirectory(0);
1999                   projectionRms->SetDirectory(0);
2000                   
2001                   Double_t entries = 0;
2002                   Int_t    nbins   = 0;   // counts, how many bins were accumulated
2003                   zMean     = 0;
2004                   angleMean = 0;
2005                   entries   = 0;
2006                   
2007                   // fill projectionRes and projectionRms for given dim, ipad and iq, 
2008                   // as well as for given angleBin and driftlengthBin
2009                   // if this gives not enough statistic, include neighbourhood 
2010                   // (angle and driftlength) successifely
2011                   for (Int_t dbin = 0; dbin <= 8; dbin++){              // delta-bins around centered angleBin and driftlengthBin
2012                      for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) {   // delta-bins in angle direction
2013                         for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
2014                            if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue;   // add each bin only one time !
2015                            Int_t binx2 = ibinxDL + dbinx2;                       // position variable in x (driftlength) direction
2016                            Int_t biny2 = ibinyAngle + dbiny2;                    // position variable in y (angle)  direction
2017                            if (binx2 < 1 || biny2 < 1) continue;                 // don't go out of the histogram!
2018                            if (binx2 >= xAxisDriftLength->GetNbins()) continue;  // don't go out of the histogram!
2019                            if (biny2 >= yAxisAngle->GetNbins()) continue;        // don't go out of the histogram!
2020                            nbins++;                                              // count the number of accumulated bins
2021                            // Fill resolution histo
2022                            for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
2023                               // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3);     // unused variable
2024                               projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
2025                               entries   += hres->GetBinContent(binx2, biny2, ibin3);
2026                               zMean     += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
2027                               angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
2028                            }  // ibin3 loop
2029                            // fill RMS histo
2030                            for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
2031                               projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
2032                            }
2033                         }  //dbinx2 loop
2034                         if (entries > minEntries) break; // enough statistic accumulated
2035                      }  // dbiny2 loop
2036                      if (entries > minEntries) break;    // enough statistic accumulated
2037                   }  // dbin loop
2038                   if ( entries< minEntries) continue;  // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree  
2039                   zMean /= entries;
2040                   angleMean /= entries;
2041                   
2042                   if (entries > minEntries) {
2043                      //  when enough statistic is accumulated
2044                      //  fit Delta histograms with a gausian
2045                      //  of the gausian is the resolution (resol), its fit error is sigma
2046                      //  store also mean and RMS of the histogram
2047                      Float_t xmin     = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2048                      Float_t xmax     = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2049                      
2050 //                      projectionRes->Fit("gaus", "q0", "", xmin, xmax);
2051 //                      Float_t resol    = projectionRes->GetFunction("gaus")->GetParameter(2);
2052 //                      Float_t sigma    = projectionRes->GetFunction("gaus")->GetParError(2);
2053                      fitFunction->SetMaximum(xmax);
2054                      fitFunction->SetMinimum(xmin);
2055                      projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
2056                      Float_t resol    = fitFunction->GetParameter(2);
2057                      Float_t sigma    = fitFunction->GetParError(2);
2058                      
2059                      Float_t meanR    = projectionRes->GetMean();
2060                      Float_t sigmaR   = projectionRes->GetRMS();
2061                      // fit also RMS histograms with a gausian
2062                      // store mean and sigma of the gausian in rmsMean and rmsSigma
2063                      // store also the fit errors in errorRMS and errorSigma
2064                      xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2065                      xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2066                      
2067 //                      projectionRms->Fit("gaus","q0","",xmin,xmax);
2068 //                      Float_t rmsMean    = projectionRms->GetFunction("gaus")->GetParameter(1);
2069 //                      Float_t rmsSigma   = projectionRms->GetFunction("gaus")->GetParameter(2);
2070 //                      Float_t errorRMS   = projectionRms->GetFunction("gaus")->GetParError(1);
2071 //                      Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
2072                      projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
2073                      Float_t rmsMean    = fitFunction->GetParameter(1);
2074                      Float_t rmsSigma   = fitFunction->GetParameter(2);
2075                      Float_t errorRMS   = fitFunction->GetParError(1);
2076                      Float_t errorSigma = fitFunction->GetParError(2);
2077                     
2078                      Float_t length = 0.75;
2079                      if (ipad == 1) length = 1;
2080                      if (ipad == 2) length = 1.5;
2081                      
2082                      fTreeResol<<"Resol"<<
2083                         "Entries="<<entries<<      // number of entries for this resolution point
2084                         "nbins="<<nbins<<          // number of bins that were accumulated
2085                         "Dim="<<idim<<             // direction, Dim==0: y-direction, Dim==1: z-direction
2086                         "Pad="<<ipad<<             // padSize; short, medium and long
2087                         "Length="<<length<<        // pad length, 0.75, 1, 1.5
2088                         "QMean="<<qMean<<          // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
2089                         "Zc="<<zCenter<<           // center of middle bin in drift direction
2090                         "Zm="<<zMean<<             // mean dirftlength for accumulated Delta-Histograms
2091                         "Zs="<<zSigma<<            // width of driftlength bin
2092                         "AngleC="<<angleCenter<<   // center of middle bin in Angle-Direction
2093                         "AngleM="<<angleMean<<     // mean angle for accumulated Delta-Histograms
2094                         "AngleS="<<angleSigma<<    // width of Angle-bin
2095                         "Resol="<<resol<<          // sigma for gaus fit through Delta-Histograms
2096                         "Sigma="<<sigma<<          // error of sigma for gaus fit through Delta Histograms
2097                         "MeanR="<<meanR<<          // mean of the Delta-Histogram
2098                         "SigmaR="<<sigmaR<<        // rms of the Delta-Histogram
2099                         "RMSm="<<rmsMean<<         // mean of the gaus fit through RMS-Histogram
2100                         "RMSs="<<rmsSigma<<        // sigma of the gaus fit through RMS-Histogram
2101                         "RMSe0="<<errorRMS<<       // error of mean of gaus fit in RMS-Histogram
2102                         "RMSe1="<<errorSigma<<     // error of sigma of gaus fit in RMS-Histogram
2103                         "\n";
2104                      if (GetDebugLevel() > 5) {
2105                         projectionRes->SetDirectory(fTreeResol.GetFile());
2106                         projectionRes->Write(projectionRes->GetName());
2107                         projectionRes->SetDirectory(0);
2108                         projectionRms->SetDirectory(fTreeResol.GetFile());
2109                         projectionRms->Write(projectionRms->GetName());
2110                         projectionRes->SetDirectory(0);
2111                      }
2112                   }  // if (projectionRes->GetSum() > minEntries)
2113                }  // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
2114             }  // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
2115             
2116          }  // iq-loop
2117       }  // ipad-loop
2118    }  // idim-loop
2119    delete projectionRes;
2120    delete projectionRms;
2121    
2122 //    TFile resolFile(fTreeResol.GetFile());
2123    TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
2124    fileInfo.Write("fileInfo");
2125 //    resolFile.Close();
2126 //    fTreeResol.GetFile()->Close();
2127    if (GetDebugLevel() > -1) cout << endl;
2128    if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
2129    gSystem->ChangeDirectory("..");
2130 }
2131
2132
2133
2134
2135
2136 Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
2137    // 
2138    // function to merge several AliTPCcalibTracks objects after PROOF calculation
2139    // The object's histograms are merged via their merge functions
2140    // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
2141    // 
2142    
2143   if (GetDebugLevel() > 0) cout << " *****  this is AliTPCcalibTracks::Merge(TCollection *collectionList)  *****"<< endl;  
2144    if (!collectionList) return 0;
2145    if (collectionList->IsEmpty()) return -1;
2146    
2147    if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl;     //    REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
2148    if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
2149    collectionList->Print();
2150    
2151    // create a list for each data member
2152    TList* deltaYList = new TList;
2153    TList* deltaZList = new TList;
2154    TList* arrayAmpRowList = new TList;
2155    TList* rejectedTracksList = new TList;
2156    TList* hclusList = new TList;
2157    TList* clusterCutHistoList = new TList;
2158    TList* arrayAmpList = new TList;
2159    TList* arrayQDYList = new TList;
2160    TList* arrayQDZList = new TList;
2161    TList* arrayQRMSYList = new TList;
2162    TList* arrayQRMSZList = new TList;
2163    TList* arrayChargeVsDriftlengthList = new TList;
2164    TList* calPadRegionChargeVsDriftlengthList = new TList;
2165    TList* hclusterPerPadrowList = new TList;
2166    TList* hclusterPerPadrowRawList = new TList;
2167    TList* resolYList = new TList;
2168    TList* resolZList = new TList;
2169    TList* rMSYList = new TList;
2170    TList* rMSZList = new TList;
2171    
2172 //    TList* nRowsList = new TList;
2173 //    TList* nSectList = new TList;
2174 //    TList* fileNoList = new TList;
2175    
2176    TIterator *listIterator = collectionList->MakeIterator();
2177    AliTPCcalibTracks *calibTracks = 0;
2178    if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;    
2179    Int_t counter = 0;
2180    while ( (calibTracks = (AliTPCcalibTracks*)listIterator->Next()) ){
2181       // loop over all entries in the collectionList and get dataMembers into lists
2182       if (!calibTracks) continue;
2183       
2184       deltaYList->Add( calibTracks->GetfDeltaY() );
2185       deltaZList->Add( calibTracks->GetfDeltaZ() );
2186       arrayAmpRowList->Add(calibTracks->GetfArrayAmpRow());
2187       arrayAmpList->Add(calibTracks->GetfArrayAmp());
2188       arrayQDYList->Add(calibTracks->GetfArrayQDY());
2189       arrayQDZList->Add(calibTracks->GetfArrayQDZ());
2190       arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
2191       arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
2192       resolYList->Add(calibTracks->GetfResolY());
2193       resolZList->Add(calibTracks->GetfResolZ());
2194       rMSYList->Add(calibTracks->GetfRMSY());
2195       rMSZList->Add(calibTracks->GetfRMSZ());
2196       arrayChargeVsDriftlengthList->Add(calibTracks->GetfArrayChargeVsDriftlength());
2197       calPadRegionChargeVsDriftlengthList->Add(calibTracks->GetCalPadRegionchargeVsDriftlength());
2198       hclusList->Add(calibTracks->GetfHclus());
2199       rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
2200       clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
2201       hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow());
2202       hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw());
2203       fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
2204       fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
2205       counter++;
2206       if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
2207    }
2208    
2209    
2210    // merge data members
2211    if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl; 
2212    fDeltaY->Merge(deltaYList);
2213    fDeltaZ->Merge(deltaZList);
2214    fHclus->Merge(hclusList);
2215    fClusterCutHisto->Merge(clusterCutHistoList);
2216    fRejectedTracksHisto->Merge(rejectedTracksList);
2217    fHclusterPerPadrow->Merge(hclusterPerPadrowList);
2218    fHclusterPerPadrowRaw->Merge(hclusterPerPadrowRawList);
2219    
2220    TObjArray* objarray = 0;
2221    TH1* hist = 0;
2222    TList* histList = 0;
2223    TIterator *objListIterator = 0;
2224    
2225    if (GetDebugLevel() > 0) cout << "merging fArrayAmpRows..." << endl;
2226    // merge fArrayAmpRows
2227    for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) {  // loop over data member, i<72
2228       objListIterator = arrayAmpRowList->MakeIterator();
2229       histList = new TList;
2230       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2231          // loop over arrayAmpRowList, get TObjArray, get object at position i, cast it into TProfile
2232          if (!objarray) continue;
2233          hist = (TProfile*)(objarray->At(i));
2234          histList->Add(hist);
2235       }
2236       ((TProfile*)(fArrayAmpRow->At(i)))->Merge(histList);
2237       delete histList;
2238       delete objListIterator;
2239    }
2240    
2241    if (GetDebugLevel() > 0) cout << "merging fArrayAmps..." << endl;
2242    // merge fArrayAmps
2243    for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) {  // loop over data member, i<72
2244       TIterator *objListIterator = arrayAmpList->MakeIterator();
2245       histList = new TList;
2246       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2247          // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F
2248          if (!objarray) continue;
2249          hist = (TH1F*)(objarray->At(i));
2250          histList->Add(hist);
2251       }
2252       ((TH1F*)(fArrayAmp->At(i)))->Merge(histList);
2253       delete histList;
2254       delete objListIterator;
2255    }
2256    
2257    if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
2258    // merge fArrayQDY
2259    for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
2260       objListIterator = arrayQDYList->MakeIterator();
2261       histList = new TList;
2262       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2263          // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
2264          if (!objarray) continue;
2265          hist = (TH3F*)(objarray->At(i));
2266          histList->Add(hist);
2267       }
2268       ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
2269       delete histList;
2270       delete objListIterator;
2271    }
2272
2273    if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
2274    // merge fArrayQDZ
2275    for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2276       objListIterator = arrayQDZList->MakeIterator();
2277       histList = new TList;
2278       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2279          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2280          if (!objarray) continue;
2281          hist = (TH3F*)(objarray->At(i));
2282          histList->Add(hist);
2283       }
2284       ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
2285       delete histList;
2286       delete objListIterator;
2287    }
2288
2289    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
2290    // merge fArrayQRMSY
2291    for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
2292       objListIterator = arrayQRMSYList->MakeIterator();
2293       histList = new TList;
2294       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2295          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2296          if (!objarray) continue;
2297          hist = (TH3F*)(objarray->At(i));
2298          histList->Add(hist);
2299       }
2300       ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
2301       delete histList;
2302       delete objListIterator;
2303    }   
2304
2305    if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
2306    // merge fArrayQRMSZ
2307    for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2308       objListIterator = arrayQRMSZList->MakeIterator();
2309       histList = new TList;
2310       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2311          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2312          if (!objarray) continue;
2313          hist = (TH3F*)(objarray->At(i));
2314          histList->Add(hist);
2315       }
2316       ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
2317       delete histList;
2318       delete objListIterator;
2319    } 
2320    
2321    if (GetDebugLevel() > 0) cout << "merging fArrayChargeVsDriftlength..." << endl;
2322    // merge fArrayChargeVsDriftlength
2323    for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { // loop over data member, i < 300
2324       objListIterator = arrayChargeVsDriftlengthList->MakeIterator();
2325       histList = new TList;
2326       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2327          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TProfile
2328          if (!objarray) continue;
2329          hist = (TProfile*)(objarray->At(i));
2330          histList->Add(hist);
2331       }
2332       ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Merge(histList);
2333       delete histList;
2334       delete objListIterator;
2335    }    
2336    
2337    if (GetDebugLevel() > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl;
2338    // merge fcalPadRegionChargeVsDriftlength
2339    AliTPCCalPadRegion *cpr = 0x0;
2340    
2341    /*
2342    TIterator *regionIterator = fcalPadRegionChargeVsDriftlength->MakeIterator();
2343    while (hist = (TProfile*)regionIterator->Next()) {
2344       // loop over all calPadRegion's in destination calibTracks object
2345          objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2346          while (( cpr =  (AliTPCCalPadRegion*)objListIterator->Next() )) { 
2347
2348       
2349       hist->Merge(...);
2350    }
2351    */
2352    
2353    for (UInt_t isec = 0; isec < 36; isec++) {
2354       for (UInt_t padSize = 0; padSize < 3; padSize++){
2355          objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2356          histList = new TList;
2357          while (( cpr =  (AliTPCCalPadRegion*)objListIterator->Next() )) { 
2358             // loop over calPadRegionChargeVsDriftlengthList, get AliTPCCalPadRegion, get object 
2359             if (!cpr) continue;
2360             hist = (TProfile*)cpr->GetObject(isec, padSize);
2361             histList->Add(hist);
2362          }
2363          ((TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isec, padSize)))->Merge(histList);
2364          delete histList;
2365          delete objListIterator;
2366       }
2367    }
2368   
2369    
2370         
2371    
2372    if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
2373    // merge fResolY
2374    for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
2375       objListIterator = resolYList->MakeIterator();
2376       histList = new TList;
2377       while (( objarray =  (TObjArray*)objListIterator->Next() )) { 
2378          // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2379          if (!objarray) continue;
2380          hist = (TH3F*)(objarray->At(i));
2381          histList->Add(hist);
2382       }
2383       ((TH3F*)(fResolY->At(i)))->Merge(histList);
2384       delete histList;
2385       delete objListIterator;
2386    }
2387    
2388    // merge fResolZ
2389    for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2390       objListIterator = resolZList->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*)(fResolZ->At(i)))->Merge(histList);
2399       delete histList;
2400       delete objListIterator;
2401    }
2402
2403    // merge fRMSY
2404    for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
2405       objListIterator = rMSYList->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*)(fRMSY->At(i)))->Merge(histList);
2414       delete histList;
2415       delete objListIterator;
2416    }
2417          
2418    // merge fRMSZ
2419    for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2420       objListIterator = rMSZList->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*)(fRMSZ->At(i)))->Merge(histList);
2429       delete histList;
2430       delete objListIterator;
2431    }
2432    
2433    delete deltaYList;
2434    delete deltaZList;
2435    delete arrayAmpRowList;
2436    delete arrayAmpList;
2437    delete arrayQDYList;
2438    delete arrayQDZList;
2439    delete arrayQRMSYList;
2440    delete arrayQRMSZList;
2441    delete resolYList;
2442    delete resolZList;
2443    delete rMSYList;
2444    delete rMSZList;
2445    delete listIterator;
2446    
2447    if (GetDebugLevel() > 0) cout << "merging done!" << endl;
2448    
2449    return 1;
2450 }
2451
2452
2453 AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks){
2454    // 
2455    // function to test AliTPCcalibTrack::Merge:
2456    // in the file 'f' is a AliTPCcalibTrack object with name "calibTracks"
2457    // this object is appended 'nCalTracks' times to a TList
2458    // A new AliTPCcalibTrack object is created which merges the list
2459    // this object is returned
2460    // 
2461    /*
2462    // .L AliTPCcalibTracks.cxx+g
2463    .L libTPCcalib.so
2464    TFile f("Output.root");
2465    AliTPCcalibTracks* calTracks = (AliTPCcalibTracks*)f.Get("calibTracks");
2466    //f.Close();
2467    TFile clusterParamFile("/u/lbozyk/calibration/workdir/calibTracks/TPCClusterParam.root");
2468    AliTPCClusterParam *clusterParam  =  (AliTPCClusterParam *) clusterParamFile.Get("Param"); 
2469    clusterParamFile.Close();
2470
2471    AliTPCcalibTracks::TestMerge(calTracks, clusterParam);
2472    */
2473    TList *list = new TList();
2474    if (ct == 0 || clusterParam == 0) return 0;
2475    cout << "making list with " << nCalTracks << " AliTPCcalibTrack objects" << endl;
2476    for (Int_t i = 0; i < nCalTracks; i++) {
2477       if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl;
2478       list->Add(new AliTPCcalibTracks(*ct));
2479    }
2480    
2481    // only for check at the end
2482    AliTPCcalibTracks* cal1 = new AliTPCcalibTracks(*ct);
2483    Double_t cal1Entries = ((TH1F*)cal1->GetfArrayAmpRow()->At(5))->GetEntries();
2484 //    Double_t cal1Entries = 5; //((TH1F*)ct->GetfArrayAmpRow()->At(5))->GetEntries();
2485
2486    cout  << "The list contains " << list->GetEntries() << " entries. " << endl;
2487   
2488    
2489    AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018);
2490    AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts, 5);
2491    cal->Merge(list);
2492    
2493    cout << "cal->GetfArrayAmpRow()->At(5)->Print():" << endl;
2494    cal->GetfArrayAmpRow()->At(5)->Print();
2495    Double_t calEntries = ((TH1F*)cal->GetfArrayAmpRow()->At(5))->GetEntries();
2496    
2497    cout << "cal1->GetfArrayAmpRow()->At(5))->GetEntries() = " << cal1Entries << endl;
2498    cout << " cal->GetfArrayAmpRow()->At(5))->GetEntries() = " <<  calEntries << endl;
2499    printf("That means there were %f / %f = %f AliTPCcalibTracks-Objects merged. \n", 
2500       calEntries, cal1Entries, ((Double_t)calEntries/cal1Entries));
2501    
2502    return cal;
2503
2504 }
2505
2506
2507 void  AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints, Int_t verbose){
2508   //
2509   // Make position corrections
2510   // for the moment Only using debug streamer 
2511   // chainres  - debug tree
2512   // param     - parameters to be updated
2513   // maxPoints - maximal number of points using for fit
2514   // verbose   - print info flag
2515   //
2516   // Current implementation - only using debug streamers
2517   // 
2518   
2519   /*    
2520     //Defaults
2521     Int_t maxPoints=100000;
2522   */
2523   /*
2524     Usage: 
2525     //0. Load libraries
2526     gSystem->Load("libANALYSIS");
2527     gSystem->Load("libSTAT");
2528     gSystem->Load("libTPCcalib");
2529     
2530
2531     //1. Load Parameters to be modified:
2532     //e.g:
2533     
2534     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
2535     AliCDBManager::Instance()->SetRun(0) 
2536     AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
2537
2538     //2. Load chain from debug streamers
2539     //
2540     //e.g
2541     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2542     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2543     AliXRDPROOFtoolkit tool;  
2544     TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
2545     chainres->Lookup();
2546     //3. Do fits and store results
2547     // 
2548     AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ;
2549     TFile f("paramout.root","recreate");
2550     param->Write("clusterParam");
2551     f.Close();
2552     //4. Verification
2553     TFile f2("paramout.root");
2554     AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
2555     param2->SetInstance(param2);
2556     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 
2557     
2558    */
2559
2560
2561   TStatToolkit toolkit;
2562   Double_t chi2;
2563   TVectorD fitParamY0;
2564   TVectorD fitParamY1;
2565   TVectorD fitParamZ0;
2566   TVectorD fitParamZ1;
2567   TMatrixD covMatrix;
2568   Int_t npoints;
2569   
2570   chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
2571   chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
2572   chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
2573   chainres->SetAlias("st","(sin(dt)-dt)");
2574   //
2575   chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
2576   chainres->SetAlias("dq","sqrt(15./(5+Cl.fMax))");
2577   chainres->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
2578   chainres->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
2579   //
2580   //
2581   //
2582   TCut cutA("1");
2583   TString fstringY="";  
2584   //
2585   fstringY+="(dp)++";            //1
2586   fstringY+="(dp)*di++";         //2
2587   fstringY+="(dp)*dq++";         //3
2588   fstringY+="(dp)*sy++";         //4
2589   //
2590   fstringY+="(sp)++";            //5
2591   fstringY+="(sp)*di++";         //6
2592   fstringY+="(sp)*dq++";         //7
2593   fstringY+="(sp)*sy++";         //8
2594   //
2595
2596   TString fstringZ="";  
2597   fstringZ+="(dt)++";            //1
2598   fstringZ+="(dt)*di++";         //2
2599   fstringZ+="(dt)*dq++";         //3
2600   fstringZ+="(dt)*sz++";         //4
2601   //
2602   fstringZ+="(st)++";            //5
2603   fstringZ+="(st)*di++";         //6
2604   fstringZ+="(st)*dq++";         //7
2605   fstringZ+="(st)*sz++";         //8
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   gSystem->Load("libSTAT.so");
2643   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2644   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2645   AliXRDPROOFtoolkit tool;  
2646   TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
2647   chainres->Lookup();
2648   //
2649   //
2650   //
2651
2652
2653
2654   TStatToolkit toolkit;
2655   Double_t chi2;
2656   TVectorD fitParamY0;
2657   TVectorD fitParamY1;
2658   TVectorD fitParamZ0;
2659   TVectorD fitParamZ1;
2660   TMatrixD covMatrix;
2661   Int_t npoints;
2662
2663   chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
2664   chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
2665   chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
2666   chainres->SetAlias("st","(sin(dt)-dt)");
2667   //
2668   chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
2669   chainres->SetAlias("dq","sqrt(15./(5+Cl.fMax))");
2670   chainres->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
2671   chainres->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
2672   //
2673   //
2674   //
2675   TCut cutA("1")
2676   
2677   TString fstringY="";  
2678   //
2679   fstringY+="(dp)++";            //1
2680   fstringY+="(dp)*di++";         //2
2681   fstringY+="(dp)*dq++";         //3
2682   fstringY+="(dp)*sy++";         //4
2683   //
2684   fstringY+="(sp)++";            //5
2685   fstringY+="(sp)*di++";         //6
2686   fstringY+="(sp)*dq++";         //7
2687   fstringY+="(sp)*sy++";         //8
2688   //
2689
2690   TString fstringZ="";  
2691   fstringZ+="(dt)++";            //1
2692   fstringZ+="(dt)*di++";         //2
2693   fstringZ+="(dt)*dq++";         //3
2694   fstringZ+="(dt)*sz++";         //4
2695   //
2696   fstringZ+="(st)++";            //5
2697   fstringZ+="(st)*di++";         //6
2698   fstringZ+="(st)*dq++";         //7
2699   fstringZ+="(st)*sz++";         //8
2700
2701   
2702   TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,100000);
2703   TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,100000);
2704
2705   TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,100000);
2706   TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,100000);
2707
2708   chainres->SetAlias("fitZ0",strZ0.Data());
2709   chainres->SetAlias("fitZ1",strZ1.Data());
2710   chainres->SetAlias("fitY0",strY0.Data());
2711   chainres->SetAlias("fitY1",strY1.Data());
2712   chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000)
2713
2714
2715
2716 */
2717
2718
2719