1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
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. //
33 Raw Data -> Local Reconstruction -> Tracking -> Calibration -> RefData (component itself)
34 Offline/HLT Offline/HLT OCDB entries (AliTPCClusterParam)
39 How to retrive it from file (created using calibration task):
41 gSystem->Load("libANALYSIS");
42 gSystem->Load("libTPCcalib");
43 TFile fcalib("CalibObjects.root");
44 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
45 AliTPCcalibTracks * calibTracks = ( AliTPCcalibTracks *)array->FindObject("calibTracks");
48 //USAGE of debug stream example
49 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
50 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
51 AliXRDPROOFtoolkit tool;
52 TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
58 ///////////////////////////////////////////////////////////////////////////////
74 //#include <TPDGCode.h>
76 #include "TLinearFitter.h"
77 //#include "TMatrixD.h"
78 #include "TTreeStream.h"
81 #include <TGraph2DErrors.h>
82 #include "TPostScript.h"
88 #include <TCollection.h>
90 #include <TLinearFitter.h>
97 #include "AliTracker.h"
99 #include "AliESDtrack.h"
100 #include "AliESDfriend.h"
101 #include "AliESDfriendTrack.h"
102 #include "AliTPCseed.h"
103 #include "AliTPCclusterMI.h"
104 #include "AliTPCROC.h"
106 #include "AliTPCParamSR.h"
107 #include "AliTrackPointArray.h"
108 #include "AliTPCcalibTracks.h"
109 #include "AliTPCClusterParam.h"
110 #include "AliTPCcalibTracksCuts.h"
111 #include "AliTPCCalPadRegion.h"
112 #include "AliTPCCalPad.h"
113 #include "AliTPCCalROC.h"
115 #include "TPaveText.h"
117 #include "TStatToolkit.h"
122 ClassImp(AliTPCcalibTracks)
125 AliTPCcalibTracks::AliTPCcalibTracks():
135 fArrayChargeVsDriftlength(0),
136 fcalPadRegionChargeVsDriftlength(0),
145 fRejectedTracksHisto(0),
146 fHclusterPerPadrow(0),
147 fHclusterPerPadrowRaw(0),
149 fCalPadClusterPerPad(0),
150 fCalPadClusterPerPadRaw(0)
153 // AliTPCcalibTracks default constructor
156 if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl;
161 AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
162 AliTPCcalibBase(calibTracks),
171 fArrayChargeVsDriftlength(0),
172 fcalPadRegionChargeVsDriftlength(0),
181 fRejectedTracksHisto(0),
182 fHclusterPerPadrow(0),
183 fHclusterPerPadrowRaw(0),
185 fCalPadClusterPerPad(0),
186 fCalPadClusterPerPadRaw(0)
189 // AliTPCcalibTracks copy constructor
191 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl;
193 Bool_t dirStatus = TH1::AddDirectoryStatus();
194 TH1::AddDirectory(kFALSE);
197 // backward compatibility: if the data member doesn't yet exist, it will not be merged
198 (calibTracks.fArrayAmpRow) ? length = calibTracks.fArrayAmpRow->GetEntriesFast() : length = -1;
199 fArrayAmpRow = new TObjArray(length);
200 fArrayAmp = new TObjArray(length);
201 for (Int_t i = 0; i < length; i++) {
202 fArrayAmpRow->AddAt( (TProfile*)calibTracks.fArrayAmpRow->At(i)->Clone(), i);
203 fArrayAmp->AddAt( ((TProfile*)calibTracks.fArrayAmp->At(i)->Clone()), i);
206 (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
207 fArrayQDY= new TObjArray(length);
208 fArrayQDZ= new TObjArray(length);
209 fArrayQRMSY= new TObjArray(length);
210 fArrayQRMSZ= new TObjArray(length);
211 for (Int_t i = 0; i < length; i++) {
212 fArrayQDY->AddAt( ((TH1F*)calibTracks.fArrayQDY->At(i)->Clone()), i);
213 fArrayQDZ->AddAt( ((TH1F*)calibTracks.fArrayQDZ->At(i)->Clone()), i);
214 fArrayQRMSY->AddAt( ((TH1F*)calibTracks.fArrayQRMSY->At(i)->Clone()), i);
215 fArrayQRMSZ->AddAt( ((TH1F*)calibTracks.fArrayQRMSZ->At(i)->Clone()), i);
218 (calibTracks.fResolY) ? length = calibTracks.fResolY->GetEntriesFast() : length = -1;
219 fResolY = new TObjArray(length);
220 fResolZ = new TObjArray(length);
221 fRMSY = new TObjArray(length);
222 fRMSZ = new TObjArray(length);
223 for (Int_t i = 0; i < length; i++) {
224 fResolY->AddAt( ((TH1F*)calibTracks.fResolY->At(i)->Clone()), i);
225 fResolZ->AddAt( ((TH1F*)calibTracks.fResolZ->At(i)->Clone()), i);
226 fRMSY->AddAt( ((TH1F*)calibTracks.fRMSY->At(i)->Clone()), i);
227 fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
230 (calibTracks.fArrayChargeVsDriftlength) ? length = calibTracks.fArrayChargeVsDriftlength->GetEntriesFast() : length = -1;
231 (calibTracks.fArrayChargeVsDriftlength) ? fArrayChargeVsDriftlength = new TObjArray(length) : fArrayChargeVsDriftlength = 0;
232 for (Int_t i = 0; i < length; i++) {
233 fArrayChargeVsDriftlength->AddAt( ((TProfile*)calibTracks.fArrayChargeVsDriftlength->At(i)->Clone()), i);
236 fDeltaY = (TH1F*)calibTracks.fDeltaY->Clone();
237 fDeltaZ = (TH1F*)calibTracks.fDeltaZ->Clone();
238 fHclus = (TH1I*)calibTracks.fHclus->Clone();
239 fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
240 fRejectedTracksHisto = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
241 fHclusterPerPadrow = (TH1I*)calibTracks.fHclusterPerPadrow->Clone();
242 fHclusterPerPadrowRaw = (TH1I*)calibTracks.fHclusterPerPadrowRaw->Clone();
243 fcalPadRegionChargeVsDriftlength = (AliTPCCalPadRegion*)calibTracks.fcalPadRegionChargeVsDriftlength->Clone();
244 fCalPadClusterPerPad = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
245 fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
247 fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(),
248 calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise());
249 SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle());
250 TH1::AddDirectory(dirStatus); // set status back to original status
251 // cout << "+++++ end of copy constructor +++++" << endl; // TO BE REMOVED
255 AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibTracks){
257 // assgnment operator
259 if (this != &calibTracks) {
260 new (this) AliTPCcalibTracks(calibTracks);
267 AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam, AliTPCcalibTracksCuts* cuts, Int_t logLevel) :
277 fArrayChargeVsDriftlength(0),
278 fcalPadRegionChargeVsDriftlength(0),
287 fRejectedTracksHisto(0),
288 fHclusterPerPadrow(0),
289 fHclusterPerPadrowRaw(0),
291 fCalPadClusterPerPad(0),
292 fCalPadClusterPerPadRaw(0)
295 // AliTPCcalibTracks constructor
296 // specify 'name' and 'title' of your object
297 // specify 'clusterParam', (needed for TPC cluster error and shape parameterization)
298 // In the parameter 'cuts' the cuts are specified, that decide
299 // weather a track will be accepted for calibration or not.
301 // fDebugLevel - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStreamer, 6: waste your screen
303 // All histograms are instatiated in this constructor.
306 this->SetTitle(title);
308 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl;
309 G__SetCatchException(0);
311 fClusterParam = clusterParam;
313 fClusterParam->SetInstance(fClusterParam);
316 Error("AliTPCcalibTracks","No cluster parametrization found! A valid clusterParam object is needed in the constructor. (To be found in 'TPCClusterParam.root'.)");
319 SetDebugLevel(logLevel);
321 TH1::AddDirectory(kFALSE);
326 fHclus = new TH1I("hclus","Number of clusters per track",160, 0, 160); // valgrind 3
327 fRejectedTracksHisto = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
328 fHclusterPerPadrow = new TH1I("fHclusterPerPadrow", " clusters per padRow, used for the resolution tree", 160, 0, 160);
329 fHclusterPerPadrowRaw = new TH1I("fHclusterPerPadrowRaw", " clusters per padRow, before cutting clusters", 160, 0, 160);
330 fCalPadClusterPerPad = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
331 fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
332 fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
334 // Amplitude - sector - row histograms
335 fArrayAmpRow = new TObjArray(72);
336 fArrayAmp = new TObjArray(72);
337 fArrayChargeVsDriftlength = new TObjArray(72);
339 for (Int_t i = 0; i < 36; i++){
340 sprintf(chname,"Amp_row_Sector%d",i);
341 prof1 = new TProfile(chname,chname,63,0,64); // valgrind 3 193,536 bytes in 354 blocks are still reachable
342 prof1->SetXTitle("Pad row");
343 prof1->SetYTitle("Mean Max amplitude");
344 fArrayAmpRow->AddAt(prof1,i);
345 sprintf(chname,"Amp_row_Sector%d",i+36);
346 prof1 = new TProfile(chname,chname,96,0,97); // valgrind 3 3,912 bytes in 6 blocks are possibly lost
347 prof1->SetXTitle("Pad row");
348 prof1->SetYTitle("Mean Max amplitude");
349 fArrayAmpRow->AddAt(prof1,i+36);
352 sprintf(chname,"Amp_Sector%d",i);
353 his1 = new TH1F(chname,chname,250,0,500); // valgrind
354 his1->SetXTitle("Max Amplitude (ADC)");
355 fArrayAmp->AddAt(his1,i);
356 sprintf(chname,"Amp_Sector%d",i+36);
357 his1 = new TH1F(chname,chname,200,0,600); // valgrind 3 13,408,208 bytes in 229 blocks are still reachable
358 his1->SetXTitle("Max Amplitude (ADC)");
359 fArrayAmp->AddAt(his1,i+36);
362 sprintf(chname, "driftlengt vs. charge, ROC %i", i);
363 prof1 = new TProfile(chname, chname, 500, 0, 250);
364 prof1->SetYTitle("Charge");
365 prof1->SetXTitle("Driftlength");
366 fArrayChargeVsDriftlength->AddAt(prof1,i);
367 sprintf(chname, "driftlengt vs. charge, ROC %i", i+36);
368 prof1 = new TProfile(chname, chname, 500, 0, 250);
369 prof1->SetYTitle("Charge");
370 prof1->SetXTitle("Driftlength");
371 fArrayChargeVsDriftlength->AddAt(prof1,i+36);
374 TH1::AddDirectory(kFALSE);
376 fDeltaY = new TH1F("DeltaY","DeltaY",100,-1,1);
377 fDeltaZ = new TH1F("DeltaZ","DeltaZ",100,-1,1);
379 fResolY = new TObjArray(3);
380 fResolZ = new TObjArray(3);
381 fRMSY = new TObjArray(3);
382 fRMSZ = new TObjArray(3);
385 his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
386 fResolY->AddAt(his3D,0);
387 his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
388 fResolY->AddAt(his3D,1);
389 his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
390 fResolY->AddAt(his3D,2);
392 his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
393 fResolZ->AddAt(his3D,0);
394 his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
395 fResolZ->AddAt(his3D,1);
396 his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
397 fResolZ->AddAt(his3D,2);
399 his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
400 fRMSY->AddAt(his3D,0);
401 his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
402 fRMSY->AddAt(his3D,1);
403 his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
404 fRMSY->AddAt(his3D,2);
406 his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
407 fRMSZ->AddAt(his3D,0);
408 his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
409 fRMSZ->AddAt(his3D,1);
410 his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
411 fRMSZ->AddAt(his3D,2);
414 TH1::AddDirectory(kFALSE);
416 fArrayQDY = new TObjArray(300);
417 fArrayQDZ = new TObjArray(300);
418 fArrayQRMSY = new TObjArray(300);
419 fArrayQRMSZ = new TObjArray(300);
420 for (Int_t iq = 0; iq <= 10; iq++){
421 for (Int_t ipad = 0; ipad < 3; ipad++){
422 Int_t bin = GetBin(iq, ipad);
423 Float_t qmean = GetQ(bin);
425 sprintf(name,"ResolY Pad%d Qmiddle%f",ipad, qmean);
426 his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1);
427 fArrayQDY->AddAt(his3D, bin);
428 sprintf(name,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
429 his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1);
430 fArrayQDZ->AddAt(his3D, bin);
431 sprintf(name,"RMSY Pad%d Qmiddle%f",ipad, qmean);
432 his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1);
433 fArrayQRMSY->AddAt(his3D, bin);
434 sprintf(name,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
435 his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1);
436 fArrayQRMSZ->AddAt(his3D, bin);
440 fcalPadRegionChargeVsDriftlength = new AliTPCCalPadRegion("fcalPadRegionChargeVsDriftlength", "TProfiles with charge vs driftlength for each pad region");
442 for (UInt_t padSize = 0; padSize < 3; padSize++) {
443 for (UInt_t isector = 0; isector < 36; isector++) {
444 if (padSize == 0) sprintf(chname, "driftlengt vs. charge, sector %i, short pads", isector);
445 if (padSize == 1) sprintf(chname, "driftlengt vs. charge, sector %i, medium pads", isector);
446 if (padSize == 2) sprintf(chname, "driftlengt vs. charge, sector %i, long pads", isector);
447 tempProf = new TProfile(chname, chname, 500, 0, 250);
448 tempProf->SetYTitle("Charge");
449 tempProf->SetXTitle("Driftlength");
450 fcalPadRegionChargeVsDriftlength->SetObject(tempProf, isector, padSize);
455 if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl;
456 cout << "end of main constructor" << endl; // TO BE REMOVED
460 AliTPCcalibTracks::~AliTPCcalibTracks() {
462 // AliTPCcalibTracks destructor
465 if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
467 if (fArrayAmpRow) length = fArrayAmpRow->GetEntriesFast();
468 for (Int_t i = 0; i < length; i++){
469 delete fArrayAmpRow->At(i);
470 delete fArrayAmp->At(i);
478 if (fResolY) length = fResolY->GetEntriesFast();
479 for (Int_t i = 0; i < length; i++){
480 delete fResolY->At(i);
481 delete fResolZ->At(i);
490 if (fArrayQDY) length = fArrayQDY->GetEntriesFast();
491 for (Int_t i = 0; i < length; i++){
492 delete fArrayQDY->At(i);
493 delete fArrayQDZ->At(i);
494 delete fArrayQRMSY->At(i);
495 delete fArrayQRMSZ->At(i);
498 if (fArrayChargeVsDriftlength) length = fArrayChargeVsDriftlength->GetEntriesFast();
499 for (Int_t i = 0; i < length; i++){
500 delete fArrayChargeVsDriftlength->At(i);
508 delete fArrayChargeVsDriftlength;
511 delete fRejectedTracksHisto;
512 delete fClusterCutHisto;
513 delete fHclusterPerPadrow;
514 delete fHclusterPerPadrowRaw;
515 if (fCalPadClusterPerPad) delete fCalPadClusterPerPad;
516 if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
517 if(fcalPadRegionChargeVsDriftlength) {
518 fcalPadRegionChargeVsDriftlength->Delete();
519 delete fcalPadRegionChargeVsDriftlength;
525 void AliTPCcalibTracks::Process(AliTPCseed *track){
527 // To be called in the selector
528 // first AcceptTrack is evaluated, then calls all the following analyse functions:
529 // FillResolutionHistoLocal(track)
530 // AlignUpDown(track, esd)
532 if (GetDebugLevel() > 5) Info("Process","Starting to process the track...");
533 Int_t accpetStatus = AcceptTrack(track);
534 if (accpetStatus == 0) {
535 FillResolutionHistoLocal(track);
536 // AlignUpDown(track, esd);
538 else fRejectedTracksHisto->Fill(accpetStatus);
543 Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){
545 // calculate bins for given q and pad type
548 Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 );
555 Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){
557 // calculate bins for given iq and pad type
560 return iq * 3 + pad;;
564 Float_t AliTPCcalibTracks::GetQ(Int_t bin){
566 // returns to bin belonging charge
569 Int_t bin0 = bin / 3;
575 Float_t AliTPCcalibTracks::GetPad(Int_t bin){
577 // returns to bin belonging pad
585 Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
587 // Function, that decides wheather a given track is accepted for
588 // the analysis or not.
589 // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts'
590 // Returns 0 if a track is accepted or an integer different from 0
591 // to indicate the failed cut
593 const Int_t kMinClusters = fCuts->GetMinClusters();
594 const Float_t kMinRatio = fCuts->GetMinRatio();
595 const Float_t kMax1pt = fCuts->GetMax1pt();
596 const Float_t kEdgeYXCutNoise = fCuts->GetEdgeYXCutNoise();
597 const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise();
600 // edge induced noise tracks - NEXT RELEASE will be removed during tracking
601 if ( TMath::Abs(track->GetY() / track->GetX()) > kEdgeYXCutNoise )
602 if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return 1;
603 if (track->GetNumberOfClusters() < kMinClusters) return 2;
604 Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.);
605 if (ratio < kMinRatio) return 3;
606 // Float_t mpt = track->Get1Pt(); // Get1Pt() doesn't exist any more
607 Float_t mpt = track->GetSigned1Pt();
608 if (TMath::Abs(mpt) > kMax1pt) return 4;
609 //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
610 //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
611 //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
613 if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");
618 void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
620 // fill resolution histograms - localy - tracklet in the neighborhood
621 // write debug information to 'TPCSelectorDebug.root'
623 // _ the main function, called during track analysis _
625 // loop over all padrows along the track
626 // fit tracklets (length: 13 clusters) calculate mean chi^2 for this track-fit in Y and Z direction
628 // loop again over all padrows along the track
629 // fit tracklet (clusters in given padrow +- kDelta padrows)
630 // with polynom of 2nd order and two polynoms of 1st order
631 // take both polynoms of 1st order, calculate difference of their parameters
632 // add covariance matrixes and calculate chi2 of this difference
633 // if this chi2 is bigger than a given threshold, assume that the current cluster is
634 // a kink an goto next padrow
636 // fill fArrayAmpRow, array with amplitudes vs. row for given sector
637 // fill fArrayAmp, array with amplitude histograms for give sector
638 // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fDeltaY, fDeltaZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
640 // write debug information to 'TPCSelectorDebug.root'
641 // only for every kDeltaWriteDebugStream'th padrow to reduce data volume
642 // and to avoid redundant data
645 static TLinearFitter fFitterLinY1(2,"pol1"); //
646 static TLinearFitter fFitterLinZ1(2,"pol1"); //
647 static TLinearFitter fFitterLinY2(2,"pol1"); //
648 static TLinearFitter fFitterLinZ2(2,"pol1"); //
649 static TLinearFitter fFitterParY(3,"pol2"); //
650 static TLinearFitter fFitterParZ(3,"pol2"); //
652 fFitterLinY1.StoreData(kFALSE);
653 fFitterLinZ1.StoreData(kFALSE);
654 fFitterLinY2.StoreData(kFALSE);
655 fFitterLinZ2.StoreData(kFALSE);
656 fFitterParY.StoreData(kFALSE);
657 fFitterParZ.StoreData(kFALSE);
660 if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
661 const Int_t kDelta = 10; // delta rows to fit
662 const Float_t kMinRatio = 0.75; // minimal ratio
663 // const Float_t kCutChi2 = 6.; // cut chi2 - left right - kink removal
664 const Float_t kErrorFraction = 0.5; // use only clusters with small interpolation error - for error param
665 const Int_t kFirstLargePad = 127; // medium pads -> long pads
666 const Float_t kLargePadSize = 1.5; // factor between medium and long pads' area
667 const Int_t kDeltaWriteDebugStream = 5; // only for every kDeltaWriteDebugStream'th padrow debug information is calulated and written to debugstream
674 TMatrixD matrixY0(2,2);
675 TMatrixD matrixZ0(2,2);
676 TMatrixD matrixY1(2,2);
677 TMatrixD matrixZ1(2,2);
679 // estimate mean error
680 Int_t nTrackletsAll = 0; // number of tracklets for given track
681 Float_t csigmaY = 0; // mean sigma for tracklet refit in Y direction
682 Float_t csigmaZ = 0; // mean sigma for tracklet refit in Z direction
683 Int_t nClusters = 0; // working variable, number of clusters per tracklet
684 Int_t sectorG = -1; // working variable, sector of tracklet, has to stay constant for one tracklet
686 fHclus->Fill(track->GetNumberOfClusters()); // for statistics overview
687 // ---------------------------------------------------------------------
688 for (Int_t irow = 0; irow < 159; irow++){
689 // loop over all rows along the track
690 // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
691 // calculate mean chi^2 for this track-fit in Y and Z direction
692 AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
693 if (!cluster0) continue; // no cluster found
694 Int_t sector = cluster0->GetDetector();
695 fHclusterPerPadrowRaw->Fill(irow);
697 Int_t ipad = TMath::Nint(cluster0->GetPad());
698 Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
699 fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
701 if (sector != sectorG){
702 // track leaves sector before it crossed enough rows to fit / initialization
704 fFitterParY.ClearPoints();
705 fFitterParZ.ClearPoints();
710 Double_t x = cluster0->GetX();
711 fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
712 fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
714 if ( nClusters >= kDelta + 3 ){
715 // if more than 13 (kDelta+3) clusters were added to the fitters
716 // fit the tracklet, increase trackletCounter
720 csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
721 csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
723 fFitterParY.ClearPoints();
724 fFitterParZ.ClearPoints();
727 } // for (Int_t irow = 0; irow < 159; irow++)
728 // mean chi^2 for all tracklet fits in Y and in Z direction:
729 csigmaY = TMath::Sqrt(csigmaY / nTrackletsAll);
730 csigmaZ = TMath::Sqrt(csigmaZ / nTrackletsAll);
731 // ---------------------------------------------------------------------
733 for (Int_t irow = 0; irow < 159; irow++){
734 // loop again over all rows along the track
737 Int_t nclFound = 0; // number of clusters in the neighborhood
738 Int_t ncl0 = 0; // number of clusters in rows < rowOfCenterCluster
739 Int_t ncl1 = 0; // number of clusters in rows > rowOfCenterCluster
740 AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
741 if (!cluster0) continue;
742 Int_t sector = cluster0->GetDetector();
743 Float_t xref = cluster0->GetX();
746 fFitterParY.ClearPoints();
747 fFitterParZ.ClearPoints();
748 fFitterLinY1.ClearPoints();
749 fFitterLinZ1.ClearPoints();
750 fFitterLinY2.ClearPoints();
751 fFitterLinZ2.ClearPoints();
753 // fit tracklet (clusters in given padrow +- kDelta padrows)
754 // with polynom of 2nd order and two polynoms of 1st order
755 // take both polynoms of 1st order, calculate difference of their parameters
756 // add covariance matrixes and calculate chi2 of this difference
757 // if this chi2 is bigger than a given threshold, assume that the current cluster is
758 // a kink an goto next padrow
760 for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
761 // loop over irow +- kDelta rows (neighboured rows)
764 if (idelta == 0) continue; // don't use center cluster
765 if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
766 AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
767 if (!currentCluster) continue;
768 if (currentCluster->GetType() < 0) continue;
769 if (currentCluster->GetDetector() != sector) continue;
770 Double_t x = currentCluster->GetX() - xref; // x = differece: current cluster - cluster @ irow
774 fFitterLinY1.AddPoint(&x, currentCluster->GetY(), csigmaY);
775 fFitterLinZ1.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
779 fFitterLinY2.AddPoint(&x, currentCluster->GetY(), csigmaY);
780 fFitterLinZ2.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
782 fFitterParY.AddPoint(&x, currentCluster->GetY(), csigmaY);
783 fFitterParZ.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
784 } // loop over neighbourhood for fitter filling
788 if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
789 if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
790 if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters (7.5) found in neighbourhood goto next padrow
793 Double_t chi2 = (fFitterParY.GetChisquare() + fFitterParZ.GetChisquare()) / (2. * nclFound - 6.);
794 //if (chi2 > kCutChi2) fRejectedTracksHisto->Fill(9);
795 //if (chi2 > kCutChi2) fClusterCutHisto->Fill(2, irow);
796 //if (chi2 > kCutChi2) continue; // if chi^2 is too big goto next padrow
797 TTreeSRedirector *cstream = GetDebugStreamer();
804 // only when there are enough clusters (4) in each direction
814 if (ncl0 > 4 && ncl1 > 4){
815 fFitterLinY1.GetCovarianceMatrix(matrixY0);
816 fFitterLinY2.GetCovarianceMatrix(matrixY1);
817 fFitterLinZ1.GetCovarianceMatrix(matrixZ0);
818 fFitterLinZ2.GetCovarianceMatrix(matrixZ1);
819 fFitterLinY2.GetParameters(paramY1);
820 fFitterLinZ2.GetParameters(paramZ1);
821 fFitterLinY1.GetParameters(paramY0);
822 fFitterLinZ1.GetParameters(paramZ0);
825 matrixY0 += matrixY1;
826 matrixZ0 += matrixZ1;
829 TMatrixD difY(2, 1, paramY0.GetMatrixArray());
830 TMatrixD difYT(1, 2, paramY0.GetMatrixArray());
832 TMatrixD mulY(matrixY0, TMatrixD::kMult, difY);
833 TMatrixD chi2Y(difYT, TMatrixD::kMult, mulY);
836 TMatrixD difZ(2, 1, paramZ0.GetMatrixArray());
837 TMatrixD difZT(1, 2, paramZ0.GetMatrixArray());
839 TMatrixD mulZ(matrixZ0, TMatrixD::kMult, difZ);
840 TMatrixD chi2Z(difZT, TMatrixD::kMult, mulZ);
843 // REMOVE KINK - TO be fixed - proper chi2 calculation for curved track to be implemented
844 //if (chi2 * 0.25 > kCutChi2) fRejectedTracksHisto->Fill(8);
845 //if (chi2 * 0.25 > kCutChi2) fClusterCutHisto->Fill(3, irow);
846 //if (chi2 * 0.25 > kCutChi2) continue; // if chi2 is too big goto next padrow
847 // fit tracklet with polynom of 2nd order and two polynoms of 1st order
848 // take both polynoms of 1st order, calculate difference of their parameters
849 // add covariance matrixes and calculate chi2 of this difference
850 // if this chi2 is bigger than a given threshold, assume that the current cluster is
851 // a kink an goto next padrow
853 TTreeSRedirector *cstream = GetDebugStreamer();
861 // current padrow has no kink
863 // get fit parameters from pol2 fit:
864 Double_t paramY[4], paramZ[4];
865 paramY[0] = fFitterParY.GetParameter(0);
866 paramY[1] = fFitterParY.GetParameter(1);
867 paramY[2] = fFitterParY.GetParameter(2);
868 paramZ[0] = fFitterParZ.GetParameter(0);
869 paramZ[1] = fFitterParZ.GetParameter(1);
870 paramZ[2] = fFitterParZ.GetParameter(2);
872 Double_t tracky = paramY[0];
873 Double_t trackz = paramZ[0];
874 Float_t deltay = tracky - cluster0->GetY();
875 Float_t deltaz = trackz - cluster0->GetZ();
876 Float_t angley = paramY[1] - paramY[0] / xref;
877 Float_t anglez = paramZ[1];
879 Float_t max = cluster0->GetMax();
880 UInt_t isegment = cluster0->GetDetector() % 36;
881 Int_t padSize = 0; // short pads
882 if (cluster0->GetDetector() >= 36) {
883 padSize = 1; // medium pads
884 if (cluster0->GetRow() > 63) padSize = 2; // long pads
887 // =========================================
888 // wirte collected information to histograms
889 // =========================================
891 TProfile *profAmpRow = (TProfile*)fArrayAmpRow->At(sector);
892 if ( irow >= kFirstLargePad) max /= kLargePadSize;
893 profAmpRow->Fill( (Double_t)cluster0->GetRow(), max );
894 TH1F *hisAmp = (TH1F*)fArrayAmp->At(sector);
897 // remove the following two lines one day:
898 TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector);
899 profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max );
901 TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize));
902 profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max );
904 fHclusterPerPadrow->Fill(irow); // fill histogram showing clusters per padrow
905 Int_t ipad = TMath::Nint(cluster0->GetPad());
906 Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
907 fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
911 his3 = (TH3F*)fRMSY->At(padSize);
912 if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
913 his3 = (TH3F*)fRMSZ->At(padSize);
914 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
916 his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
917 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
918 his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
919 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
922 // Fill resolution histograms
923 Bool_t useForResol = kTRUE;
924 if (fFitterParY.GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE;
927 Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
928 Float_t sy = cluster0->GetSigmaY2();
929 Float_t sz = cluster0->GetSigmaZ2();
930 (*cstream)<<"Resol0"<<
931 "padSize="<<padSize<<
943 fDeltaY->Fill(deltay);
944 fDeltaZ->Fill(deltaz);
945 his3 = (TH3F*)fResolY->At(padSize);
946 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
947 his3 = (TH3F*)fResolZ->At(padSize);
948 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
949 his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
950 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
951 his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
952 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
955 //=============================================================================================
957 if (useForResol && nclFound > 2 * kMinRatio * kDelta
958 && irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){
959 if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
960 FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta);
961 } // if (useForResol && nclFound > 2 * kMinRatio * kDelta)
963 } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
964 } // FillResolutionHistoLocal(...)
968 void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t angley, Float_t anglez, Int_t nclFound, Int_t kDelta) {
970 // - debug part of FillResolutionHistoLocal -
971 // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data
972 // called only for GetStreamLevel() > 4
973 // fill resolution trees
976 Int_t sector = cluster0->GetDetector();
977 Float_t xref = cluster0->GetX();
978 Int_t padSize = 0; // short pads
979 if (cluster0->GetDetector() >= 36) {
980 padSize = 1; // medium pads
981 if (cluster0->GetRow() > 63) padSize = 2; // long pads
984 static TLinearFitter fitY0(3, "pol2");
985 static TLinearFitter fitZ0(3, "pol2");
986 static TLinearFitter fitY2(5, "hyp4");
987 static TLinearFitter fitZ2(5, "hyp4");
988 static TLinearFitter fitY2Q(5, "hyp4");
989 static TLinearFitter fitZ2Q(5, "hyp4");
990 static TLinearFitter fitY2S(5, "hyp4");
991 static TLinearFitter fitZ2S(5, "hyp4");
996 fitY2Q.ClearPoints();
997 fitZ2Q.ClearPoints();
998 fitY2S.ClearPoints();
999 fitZ2S.ClearPoints();
1001 for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
1002 // loop over irow +- kDelta rows (neighboured rows)
1005 if (idelta == 0) continue;
1006 if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
1007 AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta);
1008 if (!cluster) continue;
1009 if (cluster->GetType() < 0) continue;
1010 if (cluster->GetDetector() != sector) continue;
1011 Double_t x = cluster->GetX() - xref;
1012 Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) );
1013 Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) );
1015 Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1016 Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1017 Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
1018 TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1019 Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
1020 TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1021 Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1022 TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1023 TMath::Sqrt(cluster0->GetSigmaY2()), 0 );
1024 Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1025 TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1026 TMath::Sqrt(cluster0->GetSigmaZ2()),0 );
1027 sigmaYS = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.);
1028 sigmaZS = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.);
1031 fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
1032 fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
1035 xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0;
1037 xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0;
1039 fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
1040 fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
1041 fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
1042 fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
1043 fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
1044 fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
1046 } // neigbouhood-loop
1056 Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.);
1057 Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.);
1058 Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.);
1059 Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.);
1060 Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.);
1061 Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.);
1062 Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.);
1063 Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.);
1065 static TVectorD parY0(3);
1066 static TMatrixD matY0(3, 3);
1067 static TVectorD parZ0(3);
1068 static TMatrixD matZ0(3, 3);
1069 fitY0.GetParameters(parY0);
1070 fitY0.GetCovarianceMatrix(matY0);
1071 fitZ0.GetParameters(parZ0);
1072 fitZ0.GetCovarianceMatrix(matZ0);
1074 static TVectorD parY2(5);
1075 static TMatrixD matY2(5,5);
1076 static TVectorD parZ2(5);
1077 static TMatrixD matZ2(5,5);
1078 fitY2.GetParameters(parY2);
1079 fitY2.GetCovarianceMatrix(matY2);
1080 fitZ2.GetParameters(parZ2);
1081 fitZ2.GetCovarianceMatrix(matZ2);
1083 static TVectorD parY2Q(5);
1084 static TMatrixD matY2Q(5,5);
1085 static TVectorD parZ2Q(5);
1086 static TMatrixD matZ2Q(5,5);
1087 fitY2Q.GetParameters(parY2Q);
1088 fitY2Q.GetCovarianceMatrix(matY2Q);
1089 fitZ2Q.GetParameters(parZ2Q);
1090 fitZ2Q.GetCovarianceMatrix(matZ2Q);
1091 static TVectorD parY2S(5);
1092 static TMatrixD matY2S(5,5);
1093 static TVectorD parZ2S(5);
1094 static TMatrixD matZ2S(5,5);
1095 fitY2S.GetParameters(parY2S);
1096 fitY2S.GetCovarianceMatrix(matY2S);
1097 fitZ2S.GetParameters(parZ2S);
1098 fitZ2S.GetCovarianceMatrix(matZ2S);
1099 Float_t sigmaY0 = TMath::Sqrt(matY0(0,0));
1100 Float_t sigmaZ0 = TMath::Sqrt(matZ0(0,0));
1101 Float_t sigmaDY0 = TMath::Sqrt(matY0(1,1));
1102 Float_t sigmaDZ0 = TMath::Sqrt(matZ0(1,1));
1103 Float_t sigmaY2 = TMath::Sqrt(matY2(1,1));
1104 Float_t sigmaZ2 = TMath::Sqrt(matZ2(1,1));
1105 Float_t sigmaDY2 = TMath::Sqrt(matY2(3,3));
1106 Float_t sigmaDZ2 = TMath::Sqrt(matZ2(3,3));
1107 Float_t sigmaY2Q = TMath::Sqrt(matY2Q(1,1));
1108 Float_t sigmaZ2Q = TMath::Sqrt(matZ2Q(1,1));
1109 Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
1110 Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
1111 Float_t sigmaY2S = TMath::Sqrt(matY2S(1,1));
1112 Float_t sigmaZ2S = TMath::Sqrt(matZ2S(1,1));
1113 Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
1114 Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
1117 Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
1118 Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
1120 Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1121 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1122 Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1123 TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1124 Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1125 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1126 Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1127 TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1131 Float_t meanRMSY = 0;
1132 Float_t meanRMSZ = 0;
1134 for (Int_t idelta = -2; idelta <= 2; idelta++){
1135 // loop over neighbourhood
1136 if (idelta+irow < 0 || idelta+irow > 159) continue;
1137 // if (idelta+irow>159) continue;
1138 AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
1139 if (!cluster) continue;
1140 meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
1141 meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
1147 Float_t rmsY = TMath::Sqrt(cluster0->GetSigmaY2());
1148 Float_t rmsZ = TMath::Sqrt(cluster0->GetSigmaZ2());
1149 Float_t rmsYT = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1150 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1151 Float_t rmsZT = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1152 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1153 Float_t rmsYT0 = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1154 TMath::Abs(angley));
1155 Float_t rmsZT0 = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1156 TMath::Abs(anglez));
1157 Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1158 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1159 Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1160 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1161 Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1162 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1164 Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1165 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1169 TTreeSRedirector *cstream = GetDebugStreamer();
1171 (*cstream)<<"ResolCl"<< // valgrind 3 40,000 bytes in 1 blocks are possibly
1174 "CSigmaY0="<<csigmaY0<< // cluster errorY
1175 "CSigmaYQ="<<csigmaYQ<< // cluster errorY - q scaled
1176 "CSigmaYS="<<csigmaYS<< // cluster errorY - q scaled
1177 "CSigmaZ0="<<csigmaZ0<< //
1178 "CSigmaZQ="<<csigmaZQ<<
1179 "CSigmaZS="<<csigmaZS<<
1180 "shapeYF="<<rmsYFactor<<
1181 "shapeZF="<<rmsZFactor<<
1184 "rmsYM="<<meanRMSY<<
1185 "rmsZM="<<meanRMSZ<<
1190 "rmsYS="<<rmsYSigma<<
1191 "rmsZS="<<rmsZSigma<<
1192 "padSize="<<padSize<<
1196 "SigmaY0="<<sigmaY0<<
1197 "SigmaZ0="<<sigmaZ0<<
1203 (*cstream)<<"ResolTr"<<
1204 "padSize="<<padSize<<
1212 "chi2Y2Q="<<chi2Y2Q<<
1213 "chi2Z2Q="<<chi2Z2Q<<
1214 "chi2Y2S="<<chi2Y2S<<
1215 "chi2Z2S="<<chi2Z2S<<
1224 "SigmaY0="<<sigmaY0<<
1225 "SigmaZ0="<<sigmaZ0<<
1226 "SigmaDY0="<<sigmaDY0<<
1227 "SigmaDZ0="<<sigmaDZ0<<
1228 "SigmaY2="<<sigmaY2<<
1229 "SigmaZ2="<<sigmaZ2<<
1230 "SigmaDY2="<<sigmaDY2<<
1231 "SigmaDZ2="<<sigmaDZ2<<
1232 "SigmaY2Q="<<sigmaY2Q<<
1233 "SigmaZ2Q="<<sigmaZ2Q<<
1234 "SigmaDY2Q="<<sigmaDY2Q<<
1235 "SigmaDZ2Q="<<sigmaDZ2Q<<
1236 "SigmaY2S="<<sigmaY2S<<
1237 "SigmaZ2S="<<sigmaZ2S<<
1238 "SigmaDY2S="<<sigmaDY2S<<
1239 "SigmaDZ2S="<<sigmaDZ2S<<
1250 TH2D * AliTPCcalibTracks::MakeDiff(TH2D * hfit, TF2 * func){
1252 // creates a new histogram which contains the difference between
1253 // the histogram hfit and the function func
1255 TH2D * result = (TH2D*)hfit->Clone(); // valgrind 3 40,139 bytes in 11 blocks are still reachable
1256 result->SetTitle(Form("%s fit residuals",result->GetTitle()));
1257 result->SetName(Form("%s fit residuals",result->GetName()));
1258 TAxis *xaxis = hfit->GetXaxis();
1259 TAxis *yaxis = hfit->GetYaxis();
1261 for (Int_t biny = 0; biny <= yaxis->GetNbins(); biny++) {
1262 x[1] = yaxis->GetBinCenter(biny);
1263 for (Int_t binx = 0; binx <= xaxis->GetNbins(); binx++) {
1264 x[0] = xaxis->GetBinCenter(binx);
1265 Int_t bin = hfit->GetBin(binx, biny);
1266 Double_t val = hfit->GetBinContent(bin);
1267 // result->SetBinContent( bin, (val - func->Eval(x[0], x[1])) / func->Eval(x[0], x[1]) );
1268 result->SetBinContent( bin, (val / func->Eval(x[0], x[1])) - 1 );
1275 void AliTPCcalibTracks::SetStyle() const {
1277 // set style, can be called by all draw functions
1279 gROOT->SetStyle("Plain");
1280 gStyle->SetFillColor(10);
1281 gStyle->SetPadColor(10);
1282 gStyle->SetCanvasColor(10);
1283 gStyle->SetStatColor(10);
1284 gStyle->SetPalette(1,0);
1285 gStyle->SetNumberContours(60);
1289 void AliTPCcalibTracks::Draw(Option_t* opt){
1291 // draw-function of AliTPCcalibTracks
1292 // will draws some exemplaric pictures
1295 if (GetDebugLevel() > 6) Info("Draw", "Drawing an exemplaric picture.");
1299 TCanvas *c1 = new TCanvas();
1301 TVirtualPad *upperThird = c1->GetPad(1);
1302 TVirtualPad *middleThird = c1->GetPad(2);
1303 TVirtualPad *lowerThird = c1->GetPad(3);
1304 upperThird->Divide(2,0);
1305 TVirtualPad *upleft = upperThird->GetPad(1);
1306 TVirtualPad *upright = upperThird->GetPad(2);
1307 middleThird->Divide(2,0);
1308 TVirtualPad *middleLeft = middleThird->GetPad(1);
1309 TVirtualPad *middleRight = middleThird->GetPad(2);
1310 lowerThird->Divide(2,0);
1311 TVirtualPad *downLeft = lowerThird->GetPad(1);
1312 TVirtualPad *downRight = lowerThird->GetPad(2);
1316 min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1317 max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1318 fDeltaY->SetAxisRange(min, max);
1319 fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable
1323 max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1324 min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1325 fDeltaZ->SetAxisRange(min, max);
1326 fDeltaZ->Fit("gaus","q","",min, max);
1333 fRejectedTracksHisto->Draw(opt);
1334 TPaveText *pt = new TPaveText(0.6,0.6, 0.8,0.8, "NDC");
1335 TText *t1 = pt->AddText("1: kEdgeThetaCutNoise");
1336 TText *t2 = pt->AddText("2: kMinClusters");
1337 TText *t3 = pt->AddText("3: kMinRatio");
1338 TText *t4 = pt->AddText("4: kMax1pt");
1339 t1 = t1; t2 = t2; t3 = t3; t4 = t4; // avoid compiler warnings
1340 pt->SetToolTipText("Legend for failed cuts");
1344 fHclusterPerPadrowRaw->Draw(opt);
1347 fHclusterPerPadrow->Draw(opt);
1351 void AliTPCcalibTracks::MakeReport(Int_t stat, char* pathName){
1353 // all functions are called, that produce pictures
1354 // the histograms are written to the directory 'pathName'
1355 // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
1356 // 'stat' is also the number of minEntries for MakeResPlotsQTree
1359 if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
1360 MakeAmpPlots(stat, pathName);
1361 MakeDeltaPlots(pathName);
1362 FitResolutionNew(pathName);
1363 FitRMSNew(pathName);
1364 MakeChargeVsDriftLengthPlots(pathName);
1365 // MakeResPlotsQ(1, 1);
1366 MakeResPlotsQTree(stat, pathName);
1370 void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){
1372 // creates several plots:
1373 // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps
1374 // fArrayAmp.ps: one histogram per sector, the histogram shows the charge per cluster
1375 // fArrayAmpRow.ps: one histogram per sector, mean max. amplitude vs. pad row with landau fit
1376 // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1377 // Empty histograms (sectors without data) are not written to file
1378 // the ps-files are written to the directory 'pathName', that is created if it does not exist
1379 // 'stat': only histograms with more than 'stat' entries are written to file.
1383 gSystem->MakeDirectory(pathName);
1384 gSystem->ChangeDirectory(pathName);
1386 TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1388 // histograms with accumulated amplitude for all IROCs and OROCs
1389 TH1F *allAmpHisIROC = ((TH1F*)(fArrayAmp->At(0))->Clone());
1390 allAmpHisIROC->SetName("Amp all IROCs");
1391 allAmpHisIROC->SetTitle("Amp all IROCs");
1392 TH1F *allAmpHisOROC = ((TH1F*)(fArrayAmp->At(36))->Clone());
1393 allAmpHisOROC->SetName("Amp all OROCs");
1394 allAmpHisOROC->SetTitle("Amp all OROCs");
1396 ps = new TPostScript("fArrayAmp.ps", 112);
1397 if (GetDebugLevel() > -1) cout << "creating fArrayAmp.ps..." << endl;
1398 for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++){
1399 if ( ((TH1F*)fArrayAmp->At(i))->GetEntries() < stat ) continue;
1401 ((TH1F*)fArrayAmp->At(i))->Draw();
1402 c1->Update(); // valgrind 3
1403 if (i > 0 && i < 36) {
1404 allAmpHisIROC->Add(((TH1F*)fArrayAmp->At(i)));
1405 allAmpHisOROC->Add(((TH1F*)fArrayAmp->At(i+36)));
1409 allAmpHisIROC->Draw();
1410 c1->Update(); // valgrind
1412 allAmpHisOROC->Draw();
1420 ps = new TPostScript("fArrayAmpRow.ps", 112);
1421 if (GetDebugLevel() > -1) cout << "creating fArrayAmpRow.ps..." << endl;
1422 for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++){
1423 his = (TH1F*)fArrayAmpRow->At(i);
1424 if (his->GetEntries() < stat) continue;
1426 min = TMath::Max( his->GetBinCenter(his->GetMaximumBin() )-100., 0.);
1427 max = his->GetBinCenter(5*his->GetMaximumBin()) + 100;
1428 his->SetAxisRange(min, max);
1429 his->Fit("pol3", "q", "", min, max);
1430 // his->Draw("error"); // don't use this line when you don't want to have empty pages in the ps-file
1436 gSystem->ChangeDirectory("..");
1440 void AliTPCcalibTracks::MakeDeltaPlots(char* pathName){
1442 // creates several plots:
1443 // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1444 // the ps-files are written to the directory 'pathName', that is created if it does not exist
1448 gSystem->MakeDirectory(pathName);
1449 gSystem->ChangeDirectory(pathName);
1451 TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1456 ps = new TPostScript("DeltaYZ.ps", 112);
1457 if (GetDebugLevel() > -1) cout << "creating DeltaYZ.ps..." << endl;
1458 min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1459 max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1460 fDeltaY->SetAxisRange(min, max);
1462 fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable
1465 max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1466 min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1467 fDeltaZ->SetAxisRange(min, max);
1468 fDeltaZ->Fit("gaus","q","",min, max);
1473 gSystem->ChangeDirectory("..");
1477 void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(char* pathName){
1479 // creates charge vs. driftlength plots, one TProfile for each ROC
1480 // is not correct like this, should be one TProfile for each sector and padsize
1484 gSystem->MakeDirectory(pathName);
1485 gSystem->ChangeDirectory(pathName);
1487 TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1489 ps = new TPostScript("chargeVsDriftlengthOld.ps", 112);
1490 if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlength.ps..." << endl;
1491 TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone());
1492 TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone());
1493 chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC");
1494 chargeVsDriftlengthAllIROCs->SetTitle("charge vs. driftlength, all IROCs");
1495 chargeVsDriftlengthAllOROCs->SetName("allAmpHisOROC");
1496 chargeVsDriftlengthAllOROCs->SetTitle("charge vs. driftlength, all OROCs");
1498 for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) {
1499 ((TProfile*)fArrayChargeVsDriftlength->At(i))->Draw();
1501 if (i > 0 && i < 36) {
1502 chargeVsDriftlengthAllIROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i)));
1503 chargeVsDriftlengthAllOROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i+36)));
1507 chargeVsDriftlengthAllIROCs->Draw();
1508 c1->Update(); // valgrind
1510 chargeVsDriftlengthAllOROCs->Draw();
1515 gSystem->ChangeDirectory("..");
1519 void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(char* pathName){
1521 // creates charge vs. driftlength plots, one TProfile for each ROC
1522 // under development....
1526 gSystem->MakeDirectory(pathName);
1527 gSystem->ChangeDirectory(pathName);
1529 TCanvas* c1 = new TCanvas("c1", "c1", 700,(Int_t)(TMath::Sqrt(2)*700)); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1530 // TCanvas c1("c1", "c1", 500,(sqrt(2)*500))
1533 ps = new TPostScript("chargeVsDriftlength.ps", 111);
1534 if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl;
1536 TProfile *chargeVsDriftlengthAllShortPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone());
1537 TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone());
1538 TProfile *chargeVsDriftlengthAllLongPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)->Clone());
1539 chargeVsDriftlengthAllShortPads->SetName("allAmpHisShortPads");
1540 chargeVsDriftlengthAllShortPads->SetTitle("charge vs. driftlength, all sectors, short pads");
1541 chargeVsDriftlengthAllMediumPads->SetName("allAmpHisMediumPads");
1542 chargeVsDriftlengthAllMediumPads->SetTitle("charge vs. driftlength, all sectors, medium pads");
1543 chargeVsDriftlengthAllLongPads->SetName("allAmpHisLongPads");
1544 chargeVsDriftlengthAllLongPads->SetTitle("charge vs. driftlength, all sectors, long pads");
1546 for (Int_t i = 0; i < 36; i++) {
1547 c1->cd(1)->SetGridx();
1548 c1->cd(1)->SetGridy();
1549 ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,0))->Draw();
1550 c1->cd(2)->SetGridx();
1551 c1->cd(2)->SetGridy();
1552 ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,1))->Draw();
1553 c1->cd(3)->SetGridx();
1554 c1->cd(3)->SetGridy();
1555 ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,2))->Draw();
1557 chargeVsDriftlengthAllShortPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0));
1558 chargeVsDriftlengthAllMediumPads->Add((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1));
1559 chargeVsDriftlengthAllLongPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2));
1562 c1->cd(1)->SetGridx();
1563 c1->cd(1)->SetGridy();
1564 chargeVsDriftlengthAllShortPads->Draw();
1565 c1->cd(2)->SetGridx();
1566 c1->cd(2)->SetGridy();
1567 chargeVsDriftlengthAllMediumPads->Draw();
1568 c1->cd(3)->SetGridx();
1569 c1->cd(3)->SetGridy();
1570 chargeVsDriftlengthAllLongPads->Draw();
1571 c1->Update(); // valgrind
1576 gSystem->ChangeDirectory("..");
1581 void AliTPCcalibTracks::FitResolutionNew(char* pathName){
1583 // calculates different resulution fits in Y and Z direction
1584 // the histograms are written to 'ResolutionYZ.ps'
1585 // writes calculated resolution to 'resol.txt'
1586 // all files are stored in the directory pathName
1590 gSystem->MakeDirectory(pathName);
1591 gSystem->ChangeDirectory(pathName);
1595 if (GetDebugLevel() > -1) cout << "creating ResolutionYZ.ps..." << endl;
1596 TPostScript *ps = new TPostScript("ResolutionYZ.ps", 112);
1597 TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1598 fres->SetParameter(0,0.02);
1599 fres->SetParameter(1,0.0054);
1600 fres->SetParameter(2,0.13);
1602 TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory
1604 // create histogramw for Y-resolution
1605 TH3F * hisResY0 = (TH3F*)fResolY->At(0);
1606 hisResY0->FitSlicesZ();
1607 TH2D * hisResY02 = (TH2D*)gDirectory->Get("Resol Y0_2");
1608 TH3F * hisResY1 = (TH3F*)fResolY->At(1);
1609 hisResY1->FitSlicesZ();
1610 TH2D * hisResY12 = (TH2D*)gDirectory->Get("Resol Y1_2");
1611 TH3F * hisResY2 = (TH3F*)fResolY->At(2);
1612 hisResY2->FitSlicesZ();
1613 TH2D * hisResY22 = (TH2D*)gDirectory->Get("Resol Y2_2");
1617 hisResY02->Fit(fres, "q"); // valgrind 132,072 bytes in 6 blocks are indirectly lost
1618 hisResY02->Draw("surf1");
1620 MakeDiff(hisResY02,fres)->Draw("surf1");
1622 // c.SaveAs("ResolutionYPad0.eps");
1625 hisResY12->Fit(fres, "q");
1626 hisResY12->Draw("surf1");
1628 MakeDiff(hisResY12,fres)->Draw("surf1");
1630 // c.SaveAs("ResolutionYPad1.eps");
1633 hisResY22->Fit(fres, "q");
1634 hisResY22->Draw("surf1");
1636 MakeDiff(hisResY22,fres)->Draw("surf1");
1638 // c.SaveAs("ResolutionYPad2.eps");
1640 // create histogramw for Z-resolution
1641 TH3F * hisResZ0 = (TH3F*)fResolZ->At(0);
1642 hisResZ0->FitSlicesZ();
1643 TH2D * hisResZ02 = (TH2D*)gDirectory->Get("Resol Z0_2");
1644 TH3F * hisResZ1 = (TH3F*)fResolZ->At(1);
1645 hisResZ1->FitSlicesZ();
1646 TH2D * hisResZ12 = (TH2D*)gDirectory->Get("Resol Z1_2");
1647 TH3F * hisResZ2 = (TH3F*)fResolZ->At(2);
1648 hisResZ2->FitSlicesZ();
1649 TH2D * hisResZ22 = (TH2D*)gDirectory->Get("Resol Z2_2");
1653 hisResZ02->Fit(fres, "q");
1654 hisResZ02->Draw("surf1");
1656 MakeDiff(hisResZ02,fres)->Draw("surf1");
1658 // c.SaveAs("ResolutionZPad0.eps");
1661 hisResZ12->Fit(fres, "q");
1662 hisResZ12->Draw("surf1");
1664 MakeDiff(hisResZ12,fres)->Draw("surf1");
1666 // c.SaveAs("ResolutionZPad1.eps");
1669 hisResZ22->Fit(fres, "q");
1670 hisResZ22->Draw("surf1");
1672 MakeDiff(hisResZ22,fres)->Draw("surf1");
1674 // c.SaveAs("ResolutionZPad2.eps");
1678 // write calculated resoltuions to 'resol.txt'
1679 ofstream fresol("resol.txt");
1680 fresol<<"Pad 0.75 cm"<<"\n";
1681 hisResY02->Fit(fres, "q"); // valgrind
1682 fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1683 hisResZ02->Fit(fres, "q");
1684 fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1686 fresol<<"Pad 1.00 cm"<<1<<"\n";
1687 hisResY12->Fit(fres, "q"); // valgrind
1688 fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1689 hisResZ12->Fit(fres, "q");
1690 fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1692 fresol<<"Pad 1.50 cm"<<0<<"\n";
1693 hisResY22->Fit(fres, "q");
1694 fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1695 hisResZ22->Fit(fres, "q");
1696 fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1698 TH1::AddDirectory(kFALSE);
1699 gSystem->ChangeDirectory("..");
1704 void AliTPCcalibTracks::FitRMSNew(char* pathName){
1706 // calculates different resulution-rms fits in Y and Z direction
1707 // the histograms are written to 'RMS_YZ.ps'
1708 // writes calculated resolution rms to 'rms.txt'
1709 // all files are stored in the directory pathName
1713 gSystem->MakeDirectory(pathName);
1714 gSystem->ChangeDirectory(pathName);
1716 TCanvas c; // valgrind 3 42,120 bytes in 405 blocks are still reachable 23,816 bytes in 229 blocks are still reachable
1718 if (GetDebugLevel() > -1) cout << "creating RMS_YZ.ps..." << endl;
1719 TPostScript *ps = new TPostScript("RMS_YZ.ps", 112);
1720 TF2 *frms = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1721 frms->SetParameter(0,0.02);
1722 frms->SetParameter(1,0.0054);
1723 frms->SetParameter(2,0.13);
1725 TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory
1727 // create histogramw for Y-RMS
1728 TH3F * hisResY0 = (TH3F*)fRMSY->At(0);
1729 hisResY0->FitSlicesZ();
1730 TH2D * hisResY02 = (TH2D*)gDirectory->Get("RMS Y0_1");
1731 TH3F * hisResY1 = (TH3F*)fRMSY->At(1);
1732 hisResY1->FitSlicesZ();
1733 TH2D * hisResY12 = (TH2D*)gDirectory->Get("RMS Y1_1");
1734 TH3F * hisResY2 = (TH3F*)fRMSY->At(2);
1735 hisResY2->FitSlicesZ();
1736 TH2D * hisResY22 = (TH2D*)gDirectory->Get("RMS Y2_1");
1740 hisResY02->Fit(frms, "qn0");
1741 hisResY02->Draw("surf1");
1743 MakeDiff(hisResY02,frms)->Draw("surf1");
1745 // c.SaveAs("RMSYPad0.eps");
1748 hisResY12->Fit(frms, "qn0"); // valgrind several blocks possibly lost
1749 hisResY12->Draw("surf1");
1751 MakeDiff(hisResY12,frms)->Draw("surf1");
1753 // c.SaveAs("RMSYPad1.eps");
1756 hisResY22->Fit(frms, "qn0");
1757 hisResY22->Draw("surf1");
1759 MakeDiff(hisResY22,frms)->Draw("surf1");
1761 // c.SaveAs("RMSYPad2.eps");
1763 // create histogramw for Z-RMS
1764 TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0);
1765 hisResZ0->FitSlicesZ();
1766 TH2D * hisResZ02 = (TH2D*)gDirectory->Get("RMS Z0_1");
1767 TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1);
1768 hisResZ1->FitSlicesZ();
1769 TH2D * hisResZ12 = (TH2D*)gDirectory->Get("RMS Z1_1");
1770 TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2);
1771 hisResZ2->FitSlicesZ();
1772 TH2D * hisResZ22 = (TH2D*)gDirectory->Get("RMS Z2_1");
1776 hisResZ02->Fit(frms, "qn0"); // valgrind
1777 hisResZ02->Draw("surf1");
1779 MakeDiff(hisResZ02,frms)->Draw("surf1");
1781 // c.SaveAs("RMSZPad0.eps");
1784 hisResZ12->Fit(frms, "qn0");
1785 hisResZ12->Draw("surf1");
1787 MakeDiff(hisResZ12,frms)->Draw("surf1");
1789 // c.SaveAs("RMSZPad1.eps");
1792 hisResZ22->Fit(frms, "qn0"); // valgrind 1 block possibly lost
1793 hisResZ22->Draw("surf1");
1795 MakeDiff(hisResZ22,frms)->Draw("surf1");
1797 // c.SaveAs("RMSZPad2.eps");
1799 // write calculated resoltuion rms to 'rms.txt'
1800 ofstream filerms("rms.txt");
1801 filerms<<"Pad 0.75 cm"<<"\n";
1802 hisResY02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost
1803 filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1804 hisResZ02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost
1805 filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1807 filerms<<"Pad 1.00 cm"<<1<<"\n";
1808 hisResY12->Fit(frms, "qn0"); // valgrind 3,256 bytes in 22 blocks are indirectly lost
1809 filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1810 hisResZ12->Fit(frms, "qn0"); // valgrind 66,036 bytes in 3 blocks are still reachable
1811 filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1813 filerms<<"Pad 1.50 cm"<<0<<"\n";
1814 hisResY22->Fit(frms, "qn0"); // valgrind 40,139 bytes in 11 blocks are still reachable 330,180 bytes in 15 blocks are possibly lost
1815 filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1816 hisResZ22->Fit(frms, "qn0");
1817 filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1819 TH1::AddDirectory(kFALSE);
1820 gSystem->ChangeDirectory("..");
1827 void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){
1829 // Make tree with resolution parameters
1830 // the result is written to 'resol.root' in directory 'pathname'
1831 // file information are available in fileInfo
1832 // available variables in the tree 'Resol':
1833 // Entries: number of entries for this resolution point
1834 // nbins: number of bins that were accumulated
1835 // Dim: direction, Dim==0: y-direction, Dim==1: z-direction
1836 // Pad: padSize; short, medium and long
1837 // Length: pad length, 0.75, 1, 1.5
1838 // QMean: mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1839 // Zc: center of middle bin in drift direction
1840 // Zm: mean dirftlength for accumulated Delta-Histograms
1841 // Zs: width of driftlength bin
1842 // AngleC: center of middle bin in Angle-Direction
1843 // AngleM: mean angle for accumulated Delta-Histograms
1844 // AngleS: width of Angle-bin
1845 // Resol: sigma for gaus fit through Delta-Histograms
1846 // Sigma: error of sigma for gaus fit through Delta Histograms
1847 // MeanR: mean of the Delta-Histogram
1848 // SigmaR: rms of the Delta-Histogram
1849 // RMSm: mean of the gaus fit through RMS-Histogram
1850 // RMS: sigma of the gaus fit through RMS-Histogram
1851 // RMSe0: error of mean of gaus fit in RMS-Histogram
1852 // RMSe1: error of sigma of gaus fit in RMS-Histogram
1855 if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
1856 if (GetDebugLevel() > -1) cout << " relax, the calculation will take a while..." << endl;
1858 gSystem->MakeDirectory(pathName);
1859 gSystem->ChangeDirectory(pathName);
1860 TString kFileName = "resol.root";
1861 TTreeSRedirector fTreeResol(kFileName.Data());
1863 TH3F *resArray[2][3][11];
1864 TH3F *rmsArray[2][3][11];
1866 // load histograms from fArrayQDY and fArrayQDZ
1867 // into resArray and rmsArray
1868 // that is all we need here
1869 for (Int_t idim = 0; idim < 2; idim++){
1870 for (Int_t ipad = 0; ipad < 3; ipad++){
1871 for (Int_t iq = 0; iq <= 10; iq++){
1872 rmsArray[idim][ipad][iq]=0;
1873 resArray[idim][ipad][iq]=0;
1874 Int_t bin = GetBin(iq,ipad);
1876 if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
1877 if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
1878 if (!hresl) continue;
1879 resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
1880 resArray[idim][ipad][iq]->SetDirectory(0);
1881 TH3F * hreslRMS = 0;
1882 if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
1883 if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
1884 if (!hreslRMS) continue;
1885 rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
1886 rmsArray[idim][ipad][iq]->SetDirectory(0);
1890 if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
1892 //--------------------------------------------------------------------------------------------
1896 Double_t zMean, angleMean, zCenter, angleCenter;
1897 Double_t zSigma, angleSigma;
1898 TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
1899 TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
1900 TF1 *fitFunction = new TF1("fitFunction", "gaus");
1901 Float_t entriesQ = 0;
1902 Int_t loopCounter = 1;
1904 for (Int_t idim = 0; idim < 2; idim++){
1905 // Loop y-z corrdinate
1906 for (Int_t ipad = 0; ipad < 3; ipad++){
1908 for (Int_t iq = -1; iq < 10; iq++){
1910 if (GetDebugLevel() > -1)
1911 cout << "Loop-counter, this is loop " << loopCounter << " of 66, ("
1912 << (Int_t)((loopCounter)/66.*100) << "% done), "
1913 << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << " \r" << std::flush;
1922 // integrated spectra
1923 for (Int_t iql = 0; iql < 10; iql++){
1924 Int_t bin = GetBin(iql,ipad);
1925 TH3F *hresl = resArray[idim][ipad][iql];
1926 TH3F *hrmsl = rmsArray[idim][ipad][iql];
1927 if (!hresl) continue;
1928 if (!hrmsl) continue;
1929 entriesQ += hresl->GetEntries();
1930 qMean += hresl->GetEntries() * GetQ(bin);
1932 hres = (TH3F*)hresl->Clone();
1933 hrms = (TH3F*)hrmsl->Clone();
1941 qMean *= -1.; // integral mean charge
1944 // loop over neighboured Q-bins
1945 // accumulate entries from neighboured Q-bins
1946 for (Int_t iql = iq - 1; iql <= iq + 1; iql++){
1947 if (iql < 0) continue;
1948 Int_t bin = GetBin(iql,ipad);
1949 TH3F * hresl = resArray[idim][ipad][iql];
1950 TH3F * hrmsl = rmsArray[idim][ipad][iql];
1951 if (!hresl) continue;
1952 if (!hrmsl) continue;
1953 entriesQ += hresl->GetEntries();
1954 qMean += hresl->GetEntries() * GetQ(bin);
1956 hres = (TH3F*) hresl->Clone();
1957 hrms = (TH3F*) hrmsl->Clone();
1967 TAxis *xAxisDriftLength = hres->GetXaxis(); // driftlength / z - axis
1968 TAxis *yAxisAngle = hres->GetYaxis(); // angle axis
1969 TAxis *zAxisDelta = hres->GetZaxis(); // delta axis
1970 TAxis *zAxisRms = hrms->GetZaxis(); // rms axis
1972 // loop over all angle bins
1973 for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
1974 angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
1975 // loop over all driftlength bins
1976 for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
1977 zCenter = xAxisDriftLength->GetBinCenter(ibinxDL);
1978 zSigma = xAxisDriftLength->GetBinWidth(ibinxDL);
1979 angleSigma = yAxisAngle->GetBinWidth(ibinyAngle);
1980 zMean = zCenter; // changens, when more statistic is accumulated
1981 angleMean = angleCenter; // changens, when more statistic is accumulated
1983 // create 2 1D-Histograms, projectionRes and projectionRms
1984 // these histograms are delta histograms for given direction, padSize, chargeBin,
1985 // angleBin and driftLengthBin
1986 // later on they will be fitted with a gausian, its sigma is the resoltuion...
1987 sprintf(name,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
1988 // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1989 projectionRes->SetNameTitle(name, name);
1990 sprintf(name,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
1991 // TH1D * projectionRms = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1992 projectionRms->SetNameTitle(name, name);
1994 projectionRes->Reset();
1995 projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
1996 projectionRms->Reset();
1997 projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
1998 projectionRes->SetDirectory(0);
1999 projectionRms->SetDirectory(0);
2001 Double_t entries = 0;
2002 Int_t nbins = 0; // counts, how many bins were accumulated
2007 // fill projectionRes and projectionRms for given dim, ipad and iq,
2008 // as well as for given angleBin and driftlengthBin
2009 // if this gives not enough statistic, include neighbourhood
2010 // (angle and driftlength) successifely
2011 for (Int_t dbin = 0; dbin <= 8; dbin++){ // delta-bins around centered angleBin and driftlengthBin
2012 for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) { // delta-bins in angle direction
2013 for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
2014 if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue; // add each bin only one time !
2015 Int_t binx2 = ibinxDL + dbinx2; // position variable in x (driftlength) direction
2016 Int_t biny2 = ibinyAngle + dbiny2; // position variable in y (angle) direction
2017 if (binx2 < 1 || biny2 < 1) continue; // don't go out of the histogram!
2018 if (binx2 >= xAxisDriftLength->GetNbins()) continue; // don't go out of the histogram!
2019 if (biny2 >= yAxisAngle->GetNbins()) continue; // don't go out of the histogram!
2020 nbins++; // count the number of accumulated bins
2021 // Fill resolution histo
2022 for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
2023 // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3); // unused variable
2024 projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
2025 entries += hres->GetBinContent(binx2, biny2, ibin3);
2026 zMean += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
2027 angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
2030 for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
2031 projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
2034 if (entries > minEntries) break; // enough statistic accumulated
2036 if (entries > minEntries) break; // enough statistic accumulated
2038 if ( entries< minEntries) continue; // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree
2040 angleMean /= entries;
2042 if (entries > minEntries) {
2043 // when enough statistic is accumulated
2044 // fit Delta histograms with a gausian
2045 // of the gausian is the resolution (resol), its fit error is sigma
2046 // store also mean and RMS of the histogram
2047 Float_t xmin = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2048 Float_t xmax = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2050 // projectionRes->Fit("gaus", "q0", "", xmin, xmax);
2051 // Float_t resol = projectionRes->GetFunction("gaus")->GetParameter(2);
2052 // Float_t sigma = projectionRes->GetFunction("gaus")->GetParError(2);
2053 fitFunction->SetMaximum(xmax);
2054 fitFunction->SetMinimum(xmin);
2055 projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
2056 Float_t resol = fitFunction->GetParameter(2);
2057 Float_t sigma = fitFunction->GetParError(2);
2059 Float_t meanR = projectionRes->GetMean();
2060 Float_t sigmaR = projectionRes->GetRMS();
2061 // fit also RMS histograms with a gausian
2062 // store mean and sigma of the gausian in rmsMean and rmsSigma
2063 // store also the fit errors in errorRMS and errorSigma
2064 xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2065 xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2067 // projectionRms->Fit("gaus","q0","",xmin,xmax);
2068 // Float_t rmsMean = projectionRms->GetFunction("gaus")->GetParameter(1);
2069 // Float_t rmsSigma = projectionRms->GetFunction("gaus")->GetParameter(2);
2070 // Float_t errorRMS = projectionRms->GetFunction("gaus")->GetParError(1);
2071 // Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
2072 projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
2073 Float_t rmsMean = fitFunction->GetParameter(1);
2074 Float_t rmsSigma = fitFunction->GetParameter(2);
2075 Float_t errorRMS = fitFunction->GetParError(1);
2076 Float_t errorSigma = fitFunction->GetParError(2);
2078 Float_t length = 0.75;
2079 if (ipad == 1) length = 1;
2080 if (ipad == 2) length = 1.5;
2082 fTreeResol<<"Resol"<<
2083 "Entries="<<entries<< // number of entries for this resolution point
2084 "nbins="<<nbins<< // number of bins that were accumulated
2085 "Dim="<<idim<< // direction, Dim==0: y-direction, Dim==1: z-direction
2086 "Pad="<<ipad<< // padSize; short, medium and long
2087 "Length="<<length<< // pad length, 0.75, 1, 1.5
2088 "QMean="<<qMean<< // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
2089 "Zc="<<zCenter<< // center of middle bin in drift direction
2090 "Zm="<<zMean<< // mean dirftlength for accumulated Delta-Histograms
2091 "Zs="<<zSigma<< // width of driftlength bin
2092 "AngleC="<<angleCenter<< // center of middle bin in Angle-Direction
2093 "AngleM="<<angleMean<< // mean angle for accumulated Delta-Histograms
2094 "AngleS="<<angleSigma<< // width of Angle-bin
2095 "Resol="<<resol<< // sigma for gaus fit through Delta-Histograms
2096 "Sigma="<<sigma<< // error of sigma for gaus fit through Delta Histograms
2097 "MeanR="<<meanR<< // mean of the Delta-Histogram
2098 "SigmaR="<<sigmaR<< // rms of the Delta-Histogram
2099 "RMSm="<<rmsMean<< // mean of the gaus fit through RMS-Histogram
2100 "RMSs="<<rmsSigma<< // sigma of the gaus fit through RMS-Histogram
2101 "RMSe0="<<errorRMS<< // error of mean of gaus fit in RMS-Histogram
2102 "RMSe1="<<errorSigma<< // error of sigma of gaus fit in RMS-Histogram
2104 if (GetDebugLevel() > 5) {
2105 projectionRes->SetDirectory(fTreeResol.GetFile());
2106 projectionRes->Write(projectionRes->GetName());
2107 projectionRes->SetDirectory(0);
2108 projectionRms->SetDirectory(fTreeResol.GetFile());
2109 projectionRms->Write(projectionRms->GetName());
2110 projectionRes->SetDirectory(0);
2112 } // if (projectionRes->GetSum() > minEntries)
2113 } // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
2114 } // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
2119 delete projectionRes;
2120 delete projectionRms;
2122 // TFile resolFile(fTreeResol.GetFile());
2123 TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
2124 fileInfo.Write("fileInfo");
2125 // resolFile.Close();
2126 // fTreeResol.GetFile()->Close();
2127 if (GetDebugLevel() > -1) cout << endl;
2128 if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
2129 gSystem->ChangeDirectory("..");
2136 Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
2138 // function to merge several AliTPCcalibTracks objects after PROOF calculation
2139 // The object's histograms are merged via their merge functions
2140 // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
2143 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks::Merge(TCollection *collectionList) *****"<< endl;
2144 if (!collectionList) return 0;
2145 if (collectionList->IsEmpty()) return -1;
2147 if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
2148 if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
2149 collectionList->Print();
2151 // create a list for each data member
2152 TList* deltaYList = new TList;
2153 TList* deltaZList = new TList;
2154 TList* arrayAmpRowList = new TList;
2155 TList* rejectedTracksList = new TList;
2156 TList* hclusList = new TList;
2157 TList* clusterCutHistoList = new TList;
2158 TList* arrayAmpList = new TList;
2159 TList* arrayQDYList = new TList;
2160 TList* arrayQDZList = new TList;
2161 TList* arrayQRMSYList = new TList;
2162 TList* arrayQRMSZList = new TList;
2163 TList* arrayChargeVsDriftlengthList = new TList;
2164 TList* calPadRegionChargeVsDriftlengthList = new TList;
2165 TList* hclusterPerPadrowList = new TList;
2166 TList* hclusterPerPadrowRawList = new TList;
2167 TList* resolYList = new TList;
2168 TList* resolZList = new TList;
2169 TList* rMSYList = new TList;
2170 TList* rMSZList = new TList;
2172 // TList* nRowsList = new TList;
2173 // TList* nSectList = new TList;
2174 // TList* fileNoList = new TList;
2176 TIterator *listIterator = collectionList->MakeIterator();
2177 AliTPCcalibTracks *calibTracks = 0;
2178 if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;
2180 while ( (calibTracks = (AliTPCcalibTracks*)listIterator->Next()) ){
2181 // loop over all entries in the collectionList and get dataMembers into lists
2182 if (!calibTracks) continue;
2184 deltaYList->Add( calibTracks->GetfDeltaY() );
2185 deltaZList->Add( calibTracks->GetfDeltaZ() );
2186 arrayAmpRowList->Add(calibTracks->GetfArrayAmpRow());
2187 arrayAmpList->Add(calibTracks->GetfArrayAmp());
2188 arrayQDYList->Add(calibTracks->GetfArrayQDY());
2189 arrayQDZList->Add(calibTracks->GetfArrayQDZ());
2190 arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
2191 arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
2192 resolYList->Add(calibTracks->GetfResolY());
2193 resolZList->Add(calibTracks->GetfResolZ());
2194 rMSYList->Add(calibTracks->GetfRMSY());
2195 rMSZList->Add(calibTracks->GetfRMSZ());
2196 arrayChargeVsDriftlengthList->Add(calibTracks->GetfArrayChargeVsDriftlength());
2197 calPadRegionChargeVsDriftlengthList->Add(calibTracks->GetCalPadRegionchargeVsDriftlength());
2198 hclusList->Add(calibTracks->GetfHclus());
2199 rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
2200 clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
2201 hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow());
2202 hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw());
2203 fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
2204 fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
2206 if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
2210 // merge data members
2211 if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl;
2212 fDeltaY->Merge(deltaYList);
2213 fDeltaZ->Merge(deltaZList);
2214 fHclus->Merge(hclusList);
2215 fClusterCutHisto->Merge(clusterCutHistoList);
2216 fRejectedTracksHisto->Merge(rejectedTracksList);
2217 fHclusterPerPadrow->Merge(hclusterPerPadrowList);
2218 fHclusterPerPadrowRaw->Merge(hclusterPerPadrowRawList);
2220 TObjArray* objarray = 0;
2222 TList* histList = 0;
2223 TIterator *objListIterator = 0;
2225 if (GetDebugLevel() > 0) cout << "merging fArrayAmpRows..." << endl;
2226 // merge fArrayAmpRows
2227 for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) { // loop over data member, i<72
2228 objListIterator = arrayAmpRowList->MakeIterator();
2229 histList = new TList;
2230 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2231 // loop over arrayAmpRowList, get TObjArray, get object at position i, cast it into TProfile
2232 if (!objarray) continue;
2233 hist = (TProfile*)(objarray->At(i));
2234 histList->Add(hist);
2236 ((TProfile*)(fArrayAmpRow->At(i)))->Merge(histList);
2238 delete objListIterator;
2241 if (GetDebugLevel() > 0) cout << "merging fArrayAmps..." << endl;
2243 for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) { // loop over data member, i<72
2244 TIterator *objListIterator = arrayAmpList->MakeIterator();
2245 histList = new TList;
2246 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2247 // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F
2248 if (!objarray) continue;
2249 hist = (TH1F*)(objarray->At(i));
2250 histList->Add(hist);
2252 ((TH1F*)(fArrayAmp->At(i)))->Merge(histList);
2254 delete objListIterator;
2257 if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
2259 for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
2260 objListIterator = arrayQDYList->MakeIterator();
2261 histList = new TList;
2262 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2263 // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
2264 if (!objarray) continue;
2265 hist = (TH3F*)(objarray->At(i));
2266 histList->Add(hist);
2268 ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
2270 delete objListIterator;
2273 if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
2275 for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2276 objListIterator = arrayQDZList->MakeIterator();
2277 histList = new TList;
2278 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2279 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2280 if (!objarray) continue;
2281 hist = (TH3F*)(objarray->At(i));
2282 histList->Add(hist);
2284 ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
2286 delete objListIterator;
2289 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
2290 // merge fArrayQRMSY
2291 for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
2292 objListIterator = arrayQRMSYList->MakeIterator();
2293 histList = new TList;
2294 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2295 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2296 if (!objarray) continue;
2297 hist = (TH3F*)(objarray->At(i));
2298 histList->Add(hist);
2300 ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
2302 delete objListIterator;
2305 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
2306 // merge fArrayQRMSZ
2307 for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2308 objListIterator = arrayQRMSZList->MakeIterator();
2309 histList = new TList;
2310 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2311 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2312 if (!objarray) continue;
2313 hist = (TH3F*)(objarray->At(i));
2314 histList->Add(hist);
2316 ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
2318 delete objListIterator;
2321 if (GetDebugLevel() > 0) cout << "merging fArrayChargeVsDriftlength..." << endl;
2322 // merge fArrayChargeVsDriftlength
2323 for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { // loop over data member, i < 300
2324 objListIterator = arrayChargeVsDriftlengthList->MakeIterator();
2325 histList = new TList;
2326 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2327 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TProfile
2328 if (!objarray) continue;
2329 hist = (TProfile*)(objarray->At(i));
2330 histList->Add(hist);
2332 ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Merge(histList);
2334 delete objListIterator;
2337 if (GetDebugLevel() > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl;
2338 // merge fcalPadRegionChargeVsDriftlength
2339 AliTPCCalPadRegion *cpr = 0x0;
2342 TIterator *regionIterator = fcalPadRegionChargeVsDriftlength->MakeIterator();
2343 while (hist = (TProfile*)regionIterator->Next()) {
2344 // loop over all calPadRegion's in destination calibTracks object
2345 objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2346 while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) {
2353 for (UInt_t isec = 0; isec < 36; isec++) {
2354 for (UInt_t padSize = 0; padSize < 3; padSize++){
2355 objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2356 histList = new TList;
2357 while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) {
2358 // loop over calPadRegionChargeVsDriftlengthList, get AliTPCCalPadRegion, get object
2360 hist = (TProfile*)cpr->GetObject(isec, padSize);
2361 histList->Add(hist);
2363 ((TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isec, padSize)))->Merge(histList);
2365 delete objListIterator;
2372 if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
2374 for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
2375 objListIterator = resolYList->MakeIterator();
2376 histList = new TList;
2377 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2378 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2379 if (!objarray) continue;
2380 hist = (TH3F*)(objarray->At(i));
2381 histList->Add(hist);
2383 ((TH3F*)(fResolY->At(i)))->Merge(histList);
2385 delete objListIterator;
2389 for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2390 objListIterator = resolZList->MakeIterator();
2391 histList = new TList;
2392 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2393 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2394 if (!objarray) continue;
2395 hist = (TH3F*)(objarray->At(i));
2396 histList->Add(hist);
2398 ((TH3F*)(fResolZ->At(i)))->Merge(histList);
2400 delete objListIterator;
2404 for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
2405 objListIterator = rMSYList->MakeIterator();
2406 histList = new TList;
2407 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2408 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2409 if (!objarray) continue;
2410 hist = (TH3F*)(objarray->At(i));
2411 histList->Add(hist);
2413 ((TH3F*)(fRMSY->At(i)))->Merge(histList);
2415 delete objListIterator;
2419 for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2420 objListIterator = rMSZList->MakeIterator();
2421 histList = new TList;
2422 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2423 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2424 if (!objarray) continue;
2425 hist = (TH3F*)(objarray->At(i));
2426 histList->Add(hist);
2428 ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
2430 delete objListIterator;
2435 delete arrayAmpRowList;
2436 delete arrayAmpList;
2437 delete arrayQDYList;
2438 delete arrayQDZList;
2439 delete arrayQRMSYList;
2440 delete arrayQRMSZList;
2445 delete listIterator;
2447 if (GetDebugLevel() > 0) cout << "merging done!" << endl;
2453 AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks){
2455 // function to test AliTPCcalibTrack::Merge:
2456 // in the file 'f' is a AliTPCcalibTrack object with name "calibTracks"
2457 // this object is appended 'nCalTracks' times to a TList
2458 // A new AliTPCcalibTrack object is created which merges the list
2459 // this object is returned
2462 // .L AliTPCcalibTracks.cxx+g
2464 TFile f("Output.root");
2465 AliTPCcalibTracks* calTracks = (AliTPCcalibTracks*)f.Get("calibTracks");
2467 TFile clusterParamFile("/u/lbozyk/calibration/workdir/calibTracks/TPCClusterParam.root");
2468 AliTPCClusterParam *clusterParam = (AliTPCClusterParam *) clusterParamFile.Get("Param");
2469 clusterParamFile.Close();
2471 AliTPCcalibTracks::TestMerge(calTracks, clusterParam);
2473 TList *list = new TList();
2474 if (ct == 0 || clusterParam == 0) return 0;
2475 cout << "making list with " << nCalTracks << " AliTPCcalibTrack objects" << endl;
2476 for (Int_t i = 0; i < nCalTracks; i++) {
2477 if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl;
2478 list->Add(new AliTPCcalibTracks(*ct));
2481 // only for check at the end
2482 AliTPCcalibTracks* cal1 = new AliTPCcalibTracks(*ct);
2483 Double_t cal1Entries = ((TH1F*)cal1->GetfArrayAmpRow()->At(5))->GetEntries();
2484 // Double_t cal1Entries = 5; //((TH1F*)ct->GetfArrayAmpRow()->At(5))->GetEntries();
2486 cout << "The list contains " << list->GetEntries() << " entries. " << endl;
2489 AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018);
2490 AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts, 5);
2493 cout << "cal->GetfArrayAmpRow()->At(5)->Print():" << endl;
2494 cal->GetfArrayAmpRow()->At(5)->Print();
2495 Double_t calEntries = ((TH1F*)cal->GetfArrayAmpRow()->At(5))->GetEntries();
2497 cout << "cal1->GetfArrayAmpRow()->At(5))->GetEntries() = " << cal1Entries << endl;
2498 cout << " cal->GetfArrayAmpRow()->At(5))->GetEntries() = " << calEntries << endl;
2499 printf("That means there were %f / %f = %f AliTPCcalibTracks-Objects merged. \n",
2500 calEntries, cal1Entries, ((Double_t)calEntries/cal1Entries));
2507 void AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints, Int_t verbose){
2509 // Make position corrections
2510 // for the moment Only using debug streamer
2511 // chainres - debug tree
2512 // param - parameters to be updated
2513 // maxPoints - maximal number of points using for fit
2514 // verbose - print info flag
2516 // Current implementation - only using debug streamers
2521 Int_t maxPoints=100000;
2526 gSystem->Load("libANALYSIS");
2527 gSystem->Load("libSTAT");
2528 gSystem->Load("libTPCcalib");
2531 //1. Load Parameters to be modified:
2534 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
2535 AliCDBManager::Instance()->SetRun(0)
2536 AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
2538 //2. Load chain from debug streamers
2541 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2542 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2543 AliXRDPROOFtoolkit tool;
2544 TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
2546 //3. Do fits and store results
2548 AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ;
2549 TFile f("paramout.root","recreate");
2550 param->Write("clusterParam");
2553 TFile f2("paramout.root");
2554 AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
2555 param2->SetInstance(param2);
2556 chainres->Draw("fitZ0:AliTPCClusterParam::SPosCorrection(1,0,Cl.fPad,Cl.fTimeBin,Cl.fZ,Cl.fSigmaY2,Cl.fSigmaZ2,Cl.fMax)","Cl.fDetector<36","",10000); // should be line
2561 TStatToolkit toolkit;
2563 TVectorD fitParamY0;
2564 TVectorD fitParamY1;
2565 TVectorD fitParamZ0;
2566 TVectorD fitParamZ1;
2570 chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
2571 chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
2572 chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
2573 chainres->SetAlias("st","(sin(dt)-dt)");
2575 chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
2576 chainres->SetAlias("dq","sqrt(15./(5+Cl.fMax))");
2577 chainres->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
2578 chainres->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
2583 TString fstringY="";
2585 fstringY+="(dp)++"; //1
2586 fstringY+="(dp)*di++"; //2
2587 fstringY+="(dp)*dq++"; //3
2588 fstringY+="(dp)*sy++"; //4
2590 fstringY+="(sp)++"; //5
2591 fstringY+="(sp)*di++"; //6
2592 fstringY+="(sp)*dq++"; //7
2593 fstringY+="(sp)*sy++"; //8
2596 TString fstringZ="";
2597 fstringZ+="(dt)++"; //1
2598 fstringZ+="(dt)*di++"; //2
2599 fstringZ+="(dt)*dq++"; //3
2600 fstringZ+="(dt)*sz++"; //4
2602 fstringZ+="(st)++"; //5
2603 fstringZ+="(st)*di++"; //6
2604 fstringZ+="(st)*dq++"; //7
2605 fstringZ+="(st)*sz++"; //8
2609 TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
2610 printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2611 param->fPosZcor[0] = (TVectorD*) fitParamZ0.Clone();
2613 TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
2615 printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2616 param->fPosZcor[1] = (TVectorD*) fitParamZ1.Clone();
2617 param->fPosZcor[2] = (TVectorD*) fitParamZ1.Clone();
2621 TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
2622 printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2623 param->fPosYcor[0] = (TVectorD*) fitParamY0.Clone();
2626 TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
2628 printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2629 param->fPosYcor[1] = (TVectorD*) fitParamY1.Clone();
2630 param->fPosYcor[2] = (TVectorD*) fitParamY1.Clone();
2634 chainres->SetAlias("fitZ0",strZ0->Data());
2635 chainres->SetAlias("fitZ1",strZ1->Data());
2636 chainres->SetAlias("fitY0",strY0->Data());
2637 chainres->SetAlias("fitY1",strY1->Data());
2638 // chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000);
2642 gSystem->Load("libSTAT.so");
2643 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2644 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2645 AliXRDPROOFtoolkit tool;
2646 TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
2654 TStatToolkit toolkit;
2656 TVectorD fitParamY0;
2657 TVectorD fitParamY1;
2658 TVectorD fitParamZ0;
2659 TVectorD fitParamZ1;
2663 chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
2664 chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
2665 chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
2666 chainres->SetAlias("st","(sin(dt)-dt)");
2668 chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
2669 chainres->SetAlias("dq","sqrt(15./(5+Cl.fMax))");
2670 chainres->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
2671 chainres->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
2677 TString fstringY="";
2679 fstringY+="(dp)++"; //1
2680 fstringY+="(dp)*di++"; //2
2681 fstringY+="(dp)*dq++"; //3
2682 fstringY+="(dp)*sy++"; //4
2684 fstringY+="(sp)++"; //5
2685 fstringY+="(sp)*di++"; //6
2686 fstringY+="(sp)*dq++"; //7
2687 fstringY+="(sp)*sy++"; //8
2690 TString fstringZ="";
2691 fstringZ+="(dt)++"; //1
2692 fstringZ+="(dt)*di++"; //2
2693 fstringZ+="(dt)*dq++"; //3
2694 fstringZ+="(dt)*sz++"; //4
2696 fstringZ+="(st)++"; //5
2697 fstringZ+="(st)*di++"; //6
2698 fstringZ+="(st)*dq++"; //7
2699 fstringZ+="(st)*sz++"; //8
2702 TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,100000);
2703 TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,100000);
2705 TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,100000);
2706 TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,100000);
2708 chainres->SetAlias("fitZ0",strZ0.Data());
2709 chainres->SetAlias("fitZ1",strZ1.Data());
2710 chainres->SetAlias("fitY0",strY0.Data());
2711 chainres->SetAlias("fitY1",strY1.Data());
2712 chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000)