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