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"
119 #include "THnSparse.h"
123 ClassImp(AliTPCcalibTracks)
126 AliTPCcalibTracks::AliTPCcalibTracks():
130 fHisDeltaY(0), // THnSparse - delta Y
131 fHisDeltaZ(0), // THnSparse - delta Z
132 fHisRMSY(0), // THnSparse - rms Y
133 fHisRMSZ(0), // THnSparse - rms Z
134 fHisQmax(0), // THnSparse - qmax
135 fHisQtot(0), // THnSparse - qtot
142 fArrayChargeVsDriftlength(0),
143 fcalPadRegionChargeVsDriftlength(0),
152 fRejectedTracksHisto(0),
153 fHclusterPerPadrow(0),
154 fHclusterPerPadrowRaw(0),
156 fCalPadClusterPerPad(0),
157 fCalPadClusterPerPadRaw(0)
160 // AliTPCcalibTracks default constructor
162 if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' default constructor called" << endl;
167 AliTPCcalibTracks::AliTPCcalibTracks(const AliTPCcalibTracks& calibTracks):
168 AliTPCcalibBase(calibTracks),
171 fHisDeltaY(0), // THnSparse - delta Y
172 fHisDeltaZ(0), // THnSparse - delta Z
173 fHisRMSY(0), // THnSparse - rms Y
174 fHisRMSZ(0), // THnSparse - rms Z
175 fHisQmax(0), // THnSparse - qmax
176 fHisQtot(0), // THnSparse - qtot
183 fArrayChargeVsDriftlength(0),
184 fcalPadRegionChargeVsDriftlength(0),
193 fRejectedTracksHisto(0),
194 fHclusterPerPadrow(0),
195 fHclusterPerPadrowRaw(0),
197 fCalPadClusterPerPad(0),
198 fCalPadClusterPerPadRaw(0)
201 // AliTPCcalibTracks copy constructor
203 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' copy constructor ***** " << endl;
205 Bool_t dirStatus = TH1::AddDirectoryStatus();
206 TH1::AddDirectory(kFALSE);
209 // backward compatibility: if the data member doesn't yet exist, it will not be merged
210 (calibTracks.fArrayAmpRow) ? length = calibTracks.fArrayAmpRow->GetEntriesFast() : length = -1;
211 fArrayAmpRow = new TObjArray(length);
212 fArrayAmp = new TObjArray(length);
213 for (Int_t i = 0; i < length; i++) {
214 fArrayAmpRow->AddAt( (TProfile*)calibTracks.fArrayAmpRow->At(i)->Clone(), i);
215 fArrayAmp->AddAt( ((TProfile*)calibTracks.fArrayAmp->At(i)->Clone()), i);
218 (calibTracks.fArrayQDY) ? length = calibTracks.fArrayQDY->GetEntriesFast() : length = -1;
219 fArrayQDY= new TObjArray(length);
220 fArrayQDZ= new TObjArray(length);
221 fArrayQRMSY= new TObjArray(length);
222 fArrayQRMSZ= new TObjArray(length);
223 for (Int_t i = 0; i < length; i++) {
224 fArrayQDY->AddAt( ((TH1F*)calibTracks.fArrayQDY->At(i)->Clone()), i);
225 fArrayQDZ->AddAt( ((TH1F*)calibTracks.fArrayQDZ->At(i)->Clone()), i);
226 fArrayQRMSY->AddAt( ((TH1F*)calibTracks.fArrayQRMSY->At(i)->Clone()), i);
227 fArrayQRMSZ->AddAt( ((TH1F*)calibTracks.fArrayQRMSZ->At(i)->Clone()), i);
230 (calibTracks.fResolY) ? length = calibTracks.fResolY->GetEntriesFast() : length = -1;
231 fResolY = new TObjArray(length);
232 fResolZ = new TObjArray(length);
233 fRMSY = new TObjArray(length);
234 fRMSZ = new TObjArray(length);
235 for (Int_t i = 0; i < length; i++) {
236 fResolY->AddAt( ((TH1F*)calibTracks.fResolY->At(i)->Clone()), i);
237 fResolZ->AddAt( ((TH1F*)calibTracks.fResolZ->At(i)->Clone()), i);
238 fRMSY->AddAt( ((TH1F*)calibTracks.fRMSY->At(i)->Clone()), i);
239 fRMSZ->AddAt( ((TH1F*)calibTracks.fRMSZ->At(i)->Clone()), i);
242 (calibTracks.fArrayChargeVsDriftlength) ? length = calibTracks.fArrayChargeVsDriftlength->GetEntriesFast() : length = -1;
243 (calibTracks.fArrayChargeVsDriftlength) ? fArrayChargeVsDriftlength = new TObjArray(length) : fArrayChargeVsDriftlength = 0;
244 for (Int_t i = 0; i < length; i++) {
245 fArrayChargeVsDriftlength->AddAt( ((TProfile*)calibTracks.fArrayChargeVsDriftlength->At(i)->Clone()), i);
248 fDeltaY = (TH1F*)calibTracks.fDeltaY->Clone();
249 fDeltaZ = (TH1F*)calibTracks.fDeltaZ->Clone();
250 fHclus = (TH1I*)calibTracks.fHclus->Clone();
251 fClusterCutHisto = (TH2I*)calibTracks.fClusterCutHisto->Clone();
252 fRejectedTracksHisto = (TH1I*)calibTracks.fRejectedTracksHisto->Clone();
253 fHclusterPerPadrow = (TH1I*)calibTracks.fHclusterPerPadrow->Clone();
254 fHclusterPerPadrowRaw = (TH1I*)calibTracks.fHclusterPerPadrowRaw->Clone();
255 fcalPadRegionChargeVsDriftlength = (AliTPCCalPadRegion*)calibTracks.fcalPadRegionChargeVsDriftlength->Clone();
256 fCalPadClusterPerPad = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPad->Clone();
257 fCalPadClusterPerPadRaw = (AliTPCCalPad*)calibTracks.fCalPadClusterPerPadRaw->Clone();
259 fCuts = new AliTPCcalibTracksCuts(calibTracks.fCuts->GetMinClusters(), calibTracks.fCuts->GetMinRatio(),
260 calibTracks.fCuts->GetMax1pt(), calibTracks.fCuts->GetEdgeYXCutNoise(), calibTracks.fCuts->GetEdgeThetaCutNoise());
261 SetNameTitle(calibTracks.GetName(), calibTracks.GetTitle());
262 TH1::AddDirectory(dirStatus); // set status back to original status
263 // cout << "+++++ end of copy constructor +++++" << endl; // TO BE REMOVED
267 AliTPCcalibTracks & AliTPCcalibTracks::operator=(const AliTPCcalibTracks& calibTracks){
269 // assgnment operator
271 if (this != &calibTracks) {
272 new (this) AliTPCcalibTracks(calibTracks);
279 AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, AliTPCClusterParam *clusterParam, AliTPCcalibTracksCuts* cuts, Int_t logLevel) :
283 fHisDeltaY(0), // THnSparse - delta Y
284 fHisDeltaZ(0), // THnSparse - delta Z
285 fHisRMSY(0), // THnSparse - rms Y
286 fHisRMSZ(0), // THnSparse - rms Z
287 fHisQmax(0), // THnSparse - qmax
288 fHisQtot(0), // THnSparse - qtot
295 fArrayChargeVsDriftlength(0),
296 fcalPadRegionChargeVsDriftlength(0),
305 fRejectedTracksHisto(0),
306 fHclusterPerPadrow(0),
307 fHclusterPerPadrowRaw(0),
309 fCalPadClusterPerPad(0),
310 fCalPadClusterPerPadRaw(0)
313 // AliTPCcalibTracks constructor
314 // specify 'name' and 'title' of your object
315 // specify 'clusterParam', (needed for TPC cluster error and shape parameterization)
316 // In the parameter 'cuts' the cuts are specified, that decide
317 // weather a track will be accepted for calibration or not.
319 // fDebugLevel - debug output: -1: silence, 0: default, 1: things like constructor called, 5: write fDebugStreamer, 6: waste your screen
321 // All histograms are instatiated in this constructor.
324 this->SetTitle(title);
326 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks' main constructor ***** " << endl;
327 G__SetCatchException(0);
329 fClusterParam = clusterParam;
331 fClusterParam->SetInstance(fClusterParam);
334 Error("AliTPCcalibTracks","No cluster parametrization found! A valid clusterParam object is needed in the constructor. (To be found in 'TPCClusterParam.root'.)");
337 SetDebugLevel(logLevel);
340 TH1::AddDirectory(kFALSE);
345 fHclus = new TH1I("hclus","Number of clusters per track",160, 0, 160); // valgrind 3
346 fRejectedTracksHisto = new TH1I("RejectedTracksHisto", "Rejected tracks, sorted by failed cut", 100, -1, 10);
347 fHclusterPerPadrow = new TH1I("fHclusterPerPadrow", " clusters per padRow, used for the resolution tree", 160, 0, 160);
348 fHclusterPerPadrowRaw = new TH1I("fHclusterPerPadrowRaw", " clusters per padRow, before cutting clusters", 160, 0, 160);
349 fCalPadClusterPerPad = new AliTPCCalPad("fCalPadClusterPerPad", "clusters per pad");
350 fCalPadClusterPerPadRaw = new AliTPCCalPad("fCalPadClusterPerPadRaw", "clusters per pad, before cutting clusters");
351 fClusterCutHisto = new TH2I("fClusterCutHisto", "Cutted cluster over padRow; Cut Criterium; PadRow", 5,1,5, 160,0,159);
353 // Amplitude - sector - row histograms
354 fArrayAmpRow = new TObjArray(72);
355 fArrayAmp = new TObjArray(72);
356 fArrayChargeVsDriftlength = new TObjArray(72);
358 for (Int_t i = 0; i < 36; i++){
359 snprintf(chname,100,"Amp_row_Sector%d",i);
360 prof1 = new TProfile(chname,chname,63,0,64); // valgrind 3 193,536 bytes in 354 blocks are still reachable
361 prof1->SetXTitle("Pad row");
362 prof1->SetYTitle("Mean Max amplitude");
363 fArrayAmpRow->AddAt(prof1,i);
364 snprintf(chname,100,"Amp_row_Sector%d",i+36);
365 prof1 = new TProfile(chname,chname,96,0,97); // valgrind 3 3,912 bytes in 6 blocks are possibly lost
366 prof1->SetXTitle("Pad row");
367 prof1->SetYTitle("Mean Max amplitude");
368 fArrayAmpRow->AddAt(prof1,i+36);
371 sprintf(chname,"Amp_Sector%d",i);
372 his1 = new TH1F(chname,chname,100,0,32); // valgrind
373 his1->SetXTitle("Max Amplitude (ADC)");
374 fArrayAmp->AddAt(his1,i);
375 sprintf(chname,"Amp_Sector%d",i+36);
376 his1 = new TH1F(chname,chname,100,0,32); // valgrind 3 13,408,208 bytes in 229 blocks are still reachable
377 his1->SetXTitle("Max Amplitude (ADC)");
378 fArrayAmp->AddAt(his1,i+36);
381 snprintf(chname,100, "driftlengt vs. charge, ROC %i", i);
382 prof1 = new TProfile(chname, chname, 25, 0, 250);
383 prof1->SetYTitle("Charge");
384 prof1->SetXTitle("Driftlength");
385 fArrayChargeVsDriftlength->AddAt(prof1,i);
386 snprintf(chname,100, "driftlengt vs. charge, ROC %i", i+36);
387 prof1 = new TProfile(chname, chname, 25, 0, 250);
388 prof1->SetYTitle("Charge");
389 prof1->SetXTitle("Driftlength");
390 fArrayChargeVsDriftlength->AddAt(prof1,i+36);
393 TH1::AddDirectory(kFALSE);
395 fDeltaY = new TH1F("DeltaY","DeltaY",100,-1,1);
396 fDeltaZ = new TH1F("DeltaZ","DeltaZ",100,-1,1);
398 fResolY = new TObjArray(3);
399 fResolZ = new TObjArray(3);
400 fRMSY = new TObjArray(3);
401 fRMSZ = new TObjArray(3);
404 his3D = new TH3F("Resol Y0","Resol Y0", 5,20,250, 4, 0,1., 50, -1,1);
405 fResolY->AddAt(his3D,0);
406 his3D = new TH3F("Resol Y1","Resol Y1", 5,20,250, 4, 0,1., 50, -1,1);
407 fResolY->AddAt(his3D,1);
408 his3D = new TH3F("Resol Y2","Resol Y2", 5,20,250, 4, 0,0.8, 50, -1,1);
409 fResolY->AddAt(his3D,2);
411 his3D = new TH3F("Resol Z0","Resol Z0", 5,20,250, 4, 0,1, 50, -1,1);
412 fResolZ->AddAt(his3D,0);
413 his3D = new TH3F("Resol Z1","Resol Z1", 5,20,250, 4, 0,1, 50, -1,1);
414 fResolZ->AddAt(his3D,1);
415 his3D = new TH3F("Resol Z2","Resol Z2", 5,20,250, 4, 0,1, 50, -1,1);
416 fResolZ->AddAt(his3D,2);
418 his3D = new TH3F("RMS Y0","RMS Y0", 5,20,250, 4, 0,1., 50, 0,0.8);
419 fRMSY->AddAt(his3D,0);
420 his3D = new TH3F("RMS Y1","RMS Y1", 5,20,250, 4, 0,1., 50, 0,0.8);
421 fRMSY->AddAt(his3D,1);
422 his3D = new TH3F("RMS Y2","RMS Y2", 5,20,250, 4, 0,0.8, 50, 0,0.8);
423 fRMSY->AddAt(his3D,2);
425 his3D = new TH3F("RMS Z0","RMS Z0", 5,20,250, 4, 0,1, 50, 0,0.8);
426 fRMSZ->AddAt(his3D,0);
427 his3D = new TH3F("RMS Z1","RMS Z1", 5,20,250, 4, 0,1, 50, 0,0.8);
428 fRMSZ->AddAt(his3D,1);
429 his3D = new TH3F("RMS Z2","RMS Z2", 5,20,250, 4, 0,1, 50, 0,0.8);
430 fRMSZ->AddAt(his3D,2);
433 TH1::AddDirectory(kFALSE);
435 fArrayQDY = new TObjArray(300);
436 fArrayQDZ = new TObjArray(300);
437 fArrayQRMSY = new TObjArray(300);
438 fArrayQRMSZ = new TObjArray(300);
439 for (Int_t iq = 0; iq <= 10; iq++){
440 for (Int_t ipad = 0; ipad < 3; ipad++){
441 Int_t bin = GetBin(iq, ipad);
442 Float_t qmean = GetQ(bin);
444 snprintf(hname,100,"ResolY Pad%d Qmiddle%f",ipad, qmean);
445 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
446 fArrayQDY->AddAt(his3D, bin);
447 snprintf(hname,100,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
448 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
449 fArrayQDZ->AddAt(his3D, bin);
450 snprintf(hname,100,"RMSY Pad%d Qmiddle%f",ipad, qmean);
451 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
452 fArrayQRMSY->AddAt(his3D, bin);
453 snprintf(hname,100,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
454 his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
455 fArrayQRMSZ->AddAt(his3D, bin);
459 fcalPadRegionChargeVsDriftlength = new AliTPCCalPadRegion("fcalPadRegionChargeVsDriftlength", "TProfiles with charge vs driftlength for each pad region");
461 for (UInt_t padSize = 0; padSize < 3; padSize++) {
462 for (UInt_t isector = 0; isector < 36; isector++) {
463 if (padSize == 0) sprintf(chname, "driftlengt vs. charge, sector %i, short pads", isector);
464 if (padSize == 1) sprintf(chname, "driftlengt vs. charge, sector %i, medium pads", isector);
465 if (padSize == 2) sprintf(chname, "driftlengt vs. charge, sector %i, long pads", isector);
466 tempProf = new TProfile(chname, chname, 500, 0, 250);
467 tempProf->SetYTitle("Charge");
468 tempProf->SetXTitle("Driftlength");
469 fcalPadRegionChargeVsDriftlength->SetObject(tempProf, isector, padSize);
474 if (GetDebugLevel() > 1) cout << "AliTPCcalibTracks object sucessfully constructed: " << GetName() << endl;
475 cout << "end of main constructor" << endl; // TO BE REMOVED
479 AliTPCcalibTracks::~AliTPCcalibTracks() {
481 // AliTPCcalibTracks destructor
484 if (GetDebugLevel() > 0) cout << "AliTPCcalibTracks' destuctor called." << endl;
486 if (fArrayAmpRow) length = fArrayAmpRow->GetEntriesFast();
487 for (Int_t i = 0; i < length; i++){
488 delete fArrayAmpRow->At(i);
489 delete fArrayAmp->At(i);
498 length = fResolY->GetEntriesFast();
499 for (Int_t i = 0; i < length; i++){
500 delete fResolY->At(i);
501 delete fResolZ->At(i);
512 length = fArrayQDY->GetEntriesFast();
513 for (Int_t i = 0; i < length; i++){
514 delete fArrayQDY->At(i);
515 delete fArrayQDZ->At(i);
516 delete fArrayQRMSY->At(i);
517 delete fArrayQRMSZ->At(i);
521 if (fArrayChargeVsDriftlength) {
522 length = fArrayChargeVsDriftlength->GetEntriesFast();
523 for (Int_t i = 0; i < length; i++){
524 delete fArrayChargeVsDriftlength->At(i);
533 delete fArrayChargeVsDriftlength;
536 delete fRejectedTracksHisto;
537 delete fClusterCutHisto;
538 delete fHclusterPerPadrow;
539 delete fHclusterPerPadrowRaw;
540 if (fCalPadClusterPerPad) delete fCalPadClusterPerPad;
541 if (fCalPadClusterPerPadRaw) delete fCalPadClusterPerPadRaw;
542 if(fcalPadRegionChargeVsDriftlength) {
543 fcalPadRegionChargeVsDriftlength->Delete();
544 delete fcalPadRegionChargeVsDriftlength;
546 delete fHisDeltaY; // THnSparse - delta Y
547 delete fHisDeltaZ; // THnSparse - delta Z
548 delete fHisRMSY; // THnSparse - rms Y
549 delete fHisRMSZ; // THnSparse - rms Z
550 delete fHisQmax; // THnSparse - qmax
551 delete fHisQtot; // THnSparse - qtot
557 void AliTPCcalibTracks::Process(AliTPCseed *track){
559 // To be called in the selector
560 // first AcceptTrack is evaluated, then calls all the following analyse functions:
561 // FillResolutionHistoLocal(track)
564 if (GetDebugLevel() > 5) Info("Process","Starting to process the track...");
565 Int_t accpetStatus = AcceptTrack(track);
566 if (accpetStatus == 0) {
567 FillResolutionHistoLocal(track);
569 else fRejectedTracksHisto->Fill(accpetStatus);
574 Int_t AliTPCcalibTracks::GetBin(Float_t q, Int_t pad){
576 // calculate bins for given q and pad type
579 Int_t res = TMath::Max( TMath::Nint((TMath::Sqrt(q) - 3.)), 0 );
586 Int_t AliTPCcalibTracks::GetBin(Int_t iq, Int_t pad){
588 // calculate bins for given iq and pad type
591 return iq * 3 + pad;;
595 Float_t AliTPCcalibTracks::GetQ(Int_t bin){
597 // returns to bin belonging charge
600 Int_t bin0 = bin / 3;
606 Float_t AliTPCcalibTracks::GetPad(Int_t bin){
608 // returns to bin belonging pad
616 Int_t AliTPCcalibTracks::AcceptTrack(AliTPCseed * track){
618 // Function, that decides wheather a given track is accepted for
619 // the analysis or not.
620 // The cuts are specified in the AliTPCcalibTracksCuts object 'fCuts'
621 // Returns 0 if a track is accepted or an integer different from 0
622 // to indicate the failed cut
624 const Int_t kMinClusters = fCuts->GetMinClusters();
625 const Float_t kMinRatio = fCuts->GetMinRatio();
626 const Float_t kMax1pt = fCuts->GetMax1pt();
627 const Float_t kEdgeYXCutNoise = fCuts->GetEdgeYXCutNoise();
628 const Float_t kEdgeThetaCutNoise = fCuts->GetEdgeThetaCutNoise();
631 // edge induced noise tracks - NEXT RELEASE will be removed during tracking
632 if ( TMath::Abs(track->GetY() / track->GetX()) > kEdgeYXCutNoise )
633 if ( TMath::Abs(track->GetTgl()) < kEdgeThetaCutNoise ) return 1;
634 if (track->GetNumberOfClusters() < kMinClusters) return 2;
635 Float_t ratio = track->GetNumberOfClusters() / (track->GetNFoundable() + 1.);
636 if (ratio < kMinRatio) return 3;
637 // Float_t mpt = track->Get1Pt(); // Get1Pt() doesn't exist any more
638 Float_t mpt = track->GetSigned1Pt();
639 if (TMath::Abs(mpt) > kMax1pt) return 4;
640 //if (TMath::Abs(track->GetZ())>240.) return kFALSE;
641 //if (TMath::Abs(track->GetZ())<10.) return kFALSE;
642 //if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
644 if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");
649 void AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
651 // fill resolution histograms - localy - tracklet in the neighborhood
652 // write debug information to 'TPCSelectorDebug.root'
654 // _ the main function, called during track analysis _
656 // loop over all padrows along the track
657 // fit tracklets (length: 13 clusters) calculate mean chi^2 for this track-fit in Y and Z direction
659 // loop again over all padrows along the track
660 // fit tracklet (clusters in given padrow +- kDelta padrows)
661 // with polynom of 2nd order and two polynoms of 1st order
662 // take both polynoms of 1st order, calculate difference of their parameters
663 // add covariance matrixes and calculate chi2 of this difference
664 // if this chi2 is bigger than a given threshold, assume that the current cluster is
665 // a kink an goto next padrow
667 // fill fArrayAmpRow, array with amplitudes vs. row for given sector
668 // fill fArrayAmp, array with amplitude histograms for give sector
669 // fill fRMSY, fRMSZ, fArrayQRMSY and fArrayQRMSZ, fDeltaY, fDeltaZ, fResolY, fResolZ, fArrayQDY, fArrayQDY
671 // write debug information to 'TPCSelectorDebug.root'
672 // only for every kDeltaWriteDebugStream'th padrow to reduce data volume
673 // and to avoid redundant data
676 static TLinearFitter fFitterLinY1(2,"pol1"); //
677 static TLinearFitter fFitterLinZ1(2,"pol1"); //
678 static TLinearFitter fFitterLinY2(2,"pol1"); //
679 static TLinearFitter fFitterLinZ2(2,"pol1"); //
680 static TLinearFitter fFitterParY(3,"pol2"); //
681 static TLinearFitter fFitterParZ(3,"pol2"); //
683 fFitterLinY1.StoreData(kFALSE);
684 fFitterLinZ1.StoreData(kFALSE);
685 fFitterLinY2.StoreData(kFALSE);
686 fFitterLinZ2.StoreData(kFALSE);
687 fFitterParY.StoreData(kFALSE);
688 fFitterParZ.StoreData(kFALSE);
691 if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
692 const Int_t kDelta = 10; // delta rows to fit
693 const Float_t kMinRatio = 0.75; // minimal ratio
694 // const Float_t kCutChi2 = 6.; // cut chi2 - left right - kink removal
695 const Float_t kErrorFraction = 0.5; // use only clusters with small interpolation error - for error param
696 const Int_t kFirstLargePad = 127; // medium pads -> long pads
697 const Float_t kLargePadSize = 1.5; // factor between medium and long pads' area
698 const Int_t kDeltaWriteDebugStream = 5; // only for every kDeltaWriteDebugStream'th padrow debug information is calulated and written to debugstream
705 TMatrixD matrixY0(2,2);
706 TMatrixD matrixZ0(2,2);
707 TMatrixD matrixY1(2,2);
708 TMatrixD matrixZ1(2,2);
710 // estimate mean error
711 Int_t nTrackletsAll = 0; // number of tracklets for given track
712 Float_t csigmaY = 0; // mean sigma for tracklet refit in Y direction
713 Float_t csigmaZ = 0; // mean sigma for tracklet refit in Z direction
714 Int_t nClusters = 0; // working variable, number of clusters per tracklet
715 Int_t sectorG = -1; // working variable, sector of tracklet, has to stay constant for one tracklet
717 fHclus->Fill(track->GetNumberOfClusters()); // for statistics overview
718 // ---------------------------------------------------------------------
719 for (Int_t irow = 0; irow < 159; irow++){
720 // loop over all rows along the track
721 // fit tracklets (length: 13 rows) with pol2 in Y and Z direction
722 // calculate mean chi^2 for this track-fit in Y and Z direction
723 AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
724 if (!cluster0) continue; // no cluster found
725 Int_t sector = cluster0->GetDetector();
726 fHclusterPerPadrowRaw->Fill(irow);
728 Int_t ipad = TMath::Nint(cluster0->GetPad());
729 Float_t value = fCalPadClusterPerPadRaw->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
730 fCalPadClusterPerPadRaw->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
732 if (sector != sectorG){
733 // track leaves sector before it crossed enough rows to fit / initialization
735 fFitterParY.ClearPoints();
736 fFitterParZ.ClearPoints();
741 Double_t x = cluster0->GetX();
742 fFitterParY.AddPoint(&x, cluster0->GetY(), 1);
743 fFitterParZ.AddPoint(&x, cluster0->GetZ(), 1);
745 if ( nClusters >= kDelta + 3 ){
746 // if more than 13 (kDelta+3) clusters were added to the fitters
747 // fit the tracklet, increase trackletCounter
751 csigmaY += fFitterParY.GetChisquare() / (nClusters - 3.);
752 csigmaZ += fFitterParZ.GetChisquare() / (nClusters - 3.);
754 fFitterParY.ClearPoints();
755 fFitterParZ.ClearPoints();
758 } // for (Int_t irow = 0; irow < 159; irow++)
759 // mean chi^2 for all tracklet fits in Y and in Z direction:
760 csigmaY = TMath::Sqrt(TMath::Abs(csigmaY) / (nTrackletsAll+0.1));
761 csigmaZ = TMath::Sqrt(TMath::Abs(csigmaZ) / (nTrackletsAll+0.1));
762 // ---------------------------------------------------------------------
764 for (Int_t irow = 0; irow < 159; irow++){
765 // loop again over all rows along the track
768 Int_t nclFound = 0; // number of clusters in the neighborhood
769 Int_t ncl0 = 0; // number of clusters in rows < rowOfCenterCluster
770 Int_t ncl1 = 0; // number of clusters in rows > rowOfCenterCluster
771 AliTPCclusterMI * cluster0 = track->GetClusterPointer(irow);
772 if (!cluster0) continue;
773 Int_t sector = cluster0->GetDetector();
774 Float_t xref = cluster0->GetX();
777 fFitterParY.ClearPoints();
778 fFitterParZ.ClearPoints();
779 fFitterLinY1.ClearPoints();
780 fFitterLinZ1.ClearPoints();
781 fFitterLinY2.ClearPoints();
782 fFitterLinZ2.ClearPoints();
784 // fit tracklet (clusters in given padrow +- kDelta padrows)
785 // with polynom of 2nd order and two polynoms of 1st order
786 // take both polynoms of 1st order, calculate difference of their parameters
787 // add covariance matrixes and calculate chi2 of this difference
788 // if this chi2 is bigger than a given threshold, assume that the current cluster is
789 // a kink an goto next padrow
791 for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
792 // loop over irow +- kDelta rows (neighboured rows)
795 if (idelta == 0) continue; // don't use center cluster
796 if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
797 AliTPCclusterMI * currentCluster = track->GetClusterPointer(irow + idelta);
798 if (!currentCluster) continue;
799 if (currentCluster->GetType() < 0) continue;
800 if (currentCluster->GetDetector() != sector) continue;
801 Double_t x = currentCluster->GetX() - xref; // x = differece: current cluster - cluster @ irow
805 fFitterLinY1.AddPoint(&x, currentCluster->GetY(), csigmaY);
806 fFitterLinZ1.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
810 fFitterLinY2.AddPoint(&x, currentCluster->GetY(), csigmaY);
811 fFitterLinZ2.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
813 fFitterParY.AddPoint(&x, currentCluster->GetY(), csigmaY);
814 fFitterParZ.AddPoint(&x, currentCluster->GetZ(), csigmaZ);
815 } // loop over neighbourhood for fitter filling
819 if (nclFound < kDelta * kMinRatio) fRejectedTracksHisto->Fill(10);
820 if (nclFound < kDelta * kMinRatio) fClusterCutHisto->Fill(1, irow);
821 if (nclFound < kDelta * kMinRatio) continue; // if not enough clusters (7.5) found in neighbourhood goto next padrow
824 Double_t chi2 = (fFitterParY.GetChisquare() + fFitterParZ.GetChisquare()) / (2. * nclFound - 6.);
825 TTreeSRedirector *cstream = GetDebugStreamer();
832 // only when there are enough clusters (4) in each direction
842 if (ncl0 > 4 && ncl1 > 4){
843 fFitterLinY1.GetCovarianceMatrix(matrixY0);
844 fFitterLinY2.GetCovarianceMatrix(matrixY1);
845 fFitterLinZ1.GetCovarianceMatrix(matrixZ0);
846 fFitterLinZ2.GetCovarianceMatrix(matrixZ1);
847 fFitterLinY2.GetParameters(paramY1);
848 fFitterLinZ2.GetParameters(paramZ1);
849 fFitterLinY1.GetParameters(paramY0);
850 fFitterLinZ1.GetParameters(paramZ0);
853 matrixY0 += matrixY1;
854 matrixZ0 += matrixZ1;
857 TMatrixD difY(2, 1, paramY0.GetMatrixArray());
858 TMatrixD difYT(1, 2, paramY0.GetMatrixArray());
860 TMatrixD mulY(matrixY0, TMatrixD::kMult, difY);
861 TMatrixD chi2Y(difYT, TMatrixD::kMult, mulY);
862 cchi2 += chi2Y(0, 0);
864 TMatrixD difZ(2, 1, paramZ0.GetMatrixArray());
865 TMatrixD difZT(1, 2, paramZ0.GetMatrixArray());
867 TMatrixD mulZ(matrixZ0, TMatrixD::kMult, difZ);
868 TMatrixD chi2Z(difZT, TMatrixD::kMult, mulZ);
869 cchi2 += chi2Z(0, 0);
878 // current padrow has no kink
880 // get fit parameters from pol2 fit:
881 Double_t paramY[4], paramZ[4];
882 paramY[0] = fFitterParY.GetParameter(0);
883 paramY[1] = fFitterParY.GetParameter(1);
884 paramY[2] = fFitterParY.GetParameter(2);
885 paramZ[0] = fFitterParZ.GetParameter(0);
886 paramZ[1] = fFitterParZ.GetParameter(1);
887 paramZ[2] = fFitterParZ.GetParameter(2);
889 Double_t tracky = paramY[0];
890 Double_t trackz = paramZ[0];
891 Float_t deltay = tracky - cluster0->GetY();
892 Float_t deltaz = trackz - cluster0->GetZ();
893 Float_t angley = paramY[1] - paramY[0] / xref;
894 Float_t anglez = paramZ[1];
896 Float_t max = cluster0->GetMax();
897 UInt_t isegment = cluster0->GetDetector() % 36;
898 Int_t padSize = 0; // short pads
899 if (cluster0->GetDetector() >= 36) {
900 padSize = 1; // medium pads
901 if (cluster0->GetRow() > 63) padSize = 2; // long pads
904 // =========================================
905 // wirte collected information to histograms
906 // =========================================
908 TProfile *profAmpRow = (TProfile*)fArrayAmpRow->At(sector);
909 if ( irow >= kFirstLargePad) max /= kLargePadSize;
910 Double_t smax = TMath::Sqrt(max);
911 profAmpRow->Fill( (Double_t)cluster0->GetRow(), smax );
912 TH1F *hisAmp = (TH1F*)fArrayAmp->At(sector);
915 // remove the following two lines one day:
916 TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector);
917 profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
919 TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize));
920 profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
922 fHclusterPerPadrow->Fill(irow); // fill histogram showing clusters per padrow
923 Int_t ipad = TMath::Nint(cluster0->GetPad());
924 Float_t value = fCalPadClusterPerPad->GetCalROC(sector)->GetValue((sector<36)?irow:irow-64, TMath::Nint(cluster0->GetPad()));
925 fCalPadClusterPerPad->GetCalROC(sector)->SetValue((sector<36)?irow:irow-64, ipad, value + 1 );
929 his3 = (TH3F*)fRMSY->At(padSize);
930 if (his3) his3->Fill(250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
931 his3 = (TH3F*)fRMSZ->At(padSize);
932 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
934 his3 = (TH3F*)fArrayQRMSY->At(GetBin(cluster0->GetMax(), padSize));
935 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), TMath::Sqrt(cluster0->GetSigmaY2()) );
936 his3 = (TH3F*)fArrayQRMSZ->At(GetBin(cluster0->GetMax(), padSize));
937 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), TMath::Sqrt(cluster0->GetSigmaZ2()) );
940 // Fill resolution histograms
941 Bool_t useForResol = kTRUE;
942 if (fFitterParY.GetParError(0) > kErrorFraction * csigmaY) useForResol = kFALSE;
945 Float_t zdrift = 250 - TMath::Abs(cluster0->GetZ());
946 Float_t sy = cluster0->GetSigmaY2();
947 Float_t sz = cluster0->GetSigmaZ2();
948 (*cstream)<<"Resol0"<<
949 "run="<<fRun<< // run number
950 "event="<<fEvent<< // event number
951 "time="<<fTime<< // time stamp of event
952 "trigger="<<fTrigger<< // trigger
953 "mag="<<fMagF<< // magnetic field
954 "padSize="<<padSize<<
966 fDeltaY->Fill(deltay);
967 fDeltaZ->Fill(deltaz);
968 his3 = (TH3F*)fResolY->At(padSize);
969 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(angley), deltay );
970 his3 = (TH3F*)fResolZ->At(padSize);
971 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()), TMath::Abs(anglez), deltaz );
972 his3 = (TH3F*)fArrayQDY->At(GetBin(cluster0->GetMax(), padSize));
973 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(angley), deltay );
974 his3 = (TH3F*)fArrayQDZ->At(GetBin(cluster0->GetMax(), padSize));
975 if (his3) his3->Fill( 250 - TMath::Abs(cluster0->GetZ()),TMath::Abs(anglez), deltaz );
978 //=============================================================================================
980 if (useForResol && nclFound > 2 * kMinRatio * kDelta
981 && irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){
982 if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
983 FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta);
984 } // if (useForResol && nclFound > 2 * kMinRatio * kDelta)
986 // Fill THN histograms
990 xvar[2]=1.-TMath::Abs(cluster0->GetZ())/250.;
991 xvar[3]=cluster0->GetMax();
995 xvar[4]=cluster0->GetPad()-TMath::Nint(cluster0->GetPad());
997 fHisDeltaY->Fill(xvar);
998 xvar[0]=TMath::Sqrt(cluster0->GetSigmaY2());
999 fHisRMSY->Fill(xvar);
1001 xvar[4]=cluster0->GetTimeBin()-TMath::Nint(cluster0->GetTimeBin());
1003 fHisDeltaZ->Fill(xvar);
1004 xvar[0]=TMath::Sqrt(cluster0->GetSigmaZ2());
1005 fHisRMSZ->Fill(xvar);
1007 } // loop over all padrows along the track: for (Int_t irow = 0; irow < 159; irow++)
1008 } // FillResolutionHistoLocal(...)
1012 void AliTPCcalibTracks::FillResolutionHistoLocalDebugPart(AliTPCseed *track, AliTPCclusterMI *cluster0, Int_t irow, Float_t angley, Float_t anglez, Int_t nclFound, Int_t kDelta) {
1014 // - debug part of FillResolutionHistoLocal -
1015 // called only for every kDeltaWriteDebugStream'th padrow, to avoid to much redundant data
1016 // called only for GetStreamLevel() > 4
1017 // fill resolution trees
1020 Int_t sector = cluster0->GetDetector();
1021 Float_t xref = cluster0->GetX();
1022 Int_t padSize = 0; // short pads
1023 if (cluster0->GetDetector() >= 36) {
1024 padSize = 1; // medium pads
1025 if (cluster0->GetRow() > 63) padSize = 2; // long pads
1028 static TLinearFitter fitY0(3, "pol2");
1029 static TLinearFitter fitZ0(3, "pol2");
1030 static TLinearFitter fitY2(5, "hyp4");
1031 static TLinearFitter fitZ2(5, "hyp4");
1032 static TLinearFitter fitY2Q(5, "hyp4");
1033 static TLinearFitter fitZ2Q(5, "hyp4");
1034 static TLinearFitter fitY2S(5, "hyp4");
1035 static TLinearFitter fitZ2S(5, "hyp4");
1036 fitY0.ClearPoints();
1037 fitZ0.ClearPoints();
1038 fitY2.ClearPoints();
1039 fitZ2.ClearPoints();
1040 fitY2Q.ClearPoints();
1041 fitZ2Q.ClearPoints();
1042 fitY2S.ClearPoints();
1043 fitZ2S.ClearPoints();
1045 for (Int_t idelta = -kDelta; idelta <= kDelta; idelta++){
1046 // loop over irow +- kDelta rows (neighboured rows)
1049 if (idelta == 0) continue;
1050 if (idelta + irow < 0 || idelta + irow > 159) continue; // don't go out of ROC
1051 AliTPCclusterMI * cluster = track->GetClusterPointer(irow + idelta);
1052 if (!cluster) continue;
1053 if (cluster->GetType() < 0) continue;
1054 if (cluster->GetDetector() != sector) continue;
1055 Double_t x = cluster->GetX() - xref;
1056 Double_t sigmaY0 = fClusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley) );
1057 Double_t sigmaZ0 = fClusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez) );
1059 Double_t sigmaYQ = fClusterParam->GetErrorQPar( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1060 Double_t sigmaZQ = fClusterParam->GetErrorQPar( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1061 Double_t sigmaYS = fClusterParam->GetErrorQParScaled( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
1062 TMath::Abs(angley), TMath::Abs(cluster->GetMax()) );
1063 Double_t sigmaZS = fClusterParam->GetErrorQParScaled( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())),
1064 TMath::Abs(anglez), TMath::Abs(cluster->GetMax()) );
1065 Float_t rmsYFactor = fClusterParam->GetShapeFactor( 0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1066 TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1067 TMath::Sqrt(cluster0->GetSigmaY2()), 0 );
1068 Float_t rmsZFactor = fClusterParam->GetShapeFactor(0, padSize,(250.0 - TMath::Abs(cluster->GetZ())),
1069 TMath::Abs(anglez), TMath::Abs(cluster->GetMax()),
1070 TMath::Sqrt(cluster0->GetSigmaZ2()),0 );
1071 sigmaYS = TMath::Sqrt(sigmaYS * sigmaYS + rmsYFactor * rmsYFactor / 12.);
1072 sigmaZS = TMath::Sqrt(sigmaZS * sigmaZS + rmsZFactor * rmsZFactor / 12. + rmsYFactor * rmsYFactor / 24.);
1075 fitY0.AddPoint(&x, cluster->GetY(), sigmaY0);
1076 fitZ0.AddPoint(&x, cluster->GetZ(), sigmaZ0);
1079 xxx[0] = ( (idelta+irow) % 2 == 0 ) ? 1 : 0;
1081 xxx[2] = ( (idelta+irow) % 2 == 0 ) ? x : 0;
1083 fitY2.AddPoint(xxx, cluster->GetY(), sigmaY0);
1084 fitY2Q.AddPoint(xxx, cluster->GetY(), sigmaYQ);
1085 fitY2S.AddPoint(xxx, cluster->GetY(), sigmaYS);
1086 fitZ2.AddPoint(xxx, cluster->GetZ(), sigmaZ0);
1087 fitZ2Q.AddPoint(xxx, cluster->GetZ(), sigmaZQ);
1088 fitZ2S.AddPoint(xxx, cluster->GetZ(), sigmaZS);
1090 } // neigbouhood-loop
1100 Float_t chi2Y0 = fitY0.GetChisquare() / (nclFound-3.);
1101 Float_t chi2Z0 = fitZ0.GetChisquare() / (nclFound-3.);
1102 Float_t chi2Y2 = fitY2.GetChisquare() / (nclFound-5.);
1103 Float_t chi2Z2 = fitZ2.GetChisquare() / (nclFound-5.);
1104 Float_t chi2Y2Q = fitY2Q.GetChisquare() / (nclFound-5.);
1105 Float_t chi2Z2Q = fitZ2Q.GetChisquare() / (nclFound-5.);
1106 Float_t chi2Y2S = fitY2S.GetChisquare() / (nclFound-5.);
1107 Float_t chi2Z2S = fitZ2S.GetChisquare() / (nclFound-5.);
1109 static TVectorD parY0(3);
1110 static TMatrixD matY0(3, 3);
1111 static TVectorD parZ0(3);
1112 static TMatrixD matZ0(3, 3);
1113 fitY0.GetParameters(parY0);
1114 fitY0.GetCovarianceMatrix(matY0);
1115 fitZ0.GetParameters(parZ0);
1116 fitZ0.GetCovarianceMatrix(matZ0);
1118 static TVectorD parY2(5);
1119 static TMatrixD matY2(5,5);
1120 static TVectorD parZ2(5);
1121 static TMatrixD matZ2(5,5);
1122 fitY2.GetParameters(parY2);
1123 fitY2.GetCovarianceMatrix(matY2);
1124 fitZ2.GetParameters(parZ2);
1125 fitZ2.GetCovarianceMatrix(matZ2);
1127 static TVectorD parY2Q(5);
1128 static TMatrixD matY2Q(5,5);
1129 static TVectorD parZ2Q(5);
1130 static TMatrixD matZ2Q(5,5);
1131 fitY2Q.GetParameters(parY2Q);
1132 fitY2Q.GetCovarianceMatrix(matY2Q);
1133 fitZ2Q.GetParameters(parZ2Q);
1134 fitZ2Q.GetCovarianceMatrix(matZ2Q);
1135 static TVectorD parY2S(5);
1136 static TMatrixD matY2S(5,5);
1137 static TVectorD parZ2S(5);
1138 static TMatrixD matZ2S(5,5);
1139 fitY2S.GetParameters(parY2S);
1140 fitY2S.GetCovarianceMatrix(matY2S);
1141 fitZ2S.GetParameters(parZ2S);
1142 fitZ2S.GetCovarianceMatrix(matZ2S);
1143 Float_t sigmaY0 = TMath::Sqrt(matY0(0,0));
1144 Float_t sigmaZ0 = TMath::Sqrt(matZ0(0,0));
1145 Float_t sigmaDY0 = TMath::Sqrt(matY0(1,1));
1146 Float_t sigmaDZ0 = TMath::Sqrt(matZ0(1,1));
1147 Float_t sigmaY2 = TMath::Sqrt(matY2(1,1));
1148 Float_t sigmaZ2 = TMath::Sqrt(matZ2(1,1));
1149 Float_t sigmaDY2 = TMath::Sqrt(matY2(3,3));
1150 Float_t sigmaDZ2 = TMath::Sqrt(matZ2(3,3));
1151 Float_t sigmaY2Q = TMath::Sqrt(matY2Q(1,1));
1152 Float_t sigmaZ2Q = TMath::Sqrt(matZ2Q(1,1));
1153 Float_t sigmaDY2Q = TMath::Sqrt(matY2Q(3,3));
1154 Float_t sigmaDZ2Q = TMath::Sqrt(matZ2Q(3,3));
1155 Float_t sigmaY2S = TMath::Sqrt(matY2S(1,1));
1156 Float_t sigmaZ2S = TMath::Sqrt(matZ2S(1,1));
1157 Float_t sigmaDY2S = TMath::Sqrt(matY2S(3,3));
1158 Float_t sigmaDZ2S = TMath::Sqrt(matZ2S(3,3));
1161 Float_t csigmaY0 = fClusterParam->GetError0Par(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(angley));
1162 Float_t csigmaZ0 = fClusterParam->GetError0Par(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),TMath::Abs(anglez));
1164 Float_t csigmaYQ = fClusterParam->GetErrorQPar(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1165 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1166 Float_t csigmaZQ = fClusterParam->GetErrorQPar(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1167 TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1168 Float_t csigmaYS = fClusterParam->GetErrorQParScaled(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1169 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1170 Float_t csigmaZS = fClusterParam->GetErrorQParScaled(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1171 TMath::Abs(anglez),TMath::Abs(cluster0->GetMax()));
1175 Float_t meanRMSY = 0;
1176 Float_t meanRMSZ = 0;
1178 for (Int_t idelta = -2; idelta <= 2; idelta++){
1179 // loop over neighbourhood
1180 if (idelta+irow < 0 || idelta+irow > 159) continue;
1181 // if (idelta+irow>159) continue;
1182 AliTPCclusterMI * cluster = track->GetClusterPointer(irow+idelta);
1183 if (!cluster) continue;
1184 meanRMSY += TMath::Sqrt(cluster->GetSigmaY2());
1185 meanRMSZ += TMath::Sqrt(cluster->GetSigmaZ2());
1191 Float_t rmsY = TMath::Sqrt(cluster0->GetSigmaY2());
1192 Float_t rmsZ = TMath::Sqrt(cluster0->GetSigmaZ2());
1193 Float_t rmsYT = fClusterParam->GetRMSQ(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1194 TMath::Abs(angley), TMath::Abs(cluster0->GetMax()));
1195 Float_t rmsZT = fClusterParam->GetRMSQ(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1196 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1197 Float_t rmsYT0 = fClusterParam->GetRMS0(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1198 TMath::Abs(angley));
1199 Float_t rmsZT0 = fClusterParam->GetRMS0(1,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1200 TMath::Abs(anglez));
1201 Float_t rmsYSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1202 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1203 Float_t rmsZSigma = fClusterParam->GetRMSSigma(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1204 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()));
1205 Float_t rmsYFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1206 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1208 Float_t rmsZFactor = fClusterParam->GetShapeFactor(0,padSize,(250.0-TMath::Abs(cluster0->GetZ())),
1209 TMath::Abs(anglez), TMath::Abs(cluster0->GetMax()),
1213 TTreeSRedirector *cstream = GetDebugStreamer();
1215 (*cstream)<<"ResolCl"<< // valgrind 3 40,000 bytes in 1 blocks are possibly
1216 "run="<<fRun<< // run number
1217 "event="<<fEvent<< // event number
1218 "time="<<fTime<< // time stamp of event
1219 "trigger="<<fTrigger<< // trigger
1220 "mag="<<fMagF<< // magnetic field
1223 "CSigmaY0="<<csigmaY0<< // cluster errorY
1224 "CSigmaYQ="<<csigmaYQ<< // cluster errorY - q scaled
1225 "CSigmaYS="<<csigmaYS<< // cluster errorY - q scaled
1226 "CSigmaZ0="<<csigmaZ0<< //
1227 "CSigmaZQ="<<csigmaZQ<<
1228 "CSigmaZS="<<csigmaZS<<
1229 "shapeYF="<<rmsYFactor<<
1230 "shapeZF="<<rmsZFactor<<
1233 "rmsYM="<<meanRMSY<<
1234 "rmsZM="<<meanRMSZ<<
1239 "rmsYS="<<rmsYSigma<<
1240 "rmsZS="<<rmsZSigma<<
1241 "padSize="<<padSize<<
1245 "SigmaY0="<<sigmaY0<<
1246 "SigmaZ0="<<sigmaZ0<<
1252 (*cstream)<<"ResolTr"<<
1253 "run="<<fRun<< // run number
1254 "event="<<fEvent<< // event number
1255 "time="<<fTime<< // time stamp of event
1256 "trigger="<<fTrigger<< // trigger
1257 "mag="<<fMagF<< // magnetic field
1258 "padSize="<<padSize<<
1266 "chi2Y2Q="<<chi2Y2Q<<
1267 "chi2Z2Q="<<chi2Z2Q<<
1268 "chi2Y2S="<<chi2Y2S<<
1269 "chi2Z2S="<<chi2Z2S<<
1278 "SigmaY0="<<sigmaY0<<
1279 "SigmaZ0="<<sigmaZ0<<
1280 "SigmaDY0="<<sigmaDY0<<
1281 "SigmaDZ0="<<sigmaDZ0<<
1282 "SigmaY2="<<sigmaY2<<
1283 "SigmaZ2="<<sigmaZ2<<
1284 "SigmaDY2="<<sigmaDY2<<
1285 "SigmaDZ2="<<sigmaDZ2<<
1286 "SigmaY2Q="<<sigmaY2Q<<
1287 "SigmaZ2Q="<<sigmaZ2Q<<
1288 "SigmaDY2Q="<<sigmaDY2Q<<
1289 "SigmaDZ2Q="<<sigmaDZ2Q<<
1290 "SigmaY2S="<<sigmaY2S<<
1291 "SigmaZ2S="<<sigmaZ2S<<
1292 "SigmaDY2S="<<sigmaDY2S<<
1293 "SigmaDZ2S="<<sigmaDZ2S<<
1304 TH2D * AliTPCcalibTracks::MakeDiff(TH2D * hfit, TF2 * func){
1306 // creates a new histogram which contains the difference between
1307 // the histogram hfit and the function func
1309 TH2D * result = (TH2D*)hfit->Clone(); // valgrind 3 40,139 bytes in 11 blocks are still reachable
1310 result->SetTitle(Form("%s fit residuals",result->GetTitle()));
1311 result->SetName(Form("%s fit residuals",result->GetName()));
1312 TAxis *xaxis = hfit->GetXaxis();
1313 TAxis *yaxis = hfit->GetYaxis();
1315 for (Int_t biny = 0; biny <= yaxis->GetNbins(); biny++) {
1316 x[1] = yaxis->GetBinCenter(biny);
1317 for (Int_t binx = 0; binx <= xaxis->GetNbins(); binx++) {
1318 x[0] = xaxis->GetBinCenter(binx);
1319 Int_t bin = hfit->GetBin(binx, biny);
1320 Double_t val = hfit->GetBinContent(bin);
1321 // result->SetBinContent( bin, (val - func->Eval(x[0], x[1])) / func->Eval(x[0], x[1]) );
1322 result->SetBinContent( bin, (val / func->Eval(x[0], x[1])) - 1 );
1329 void AliTPCcalibTracks::SetStyle() const {
1331 // set style, can be called by all draw functions
1333 gROOT->SetStyle("Plain");
1334 gStyle->SetFillColor(10);
1335 gStyle->SetPadColor(10);
1336 gStyle->SetCanvasColor(10);
1337 gStyle->SetStatColor(10);
1338 gStyle->SetPalette(1,0);
1339 gStyle->SetNumberContours(60);
1343 void AliTPCcalibTracks::Draw(Option_t* opt){
1345 // draw-function of AliTPCcalibTracks
1346 // will draws some exemplaric pictures
1349 if (GetDebugLevel() > 6) Info("Draw", "Drawing an exemplaric picture.");
1353 TCanvas *c1 = new TCanvas();
1355 TVirtualPad *upperThird = c1->GetPad(1);
1356 TVirtualPad *middleThird = c1->GetPad(2);
1357 TVirtualPad *lowerThird = c1->GetPad(3);
1358 upperThird->Divide(2,0);
1359 TVirtualPad *upleft = upperThird->GetPad(1);
1360 TVirtualPad *upright = upperThird->GetPad(2);
1361 middleThird->Divide(2,0);
1362 TVirtualPad *middleLeft = middleThird->GetPad(1);
1363 TVirtualPad *middleRight = middleThird->GetPad(2);
1364 lowerThird->Divide(2,0);
1365 TVirtualPad *downLeft = lowerThird->GetPad(1);
1366 TVirtualPad *downRight = lowerThird->GetPad(2);
1370 min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1371 max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1372 fDeltaY->SetAxisRange(min, max);
1373 fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable
1377 max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1378 min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1379 fDeltaZ->SetAxisRange(min, max);
1380 fDeltaZ->Fit("gaus","q","",min, max);
1387 fRejectedTracksHisto->Draw(opt);
1388 TPaveText *pt = new TPaveText(0.6,0.6, 0.8,0.8, "NDC");
1389 TText *t1 = pt->AddText("1: kEdgeThetaCutNoise");
1390 TText *t2 = pt->AddText("2: kMinClusters");
1391 TText *t3 = pt->AddText("3: kMinRatio");
1392 TText *t4 = pt->AddText("4: kMax1pt");
1393 t1 = t1; t2 = t2; t3 = t3; t4 = t4; // avoid compiler warnings
1394 pt->SetToolTipText("Legend for failed cuts");
1398 fHclusterPerPadrowRaw->Draw(opt);
1401 fHclusterPerPadrow->Draw(opt);
1405 void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){
1407 // all functions are called, that produce pictures
1408 // the histograms are written to the directory 'pathName'
1409 // 'stat' is a threshhold: only histograms with more than 'stat' entries are wirtten to file
1410 // 'stat' is also the number of minEntries for MakeResPlotsQTree
1413 if (GetDebugLevel() > 0) Info("MakeReport","Writing plots and trees to '%s'.", pathName);
1414 MakeAmpPlots(stat, pathName);
1415 MakeDeltaPlots(pathName);
1416 FitResolutionNew(pathName);
1417 FitRMSNew(pathName);
1418 MakeChargeVsDriftLengthPlots(pathName);
1419 // MakeResPlotsQ(1, 1);
1420 MakeResPlotsQTree(stat, pathName);
1424 void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, const char* pathName){
1426 // creates several plots:
1427 // fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps
1428 // fArrayAmp.ps: one histogram per sector, the histogram shows the charge per cluster
1429 // fArrayAmpRow.ps: one histogram per sector, mean max. amplitude vs. pad row with landau fit
1430 // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1431 // Empty histograms (sectors without data) are not written to file
1432 // the ps-files are written to the directory 'pathName', that is created if it does not exist
1433 // 'stat': only histograms with more than 'stat' entries are written to file.
1437 gSystem->MakeDirectory(pathName);
1438 gSystem->ChangeDirectory(pathName);
1440 TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1442 // histograms with accumulated amplitude for all IROCs and OROCs
1443 TH1F *allAmpHisIROC = ((TH1F*)(fArrayAmp->At(0))->Clone());
1444 allAmpHisIROC->SetName("Amp all IROCs");
1445 allAmpHisIROC->SetTitle("Amp all IROCs");
1446 TH1F *allAmpHisOROC = ((TH1F*)(fArrayAmp->At(36))->Clone());
1447 allAmpHisOROC->SetName("Amp all OROCs");
1448 allAmpHisOROC->SetTitle("Amp all OROCs");
1450 ps = new TPostScript("fArrayAmp.ps", 112);
1451 if (GetDebugLevel() > -1) cout << "creating fArrayAmp.ps..." << endl;
1452 for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++){
1453 if ( ((TH1F*)fArrayAmp->At(i))->GetEntries() < stat ) continue;
1455 ((TH1F*)fArrayAmp->At(i))->Draw();
1456 c1->Update(); // valgrind 3
1457 if (i > 0 && i < 36) {
1458 allAmpHisIROC->Add(((TH1F*)fArrayAmp->At(i)));
1459 allAmpHisOROC->Add(((TH1F*)fArrayAmp->At(i+36)));
1463 allAmpHisIROC->Draw();
1464 c1->Update(); // valgrind
1466 allAmpHisOROC->Draw();
1474 ps = new TPostScript("fArrayAmpRow.ps", 112);
1475 if (GetDebugLevel() > -1) cout << "creating fArrayAmpRow.ps..." << endl;
1476 for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++){
1477 his = (TH1F*)fArrayAmpRow->At(i);
1478 if (his->GetEntries() < stat) continue;
1480 min = TMath::Max( his->GetBinCenter(his->GetMaximumBin() )-100., 0.);
1481 max = his->GetBinCenter(5*his->GetMaximumBin()) + 100;
1482 his->SetAxisRange(min, max);
1483 his->Fit("pol3", "q", "", min, max);
1484 // his->Draw("error"); // don't use this line when you don't want to have empty pages in the ps-file
1490 gSystem->ChangeDirectory("..");
1494 void AliTPCcalibTracks::MakeDeltaPlots(const char* pathName){
1496 // creates several plots:
1497 // DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
1498 // the ps-files are written to the directory 'pathName', that is created if it does not exist
1502 gSystem->MakeDirectory(pathName);
1503 gSystem->ChangeDirectory(pathName);
1505 TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1510 ps = new TPostScript("DeltaYZ.ps", 112);
1511 if (GetDebugLevel() > -1) cout << "creating DeltaYZ.ps..." << endl;
1512 min = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())-20;
1513 max = fDeltaY->GetBinCenter(fDeltaY->GetMaximumBin())+20;
1514 fDeltaY->SetAxisRange(min, max);
1516 fDeltaY->Fit("gaus","q","",min, max); // valgrind 3 7 block possibly lost 2,400 bytes in 1 blocks are still reachable
1519 max = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())+20;
1520 min = fDeltaZ->GetBinCenter(fDeltaZ->GetMaximumBin())-20;
1521 fDeltaZ->SetAxisRange(min, max);
1522 fDeltaZ->Fit("gaus","q","",min, max);
1527 gSystem->ChangeDirectory("..");
1531 void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(const char* pathName){
1533 // creates charge vs. driftlength plots, one TProfile for each ROC
1534 // is not correct like this, should be one TProfile for each sector and padsize
1538 gSystem->MakeDirectory(pathName);
1539 gSystem->ChangeDirectory(pathName);
1541 TCanvas* c1 = new TCanvas(); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1543 ps = new TPostScript("chargeVsDriftlengthOld.ps", 112);
1544 if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlength.ps..." << endl;
1545 TProfile *chargeVsDriftlengthAllIROCs = ((TProfile*)fArrayChargeVsDriftlength->At(0)->Clone());
1546 TProfile *chargeVsDriftlengthAllOROCs = ((TProfile*)fArrayChargeVsDriftlength->At(36)->Clone());
1547 chargeVsDriftlengthAllIROCs->SetName("allAmpHisIROC");
1548 chargeVsDriftlengthAllIROCs->SetTitle("charge vs. driftlength, all IROCs");
1549 chargeVsDriftlengthAllOROCs->SetName("allAmpHisOROC");
1550 chargeVsDriftlengthAllOROCs->SetTitle("charge vs. driftlength, all OROCs");
1552 for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) {
1553 ((TProfile*)fArrayChargeVsDriftlength->At(i))->Draw();
1555 if (i > 0 && i < 36) {
1556 chargeVsDriftlengthAllIROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i)));
1557 chargeVsDriftlengthAllOROCs->Add(((TProfile*)fArrayChargeVsDriftlength->At(i+36)));
1561 chargeVsDriftlengthAllIROCs->Draw();
1562 c1->Update(); // valgrind
1564 chargeVsDriftlengthAllOROCs->Draw();
1569 gSystem->ChangeDirectory("..");
1573 void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(const char* pathName){
1575 // creates charge vs. driftlength plots, one TProfile for each ROC
1576 // under development....
1580 gSystem->MakeDirectory(pathName);
1581 gSystem->ChangeDirectory(pathName);
1583 TCanvas* c1 = new TCanvas("c1", "c1", 700,(Int_t)(TMath::Sqrt(2)*700)); // valgrind 3 ??? 634 bytes in 28 blocks are still reachable
1584 // TCanvas c1("c1", "c1", 500,(sqrt(2)*500))
1587 ps = new TPostScript("chargeVsDriftlength.ps", 111);
1588 if (GetDebugLevel() > -1) cout << "creating chargeVsDriftlengthNew.ps..." << endl;
1590 TProfile *chargeVsDriftlengthAllShortPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0)->Clone());
1591 TProfile *chargeVsDriftlengthAllMediumPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1)->Clone());
1592 TProfile *chargeVsDriftlengthAllLongPads = ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2)->Clone());
1593 chargeVsDriftlengthAllShortPads->SetName("allAmpHisShortPads");
1594 chargeVsDriftlengthAllShortPads->SetTitle("charge vs. driftlength, all sectors, short pads");
1595 chargeVsDriftlengthAllMediumPads->SetName("allAmpHisMediumPads");
1596 chargeVsDriftlengthAllMediumPads->SetTitle("charge vs. driftlength, all sectors, medium pads");
1597 chargeVsDriftlengthAllLongPads->SetName("allAmpHisLongPads");
1598 chargeVsDriftlengthAllLongPads->SetTitle("charge vs. driftlength, all sectors, long pads");
1600 for (Int_t i = 0; i < 36; i++) {
1601 c1->cd(1)->SetGridx();
1602 c1->cd(1)->SetGridy();
1603 ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,0))->Draw();
1604 c1->cd(2)->SetGridx();
1605 c1->cd(2)->SetGridy();
1606 ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,1))->Draw();
1607 c1->cd(3)->SetGridx();
1608 c1->cd(3)->SetGridy();
1609 ((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(i,2))->Draw();
1611 chargeVsDriftlengthAllShortPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,0));
1612 chargeVsDriftlengthAllMediumPads->Add((TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,1));
1613 chargeVsDriftlengthAllLongPads->Add( (TProfile*)fcalPadRegionChargeVsDriftlength->GetObject(0,2));
1616 c1->cd(1)->SetGridx();
1617 c1->cd(1)->SetGridy();
1618 chargeVsDriftlengthAllShortPads->Draw();
1619 c1->cd(2)->SetGridx();
1620 c1->cd(2)->SetGridy();
1621 chargeVsDriftlengthAllMediumPads->Draw();
1622 c1->cd(3)->SetGridx();
1623 c1->cd(3)->SetGridy();
1624 chargeVsDriftlengthAllLongPads->Draw();
1625 c1->Update(); // valgrind
1630 gSystem->ChangeDirectory("..");
1635 void AliTPCcalibTracks::FitResolutionNew(const char* pathName){
1637 // calculates different resulution fits in Y and Z direction
1638 // the histograms are written to 'ResolutionYZ.ps'
1639 // writes calculated resolution to 'resol.txt'
1640 // all files are stored in the directory pathName
1644 gSystem->MakeDirectory(pathName);
1645 gSystem->ChangeDirectory(pathName);
1649 if (GetDebugLevel() > -1) cout << "creating ResolutionYZ.ps..." << endl;
1650 TPostScript *ps = new TPostScript("ResolutionYZ.ps", 112);
1651 TF2 *fres = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1652 fres->SetParameter(0,0.02);
1653 fres->SetParameter(1,0.0054);
1654 fres->SetParameter(2,0.13);
1656 TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory
1658 // create histogramw for Y-resolution
1659 TH3F * hisResY0 = (TH3F*)fResolY->At(0);
1660 hisResY0->FitSlicesZ();
1661 TH2D * hisResY02 = (TH2D*)gDirectory->Get("Resol Y0_2");
1662 TH3F * hisResY1 = (TH3F*)fResolY->At(1);
1663 hisResY1->FitSlicesZ();
1664 TH2D * hisResY12 = (TH2D*)gDirectory->Get("Resol Y1_2");
1665 TH3F * hisResY2 = (TH3F*)fResolY->At(2);
1666 hisResY2->FitSlicesZ();
1667 TH2D * hisResY22 = (TH2D*)gDirectory->Get("Resol Y2_2");
1671 hisResY02->Fit(fres, "q"); // valgrind 132,072 bytes in 6 blocks are indirectly lost
1672 hisResY02->Draw("surf1");
1674 MakeDiff(hisResY02,fres)->Draw("surf1");
1676 // c.SaveAs("ResolutionYPad0.eps");
1679 hisResY12->Fit(fres, "q");
1680 hisResY12->Draw("surf1");
1682 MakeDiff(hisResY12,fres)->Draw("surf1");
1684 // c.SaveAs("ResolutionYPad1.eps");
1687 hisResY22->Fit(fres, "q");
1688 hisResY22->Draw("surf1");
1690 MakeDiff(hisResY22,fres)->Draw("surf1");
1692 // c.SaveAs("ResolutionYPad2.eps");
1694 // create histogramw for Z-resolution
1695 TH3F * hisResZ0 = (TH3F*)fResolZ->At(0);
1696 hisResZ0->FitSlicesZ();
1697 TH2D * hisResZ02 = (TH2D*)gDirectory->Get("Resol Z0_2");
1698 TH3F * hisResZ1 = (TH3F*)fResolZ->At(1);
1699 hisResZ1->FitSlicesZ();
1700 TH2D * hisResZ12 = (TH2D*)gDirectory->Get("Resol Z1_2");
1701 TH3F * hisResZ2 = (TH3F*)fResolZ->At(2);
1702 hisResZ2->FitSlicesZ();
1703 TH2D * hisResZ22 = (TH2D*)gDirectory->Get("Resol Z2_2");
1707 hisResZ02->Fit(fres, "q");
1708 hisResZ02->Draw("surf1");
1710 MakeDiff(hisResZ02,fres)->Draw("surf1");
1712 // c.SaveAs("ResolutionZPad0.eps");
1715 hisResZ12->Fit(fres, "q");
1716 hisResZ12->Draw("surf1");
1718 MakeDiff(hisResZ12,fres)->Draw("surf1");
1720 // c.SaveAs("ResolutionZPad1.eps");
1723 hisResZ22->Fit(fres, "q");
1724 hisResZ22->Draw("surf1");
1726 MakeDiff(hisResZ22,fres)->Draw("surf1");
1728 // c.SaveAs("ResolutionZPad2.eps");
1732 // write calculated resoltuions to 'resol.txt'
1733 ofstream fresol("resol.txt");
1734 fresol<<"Pad 0.75 cm"<<"\n";
1735 hisResY02->Fit(fres, "q"); // valgrind
1736 fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1737 hisResZ02->Fit(fres, "q");
1738 fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1740 fresol<<"Pad 1.00 cm"<<1<<"\n";
1741 hisResY12->Fit(fres, "q"); // valgrind
1742 fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1743 hisResZ12->Fit(fres, "q");
1744 fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1746 fresol<<"Pad 1.50 cm"<<0<<"\n";
1747 hisResY22->Fit(fres, "q");
1748 fresol<<"Y\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1749 hisResZ22->Fit(fres, "q");
1750 fresol<<"Z\t"<<fres->GetParameter(0)<<"\t"<<fres->GetParameter(1)<<"\t"<<fres->GetParameter(2)<<"\n";
1752 TH1::AddDirectory(kFALSE);
1753 gSystem->ChangeDirectory("..");
1758 void AliTPCcalibTracks::FitRMSNew(const char* pathName){
1760 // calculates different resulution-rms fits in Y and Z direction
1761 // the histograms are written to 'RMS_YZ.ps'
1762 // writes calculated resolution rms to 'rms.txt'
1763 // all files are stored in the directory pathName
1767 gSystem->MakeDirectory(pathName);
1768 gSystem->ChangeDirectory(pathName);
1770 TCanvas c; // valgrind 3 42,120 bytes in 405 blocks are still reachable 23,816 bytes in 229 blocks are still reachable
1772 if (GetDebugLevel() > -1) cout << "creating RMS_YZ.ps..." << endl;
1773 TPostScript *ps = new TPostScript("RMS_YZ.ps", 112);
1774 TF2 *frms = new TF2("fres","TMath::Sqrt([0]*[0]+[1]*[1]*x+[2]*[2]*y*y)",0,250,0,1);
1775 frms->SetParameter(0,0.02);
1776 frms->SetParameter(1,0.0054);
1777 frms->SetParameter(2,0.13);
1779 TH1::AddDirectory(kTRUE); // TH3F::FitSlicesZ() writes histograms into the current directory
1781 // create histogramw for Y-RMS
1782 TH3F * hisResY0 = (TH3F*)fRMSY->At(0);
1783 hisResY0->FitSlicesZ();
1784 TH2D * hisResY02 = (TH2D*)gDirectory->Get("RMS Y0_1");
1785 TH3F * hisResY1 = (TH3F*)fRMSY->At(1);
1786 hisResY1->FitSlicesZ();
1787 TH2D * hisResY12 = (TH2D*)gDirectory->Get("RMS Y1_1");
1788 TH3F * hisResY2 = (TH3F*)fRMSY->At(2);
1789 hisResY2->FitSlicesZ();
1790 TH2D * hisResY22 = (TH2D*)gDirectory->Get("RMS Y2_1");
1794 hisResY02->Fit(frms, "qn0");
1795 hisResY02->Draw("surf1");
1797 MakeDiff(hisResY02,frms)->Draw("surf1");
1799 // c.SaveAs("RMSYPad0.eps");
1802 hisResY12->Fit(frms, "qn0"); // valgrind several blocks possibly lost
1803 hisResY12->Draw("surf1");
1805 MakeDiff(hisResY12,frms)->Draw("surf1");
1807 // c.SaveAs("RMSYPad1.eps");
1810 hisResY22->Fit(frms, "qn0");
1811 hisResY22->Draw("surf1");
1813 MakeDiff(hisResY22,frms)->Draw("surf1");
1815 // c.SaveAs("RMSYPad2.eps");
1817 // create histogramw for Z-RMS
1818 TH3F * hisResZ0 = (TH3F*)fRMSZ->At(0);
1819 hisResZ0->FitSlicesZ();
1820 TH2D * hisResZ02 = (TH2D*)gDirectory->Get("RMS Z0_1");
1821 TH3F * hisResZ1 = (TH3F*)fRMSZ->At(1);
1822 hisResZ1->FitSlicesZ();
1823 TH2D * hisResZ12 = (TH2D*)gDirectory->Get("RMS Z1_1");
1824 TH3F * hisResZ2 = (TH3F*)fRMSZ->At(2);
1825 hisResZ2->FitSlicesZ();
1826 TH2D * hisResZ22 = (TH2D*)gDirectory->Get("RMS Z2_1");
1830 hisResZ02->Fit(frms, "qn0"); // valgrind
1831 hisResZ02->Draw("surf1");
1833 MakeDiff(hisResZ02,frms)->Draw("surf1");
1835 // c.SaveAs("RMSZPad0.eps");
1838 hisResZ12->Fit(frms, "qn0");
1839 hisResZ12->Draw("surf1");
1841 MakeDiff(hisResZ12,frms)->Draw("surf1");
1843 // c.SaveAs("RMSZPad1.eps");
1846 hisResZ22->Fit(frms, "qn0"); // valgrind 1 block possibly lost
1847 hisResZ22->Draw("surf1");
1849 MakeDiff(hisResZ22,frms)->Draw("surf1");
1851 // c.SaveAs("RMSZPad2.eps");
1853 // write calculated resoltuion rms to 'rms.txt'
1854 ofstream filerms("rms.txt");
1855 filerms<<"Pad 0.75 cm"<<"\n";
1856 hisResY02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost
1857 filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1858 hisResZ02->Fit(frms, "qn0"); // valgrind 23 blocks indirectly lost
1859 filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1861 filerms<<"Pad 1.00 cm"<<1<<"\n";
1862 hisResY12->Fit(frms, "qn0"); // valgrind 3,256 bytes in 22 blocks are indirectly lost
1863 filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1864 hisResZ12->Fit(frms, "qn0"); // valgrind 66,036 bytes in 3 blocks are still reachable
1865 filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1867 filerms<<"Pad 1.50 cm"<<0<<"\n";
1868 hisResY22->Fit(frms, "qn0"); // valgrind 40,139 bytes in 11 blocks are still reachable 330,180 bytes in 15 blocks are possibly lost
1869 filerms<<"Y\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1870 hisResZ22->Fit(frms, "qn0");
1871 filerms<<"Z\t"<<frms->GetParameter(0)<<"\t"<<frms->GetParameter(1)<<"\t"<<frms->GetParameter(2)<<"\n";
1873 TH1::AddDirectory(kFALSE);
1874 gSystem->ChangeDirectory("..");
1881 void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
1883 // Make tree with resolution parameters
1884 // the result is written to 'resol.root' in directory 'pathname'
1885 // file information are available in fileInfo
1886 // available variables in the tree 'Resol':
1887 // Entries: number of entries for this resolution point
1888 // nbins: number of bins that were accumulated
1889 // Dim: direction, Dim==0: y-direction, Dim==1: z-direction
1890 // Pad: padSize; short, medium and long
1891 // Length: pad length, 0.75, 1, 1.5
1892 // QMean: mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
1893 // Zc: center of middle bin in drift direction
1894 // Zm: mean dirftlength for accumulated Delta-Histograms
1895 // Zs: width of driftlength bin
1896 // AngleC: center of middle bin in Angle-Direction
1897 // AngleM: mean angle for accumulated Delta-Histograms
1898 // AngleS: width of Angle-bin
1899 // Resol: sigma for gaus fit through Delta-Histograms
1900 // Sigma: error of sigma for gaus fit through Delta Histograms
1901 // MeanR: mean of the Delta-Histogram
1902 // SigmaR: rms of the Delta-Histogram
1903 // RMSm: mean of the gaus fit through RMS-Histogram
1904 // RMS: sigma of the gaus fit through RMS-Histogram
1905 // RMSe0: error of mean of gaus fit in RMS-Histogram
1906 // RMSe1: error of sigma of gaus fit in RMS-Histogram
1909 if (GetDebugLevel() > -1) cout << " ***** this is MakeResPlotsQTree *****" << endl;
1910 if (GetDebugLevel() > -1) cout << " relax, the calculation will take a while..." << endl;
1912 gSystem->MakeDirectory(pathName);
1913 gSystem->ChangeDirectory(pathName);
1914 TString kFileName = "resol.root";
1915 TTreeSRedirector fTreeResol(kFileName.Data());
1917 TH3F *resArray[2][3][11];
1918 TH3F *rmsArray[2][3][11];
1920 // load histograms from fArrayQDY and fArrayQDZ
1921 // into resArray and rmsArray
1922 // that is all we need here
1923 for (Int_t idim = 0; idim < 2; idim++){
1924 for (Int_t ipad = 0; ipad < 3; ipad++){
1925 for (Int_t iq = 0; iq <= 10; iq++){
1926 rmsArray[idim][ipad][iq]=0;
1927 resArray[idim][ipad][iq]=0;
1928 Int_t bin = GetBin(iq,ipad);
1930 if (idim == 0) hresl = (TH3F*)fArrayQDY->At(bin);
1931 if (idim == 1) hresl = (TH3F*)fArrayQDZ->At(bin);
1932 if (!hresl) continue;
1933 resArray[idim][ipad][iq] = (TH3F*) hresl->Clone();
1934 resArray[idim][ipad][iq]->SetDirectory(0);
1935 TH3F * hreslRMS = 0;
1936 if (idim == 0) hreslRMS = (TH3F*)fArrayQRMSY->At(bin);
1937 if (idim == 1) hreslRMS = (TH3F*)fArrayQRMSZ->At(bin);
1938 if (!hreslRMS) continue;
1939 rmsArray[idim][ipad][iq] = (TH3F*) hreslRMS->Clone();
1940 rmsArray[idim][ipad][iq]->SetDirectory(0);
1944 if (GetDebugLevel() > -1) cout << "Histograms loaded, starting to proces..." << endl;
1946 //--------------------------------------------------------------------------------------------
1950 Double_t zMean, angleMean, zCenter, angleCenter;
1951 Double_t zSigma, angleSigma;
1952 TH1D *projectionRes = new TH1D("projectionRes", "projectionRes", 50, -1, 1);
1953 TH1D *projectionRms = new TH1D("projectionRms", "projectionRms", 50, -1, 1);
1954 TF1 *fitFunction = new TF1("fitFunction", "gaus");
1955 Float_t entriesQ = 0;
1956 Int_t loopCounter = 1;
1958 for (Int_t idim = 0; idim < 2; idim++){
1959 // Loop y-z corrdinate
1960 for (Int_t ipad = 0; ipad < 3; ipad++){
1962 for (Int_t iq = -1; iq < 10; iq++){
1964 if (GetDebugLevel() > -1)
1965 cout << "Loop-counter, this is loop " << loopCounter << " of 66, ("
1966 << (Int_t)((loopCounter)/66.*100) << "% done), "
1967 << "idim = " << idim << ", ipad = " << ipad << ", iq = " << iq << " \r" << std::flush;
1976 // integrated spectra
1977 for (Int_t iql = 0; iql < 10; iql++){
1978 Int_t bin = GetBin(iql,ipad);
1979 TH3F *hresl = resArray[idim][ipad][iql];
1980 TH3F *hrmsl = rmsArray[idim][ipad][iql];
1981 if (!hresl) continue;
1982 if (!hrmsl) continue;
1983 entriesQ += hresl->GetEntries();
1984 qMean += hresl->GetEntries() * GetQ(bin);
1986 hres = (TH3F*)hresl->Clone();
1987 hrms = (TH3F*)hrmsl->Clone();
1995 qMean *= -1.; // integral mean charge
1998 // loop over neighboured Q-bins
1999 // accumulate entries from neighboured Q-bins
2000 for (Int_t iql = iq - 1; iql <= iq + 1; iql++){
2001 if (iql < 0) continue;
2002 Int_t bin = GetBin(iql,ipad);
2003 TH3F * hresl = resArray[idim][ipad][iql];
2004 TH3F * hrmsl = rmsArray[idim][ipad][iql];
2005 if (!hresl) continue;
2006 if (!hrmsl) continue;
2007 entriesQ += hresl->GetEntries();
2008 qMean += hresl->GetEntries() * GetQ(bin);
2010 hres = (TH3F*) hresl->Clone();
2011 hrms = (TH3F*) hrmsl->Clone();
2020 if (!hres) continue;
2021 if (!hrms) continue;
2023 TAxis *xAxisDriftLength = hres->GetXaxis(); // driftlength / z - axis
2024 TAxis *yAxisAngle = hres->GetYaxis(); // angle axis
2025 TAxis *zAxisDelta = hres->GetZaxis(); // delta axis
2026 TAxis *zAxisRms = hrms->GetZaxis(); // rms axis
2028 // loop over all angle bins
2029 for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++) {
2030 angleCenter = yAxisAngle->GetBinCenter(ibinyAngle);
2031 // loop over all driftlength bins
2032 for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++) {
2033 zCenter = xAxisDriftLength->GetBinCenter(ibinxDL);
2034 zSigma = xAxisDriftLength->GetBinWidth(ibinxDL);
2035 angleSigma = yAxisAngle->GetBinWidth(ibinyAngle);
2036 zMean = zCenter; // changens, when more statistic is accumulated
2037 angleMean = angleCenter; // changens, when more statistic is accumulated
2039 // create 2 1D-Histograms, projectionRes and projectionRms
2040 // these histograms are delta histograms for given direction, padSize, chargeBin,
2041 // angleBin and driftLengthBin
2042 // later on they will be fitted with a gausian, its sigma is the resoltuion...
2043 sprintf(name,"%s, zCenter: %f, angleCenter: %f", hres->GetName(), zCenter, angleCenter);
2044 // TH1D * projectionRes = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
2045 projectionRes->SetNameTitle(name, name);
2046 sprintf(name,"%s, zCenter: %f, angleCenter: %f", hrms->GetName(),zCenter,angleCenter);
2047 // TH1D * projectionRms = new TH1D(name, name, zAxisDelta->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
2048 projectionRms->SetNameTitle(name, name);
2050 projectionRes->Reset();
2051 projectionRes->SetBins(zAxisDelta->GetNbins(), zAxisDelta->GetXmin(), zAxisDelta->GetXmax());
2052 projectionRms->Reset();
2053 projectionRms->SetBins(zAxisRms->GetNbins(), zAxisRms->GetXmin(), zAxisRms->GetXmax());
2054 projectionRes->SetDirectory(0);
2055 projectionRms->SetDirectory(0);
2057 Double_t entries = 0;
2058 Int_t nbins = 0; // counts, how many bins were accumulated
2063 // fill projectionRes and projectionRms for given dim, ipad and iq,
2064 // as well as for given angleBin and driftlengthBin
2065 // if this gives not enough statistic, include neighbourhood
2066 // (angle and driftlength) successifely
2067 for (Int_t dbin = 0; dbin <= 8; dbin++){ // delta-bins around centered angleBin and driftlengthBin
2068 for (Int_t dbiny2 = -1; dbiny2 <= 1; dbiny2++) { // delta-bins in angle direction
2069 for (Int_t dbinx2 = -3; dbinx2 <= 3; dbinx2++){ // delta-bins in driftlength direction
2070 if (TMath::Abs(dbinx2) + TMath::Abs(dbiny2) != dbin) continue; // add each bin only one time !
2071 Int_t binx2 = ibinxDL + dbinx2; // position variable in x (driftlength) direction
2072 Int_t biny2 = ibinyAngle + dbiny2; // position variable in y (angle) direction
2073 if (binx2 < 1 || biny2 < 1) continue; // don't go out of the histogram!
2074 if (binx2 >= xAxisDriftLength->GetNbins()) continue; // don't go out of the histogram!
2075 if (biny2 >= yAxisAngle->GetNbins()) continue; // don't go out of the histogram!
2076 nbins++; // count the number of accumulated bins
2077 // Fill resolution histo
2078 for (Int_t ibin3 = 1; ibin3 < zAxisDelta->GetNbins(); ibin3++) {
2079 // Int_t content = (Int_t)hres->GetBinContent(binx2, biny2, ibin3); // unused variable
2080 projectionRes->Fill(zAxisDelta->GetBinCenter(ibin3), hres->GetBinContent(binx2, biny2, ibin3));
2081 entries += hres->GetBinContent(binx2, biny2, ibin3);
2082 zMean += hres->GetBinContent(binx2, biny2, ibin3) * xAxisDriftLength->GetBinCenter(binx2);
2083 angleMean += hres->GetBinContent(binx2, biny2, ibin3) * yAxisAngle->GetBinCenter(biny2);
2086 for (Int_t ibin3 = 1; ibin3 < zAxisRms->GetNbins(); ibin3++) {
2087 projectionRms->Fill(zAxisRms->GetBinCenter(ibin3), hrms->GetBinContent(binx2, biny2, ibin3));
2090 if (entries > minEntries) break; // enough statistic accumulated
2092 if (entries > minEntries) break; // enough statistic accumulated
2094 if ( entries< minEntries) continue; // when it was absolutly impossible to get enough statistic, don't write this point into the resolution tree
2096 angleMean /= entries;
2098 if (entries > minEntries) {
2099 // when enough statistic is accumulated
2100 // fit Delta histograms with a gausian
2101 // of the gausian is the resolution (resol), its fit error is sigma
2102 // store also mean and RMS of the histogram
2103 Float_t xmin = projectionRes->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2104 Float_t xmax = projectionRes->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2106 // projectionRes->Fit("gaus", "q0", "", xmin, xmax);
2107 // Float_t resol = projectionRes->GetFunction("gaus")->GetParameter(2);
2108 // Float_t sigma = projectionRes->GetFunction("gaus")->GetParError(2);
2109 fitFunction->SetMaximum(xmax);
2110 fitFunction->SetMinimum(xmin);
2111 projectionRes->Fit("fitFunction", "qN0", "", xmin, xmax);
2112 Float_t resol = fitFunction->GetParameter(2);
2113 Float_t sigma = fitFunction->GetParError(2);
2115 Float_t meanR = projectionRes->GetMean();
2116 Float_t sigmaR = projectionRes->GetRMS();
2117 // fit also RMS histograms with a gausian
2118 // store mean and sigma of the gausian in rmsMean and rmsSigma
2119 // store also the fit errors in errorRMS and errorSigma
2120 xmin = projectionRms->GetMean() - 2. * projectionRes->GetRMS() - 0.2;
2121 xmax = projectionRms->GetMean() + 2. * projectionRes->GetRMS() + 0.2;
2123 // projectionRms->Fit("gaus","q0","",xmin,xmax);
2124 // Float_t rmsMean = projectionRms->GetFunction("gaus")->GetParameter(1);
2125 // Float_t rmsSigma = projectionRms->GetFunction("gaus")->GetParameter(2);
2126 // Float_t errorRMS = projectionRms->GetFunction("gaus")->GetParError(1);
2127 // Float_t errorSigma = projectionRms->GetFunction("gaus")->GetParError(2);
2128 projectionRms->Fit("fitFunction", "qN0", "", xmin, xmax);
2129 Float_t rmsMean = fitFunction->GetParameter(1);
2130 Float_t rmsSigma = fitFunction->GetParameter(2);
2131 Float_t errorRMS = fitFunction->GetParError(1);
2132 Float_t errorSigma = fitFunction->GetParError(2);
2134 Float_t length = 0.75;
2135 if (ipad == 1) length = 1;
2136 if (ipad == 2) length = 1.5;
2138 fTreeResol<<"Resol"<<
2139 "Entries="<<entries<< // number of entries for this resolution point
2140 "nbins="<<nbins<< // number of bins that were accumulated
2141 "Dim="<<idim<< // direction, Dim==0: y-direction, Dim==1: z-direction
2142 "Pad="<<ipad<< // padSize; short, medium and long
2143 "Length="<<length<< // pad length, 0.75, 1, 1.5
2144 "QMean="<<qMean<< // mean charge of current charge bin and its neighbours, Qmean<0: integrated spectra
2145 "Zc="<<zCenter<< // center of middle bin in drift direction
2146 "Zm="<<zMean<< // mean dirftlength for accumulated Delta-Histograms
2147 "Zs="<<zSigma<< // width of driftlength bin
2148 "AngleC="<<angleCenter<< // center of middle bin in Angle-Direction
2149 "AngleM="<<angleMean<< // mean angle for accumulated Delta-Histograms
2150 "AngleS="<<angleSigma<< // width of Angle-bin
2151 "Resol="<<resol<< // sigma for gaus fit through Delta-Histograms
2152 "Sigma="<<sigma<< // error of sigma for gaus fit through Delta Histograms
2153 "MeanR="<<meanR<< // mean of the Delta-Histogram
2154 "SigmaR="<<sigmaR<< // rms of the Delta-Histogram
2155 "RMSm="<<rmsMean<< // mean of the gaus fit through RMS-Histogram
2156 "RMSs="<<rmsSigma<< // sigma of the gaus fit through RMS-Histogram
2157 "RMSe0="<<errorRMS<< // error of mean of gaus fit in RMS-Histogram
2158 "RMSe1="<<errorSigma<< // error of sigma of gaus fit in RMS-Histogram
2160 if (GetDebugLevel() > 5) {
2161 projectionRes->SetDirectory(fTreeResol.GetFile());
2162 projectionRes->Write(projectionRes->GetName());
2163 projectionRes->SetDirectory(0);
2164 projectionRms->SetDirectory(fTreeResol.GetFile());
2165 projectionRms->Write(projectionRms->GetName());
2166 projectionRes->SetDirectory(0);
2168 } // if (projectionRes->GetSum() > minEntries)
2169 } // for (Int_t ibinxDL = 1; ibinxDL <= xAxisDriftLength->GetNbins(); ibinxDL++)
2170 } // for (Int_t ibinyAngle = 1; ibinyAngle <= yAxisAngle->GetNbins(); ibinyAngle++)
2175 delete projectionRes;
2176 delete projectionRms;
2178 // TFile resolFile(fTreeResol.GetFile());
2179 TObjString fileInfo(Form("Resolution tree, minEntries = %i", minEntries));
2180 fileInfo.Write("fileInfo");
2181 // resolFile.Close();
2182 // fTreeResol.GetFile()->Close();
2183 if (GetDebugLevel() > -1) cout << endl;
2184 if (GetDebugLevel() > -1) cout << "MakeResPlotsQTree done, results are in '"<< kFileName.Data() <<"'." << endl;
2185 gSystem->ChangeDirectory("..");
2192 Long64_t AliTPCcalibTracks::Merge(TCollection *collectionList) {
2194 // function to merge several AliTPCcalibTracks objects after PROOF calculation
2195 // The object's histograms are merged via their merge functions
2196 // Be carefull: histograms are linked to a file, switch this off by TH1::AddDirectory(kFALSE) !!!
2199 if (GetDebugLevel() > 0) cout << " ***** this is AliTPCcalibTracks::Merge(TCollection *collectionList) *****"<< endl;
2200 if (!collectionList) return 0;
2201 if (collectionList->IsEmpty()) return -1;
2203 if (GetDebugLevel() > 1) cout << "the collectionList contains " << collectionList->GetEntries() << " entries." << endl; // REMOVE THIS LINE!!!!!!!!!!!!!!!!!1
2204 if (GetDebugLevel() > 5) cout << " the list in the merge-function looks as follows: " << endl;
2205 collectionList->Print();
2207 // create a list for each data member
2208 TList* deltaYList = new TList;
2209 TList* deltaZList = new TList;
2210 TList* arrayAmpRowList = new TList;
2211 TList* rejectedTracksList = new TList;
2212 TList* hclusList = new TList;
2213 TList* clusterCutHistoList = new TList;
2214 TList* arrayAmpList = new TList;
2215 TList* arrayQDYList = new TList;
2216 TList* arrayQDZList = new TList;
2217 TList* arrayQRMSYList = new TList;
2218 TList* arrayQRMSZList = new TList;
2219 TList* arrayChargeVsDriftlengthList = new TList;
2220 TList* calPadRegionChargeVsDriftlengthList = new TList;
2221 TList* hclusterPerPadrowList = new TList;
2222 TList* hclusterPerPadrowRawList = new TList;
2223 TList* resolYList = new TList;
2224 TList* resolZList = new TList;
2225 TList* rMSYList = new TList;
2226 TList* rMSZList = new TList;
2228 // TList* nRowsList = new TList;
2229 // TList* nSectList = new TList;
2230 // TList* fileNoList = new TList;
2232 TIterator *listIterator = collectionList->MakeIterator();
2233 AliTPCcalibTracks *calibTracks = 0;
2234 if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;
2236 while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
2237 // loop over all entries in the collectionList and get dataMembers into lists
2239 deltaYList->Add( calibTracks->GetfDeltaY() );
2240 deltaZList->Add( calibTracks->GetfDeltaZ() );
2241 arrayAmpRowList->Add(calibTracks->GetfArrayAmpRow());
2242 arrayAmpList->Add(calibTracks->GetfArrayAmp());
2243 arrayQDYList->Add(calibTracks->GetfArrayQDY());
2244 arrayQDZList->Add(calibTracks->GetfArrayQDZ());
2245 arrayQRMSYList->Add(calibTracks->GetfArrayQRMSY());
2246 arrayQRMSZList->Add(calibTracks->GetfArrayQRMSZ());
2247 resolYList->Add(calibTracks->GetfResolY());
2248 resolZList->Add(calibTracks->GetfResolZ());
2249 rMSYList->Add(calibTracks->GetfRMSY());
2250 rMSZList->Add(calibTracks->GetfRMSZ());
2251 arrayChargeVsDriftlengthList->Add(calibTracks->GetfArrayChargeVsDriftlength());
2252 calPadRegionChargeVsDriftlengthList->Add(calibTracks->GetCalPadRegionchargeVsDriftlength());
2253 hclusList->Add(calibTracks->GetfHclus());
2254 rejectedTracksList->Add(calibTracks->GetfRejectedTracksHisto());
2255 clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
2256 hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow());
2257 hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw());
2259 if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
2260 fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
2261 // fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
2263 if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
2264 AddHistos(calibTracks);
2268 // merge data members
2269 if (GetDebugLevel() > 0) cout << "histogram's merge-functins are called... " << endl;
2270 fDeltaY->Merge(deltaYList);
2271 fDeltaZ->Merge(deltaZList);
2272 fHclus->Merge(hclusList);
2273 fClusterCutHisto->Merge(clusterCutHistoList);
2274 fRejectedTracksHisto->Merge(rejectedTracksList);
2275 fHclusterPerPadrow->Merge(hclusterPerPadrowList);
2276 fHclusterPerPadrowRaw->Merge(hclusterPerPadrowRawList);
2278 TObjArray* objarray = 0;
2280 TList* histList = 0;
2281 TIterator *objListIterator = 0;
2283 if (GetDebugLevel() > 0) cout << "merging fArrayAmpRows..." << endl;
2284 // merge fArrayAmpRows
2285 for (Int_t i = 0; i < fArrayAmpRow->GetEntriesFast(); i++ ) { // loop over data member, i<72
2286 objListIterator = arrayAmpRowList->MakeIterator();
2287 histList = new TList;
2288 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2289 // loop over arrayAmpRowList, get TObjArray, get object at position i, cast it into TProfile
2290 hist = (TProfile*)(objarray->At(i));
2291 histList->Add(hist);
2293 ((TProfile*)(fArrayAmpRow->At(i)))->Merge(histList);
2295 delete objListIterator;
2298 if (GetDebugLevel() > 0) cout << "merging fArrayAmps..." << endl;
2300 for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) { // loop over data member, i<72
2301 TIterator *cobjListIterator = arrayAmpList->MakeIterator();
2302 histList = new TList;
2303 while (( objarray = (TObjArray*)cobjListIterator->Next() )) {
2304 // loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F
2305 hist = (TH1F*)(objarray->At(i));
2306 histList->Add(hist);
2308 ((TH1F*)(fArrayAmp->At(i)))->Merge(histList);
2310 delete cobjListIterator;
2313 if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
2315 for (Int_t i = 0; i < fArrayQDY->GetEntriesFast(); i++) { // loop over data member, i < 300
2316 objListIterator = arrayQDYList->MakeIterator();
2317 histList = new TList;
2318 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2319 // loop over arrayQDYList, get TObjArray, get object at position i, cast it into TH3F
2320 hist = (TH3F*)(objarray->At(i));
2321 histList->Add(hist);
2323 ((TH3F*)(fArrayQDY->At(i)))->Merge(histList);
2325 delete objListIterator;
2328 if (GetDebugLevel() > 0) cout << "merging fArrayQDZ..." << endl;
2330 for (Int_t i = 0; i < fArrayQDZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2331 objListIterator = arrayQDZList->MakeIterator();
2332 histList = new TList;
2333 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2334 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2335 if (!objarray) continue;
2336 hist = (TH3F*)(objarray->At(i));
2337 histList->Add(hist);
2339 ((TH3F*)(fArrayQDZ->At(i)))->Merge(histList);
2341 delete objListIterator;
2344 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSY..." << endl;
2345 // merge fArrayQRMSY
2346 for (Int_t i = 0; i < fArrayQRMSY->GetEntriesFast(); i++) { // loop over data member, i < 300
2347 objListIterator = arrayQRMSYList->MakeIterator();
2348 histList = new TList;
2349 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2350 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2351 if (!objarray) continue;
2352 hist = (TH3F*)(objarray->At(i));
2353 histList->Add(hist);
2355 ((TH3F*)(fArrayQRMSY->At(i)))->Merge(histList);
2357 delete objListIterator;
2360 if (GetDebugLevel() > 0) cout << "merging fArrayQRMSZ..." << endl;
2361 // merge fArrayQRMSZ
2362 for (Int_t i = 0; i < fArrayQRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 300
2363 objListIterator = arrayQRMSZList->MakeIterator();
2364 histList = new TList;
2365 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2366 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2367 if (!objarray) continue;
2368 hist = (TH3F*)(objarray->At(i));
2369 histList->Add(hist);
2371 ((TH3F*)(fArrayQRMSZ->At(i)))->Merge(histList);
2373 delete objListIterator;
2376 if (GetDebugLevel() > 0) cout << "merging fArrayChargeVsDriftlength..." << endl;
2377 // merge fArrayChargeVsDriftlength
2378 for (Int_t i = 0; i < fArrayChargeVsDriftlength->GetEntriesFast(); i++) { // loop over data member, i < 300
2379 objListIterator = arrayChargeVsDriftlengthList->MakeIterator();
2380 histList = new TList;
2381 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2382 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TProfile
2383 if (!objarray) continue;
2384 hist = (TProfile*)(objarray->At(i));
2385 histList->Add(hist);
2387 ((TProfile*)(fArrayChargeVsDriftlength->At(i)))->Merge(histList);
2389 delete objListIterator;
2392 if (GetDebugLevel() > 0) cout << "merging fcalPadRegionChargeVsDriftlength..." << endl;
2393 // merge fcalPadRegionChargeVsDriftlength
2394 AliTPCCalPadRegion *cpr = 0x0;
2397 TIterator *regionIterator = fcalPadRegionChargeVsDriftlength->MakeIterator();
2398 while (hist = (TProfile*)regionIterator->Next()) {
2399 // loop over all calPadRegion's in destination calibTracks object
2400 objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2401 while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) {
2408 for (UInt_t isec = 0; isec < 36; isec++) {
2409 for (UInt_t padSize = 0; padSize < 3; padSize++){
2410 objListIterator = calPadRegionChargeVsDriftlengthList->MakeIterator();
2411 histList = new TList;
2412 while (( cpr = (AliTPCCalPadRegion*)objListIterator->Next() )) {
2413 // loop over calPadRegionChargeVsDriftlengthList, get AliTPCCalPadRegion, get object
2414 hist = (TProfile*)cpr->GetObject(isec, padSize);
2415 histList->Add(hist);
2417 ((TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isec, padSize)))->Merge(histList);
2419 delete objListIterator;
2426 if (GetDebugLevel() > 0) cout << "starting to merge the rest: fResolY, fResolZ , fRMSY, fRMSZ..." << endl;
2428 for (Int_t i = 0; i < fResolY->GetEntriesFast(); i++) { // loop over data member, i < 3
2429 objListIterator = resolYList->MakeIterator();
2430 histList = new TList;
2431 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2432 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2433 if (!objarray) continue;
2434 hist = (TH3F*)(objarray->At(i));
2435 histList->Add(hist);
2437 ((TH3F*)(fResolY->At(i)))->Merge(histList);
2439 delete objListIterator;
2443 for (Int_t i = 0; i < fResolZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2444 objListIterator = resolZList->MakeIterator();
2445 histList = new TList;
2446 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2447 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2448 if (!objarray) continue;
2449 hist = (TH3F*)(objarray->At(i));
2450 histList->Add(hist);
2452 ((TH3F*)(fResolZ->At(i)))->Merge(histList);
2454 delete objListIterator;
2458 for (Int_t i = 0; i < fRMSY->GetEntriesFast(); i++) { // loop over data member, i < 3
2459 objListIterator = rMSYList->MakeIterator();
2460 histList = new TList;
2461 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2462 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2463 if (!objarray) continue;
2464 hist = (TH3F*)(objarray->At(i));
2465 histList->Add(hist);
2467 ((TH3F*)(fRMSY->At(i)))->Merge(histList);
2469 delete objListIterator;
2473 for (Int_t i = 0; i < fRMSZ->GetEntriesFast(); i++) { // loop over data member, i < 3
2474 objListIterator = rMSZList->MakeIterator();
2475 histList = new TList;
2476 while (( objarray = (TObjArray*)objListIterator->Next() )) {
2477 // loop over arrayQDZList, get TObjArray, get object at position i, cast it into TH3F
2478 if (!objarray) continue;
2479 hist = (TH3F*)(objarray->At(i));
2480 histList->Add(hist);
2482 ((TH3F*)(fRMSZ->At(i)))->Merge(histList);
2484 delete objListIterator;
2489 delete arrayAmpRowList;
2490 delete arrayAmpList;
2491 delete arrayQDYList;
2492 delete arrayQDZList;
2493 delete arrayQRMSYList;
2494 delete arrayQRMSZList;
2499 delete listIterator;
2501 if (GetDebugLevel() > 0) cout << "merging done!" << endl;
2507 AliTPCcalibTracks* AliTPCcalibTracks::TestMerge(AliTPCcalibTracks *ct, AliTPCClusterParam *clusterParam, Int_t nCalTracks){
2509 // function to test AliTPCcalibTrack::Merge:
2510 // in the file 'f' is a AliTPCcalibTrack object with name "calibTracks"
2511 // this object is appended 'nCalTracks' times to a TList
2512 // A new AliTPCcalibTrack object is created which merges the list
2513 // this object is returned
2516 // .L AliTPCcalibTracks.cxx+g
2518 TFile f("Output.root");
2519 AliTPCcalibTracks* calTracks = (AliTPCcalibTracks*)f.Get("calibTracks");
2521 TFile clusterParamFile("/u/lbozyk/calibration/workdir/calibTracks/TPCClusterParam.root");
2522 AliTPCClusterParam *clusterParam = (AliTPCClusterParam *) clusterParamFile.Get("Param");
2523 clusterParamFile.Close();
2525 AliTPCcalibTracks::TestMerge(calTracks, clusterParam);
2527 TList *list = new TList();
2528 if (ct == 0 || clusterParam == 0) return 0;
2529 cout << "making list with " << nCalTracks << " AliTPCcalibTrack objects" << endl;
2530 for (Int_t i = 0; i < nCalTracks; i++) {
2531 if (i%10==0) cout << "Adding element " << i << " of " << nCalTracks << endl;
2532 list->Add(new AliTPCcalibTracks(*ct));
2535 // only for check at the end
2536 AliTPCcalibTracks* cal1 = new AliTPCcalibTracks(*ct);
2537 Double_t cal1Entries = ((TH1F*)cal1->GetfArrayAmpRow()->At(5))->GetEntries();
2538 // Double_t cal1Entries = 5; //((TH1F*)ct->GetfArrayAmpRow()->At(5))->GetEntries();
2540 cout << "The list contains " << list->GetEntries() << " entries. " << endl;
2543 AliTPCcalibTracksCuts *cuts = new AliTPCcalibTracksCuts(20, 0.4, 0.5, 0.13, 0.018);
2544 AliTPCcalibTracks* cal = new AliTPCcalibTracks("calTracksMerged", "calTracksMerged", clusterParam, cuts, 5);
2547 cout << "cal->GetfArrayAmpRow()->At(5)->Print():" << endl;
2548 cal->GetfArrayAmpRow()->At(5)->Print();
2549 Double_t calEntries = ((TH1F*)cal->GetfArrayAmpRow()->At(5))->GetEntries();
2551 cout << "cal1->GetfArrayAmpRow()->At(5))->GetEntries() = " << cal1Entries << endl;
2552 cout << " cal->GetfArrayAmpRow()->At(5))->GetEntries() = " << calEntries << endl;
2553 printf("That means there were %f / %f = %f AliTPCcalibTracks-Objects merged. \n",
2554 calEntries, cal1Entries, ((Double_t)calEntries/cal1Entries));
2561 void AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints){
2563 // Make position corrections
2564 // for the moment Only using debug streamer
2565 // chainres - debug tree
2566 // param - parameters to be updated
2567 // maxPoints - maximal number of points using for fit
2568 // verbose - print info flag
2570 // Current implementation - only using debug streamers
2575 Int_t maxPoints=100000;
2580 gSystem->Load("libANALYSIS");
2581 gSystem->Load("libSTAT");
2582 gSystem->Load("libTPCcalib");
2585 //1. Load Parameters to be modified:
2588 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
2589 AliCDBManager::Instance()->SetRun(0)
2590 AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
2592 //2. Load chain from debug streamers
2595 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2596 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2597 AliXRDPROOFtoolkit tool;
2598 TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
2600 //3. Do fits and store results
2602 AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ;
2603 TFile f("paramout.root","recreate");
2604 param->Write("clusterParam");
2607 TFile f2("paramout.root");
2608 AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
2609 param2->SetInstance(param2);
2610 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
2615 TStatToolkit toolkit;
2617 TVectorD fitParamY0;
2618 TVectorD fitParamY1;
2619 TVectorD fitParamZ0;
2620 TVectorD fitParamZ1;
2624 chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
2625 chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
2626 chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
2627 chainres->SetAlias("st","(sin(dt)-dt)");
2629 chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
2634 TString fstringY="";
2636 fstringY+="(dp)++"; //1
2637 fstringY+="(dp)*di++"; //2
2638 fstringY+="(sp)++"; //3
2639 fstringY+="(sp)*di++"; //4
2640 TString fstringZ="";
2641 fstringZ+="(dt)++"; //1
2642 fstringZ+="(dt)*di++"; //2
2643 fstringZ+="(st)++"; //3
2644 fstringZ+="(st)*di++"; //4
2648 TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
2649 printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2650 param->PosZcor(0) = (TVectorD*) fitParamZ0.Clone();
2652 TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
2654 printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2655 param->PosZcor(1) = (TVectorD*) fitParamZ1.Clone();
2656 param->PosZcor(2) = (TVectorD*) fitParamZ1.Clone();
2660 TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
2661 printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2662 param->PosYcor(0) = (TVectorD*) fitParamY0.Clone();
2665 TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
2667 printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
2668 param->PosYcor(1) = (TVectorD*) fitParamY1.Clone();
2669 param->PosYcor(2) = (TVectorD*) fitParamY1.Clone();
2673 chainres->SetAlias("fitZ0",strZ0->Data());
2674 chainres->SetAlias("fitZ1",strZ1->Data());
2675 chainres->SetAlias("fitY0",strY0->Data());
2676 chainres->SetAlias("fitY1",strY1->Data());
2677 // chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000);
2682 void AliTPCcalibTracks::MakeHistos(){
2686 //THnSparse *fHisDeltaY; // THnSparse - delta Y
2687 //THnSparse *fHisDeltaZ; // THnSparse - delta Z
2688 //THnSparse *fHisRMSY; // THnSparse - rms Y
2689 //THnSparse *fHisRMSZ; // THnSparse - rms Z
2690 //THnSparse *fHisQmax; // THnSparse - qmax
2691 //THnSparse *fHisQtot; // THnSparse - qtot
2692 // cluster performance bins
2693 // 0 - variable of interest
2694 // 1 - pad type - 0- short 1-medium 2-long pads
2695 // 2 - drift length - drift length -0-1
2696 // 3 - Qmax - Qmax - 2- 400
2697 // 4 - cog - COG position - 0-1
2698 // 5 - tan(phi) - local y angle
2699 // 6 - tan(theta) - local z angle
2700 // 7 - sector - sector number
2701 Double_t xminTrack[8], xmaxTrack[8];
2703 TString axisName[8];
2710 xminTrack[1] =0; xmaxTrack[1]=3;
2711 axisName[1] ="pad type";
2714 xminTrack[2] =0; xmaxTrack[2]=1;
2715 axisName[2] ="drift length";
2718 xminTrack[3] =1; xmaxTrack[3]=400;
2719 axisName[3] ="Qmax";
2722 xminTrack[4] =0; xmaxTrack[4]=1;
2726 xminTrack[5] =0; xmaxTrack[5]=2;
2727 axisName[5] ="tan(phi)";
2730 xminTrack[6] =0; xmaxTrack[6]=2;
2731 axisName[6] ="tan(theta)";
2733 xminTrack[0] =-0.5; xmaxTrack[0]=0.5;
2734 fHisDeltaY=new THnSparseS("#Delta_{y} (cm)","#Delta_{y} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
2735 xminTrack[0] =-0.5; xmaxTrack[0]=0.5;
2736 fHisDeltaZ=new THnSparseS("#Delta_{z} (cm)","#Delta_{z} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
2737 xminTrack[0] =0.; xmaxTrack[0]=0.5;
2738 fHisRMSY=new THnSparseS("#RMS_{y} (cm)","#RMS_{y} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
2739 xminTrack[0] =0.; xmaxTrack[0]=0.5;
2740 fHisRMSZ=new THnSparseS("#RMS_{z} (cm)","#RMS_{z} (cm)", 7, binsTrack,xminTrack, xmaxTrack);
2741 xminTrack[0] =0.; xmaxTrack[0]=100;
2742 fHisQmax=new THnSparseS("Qmax (ADC)","Qmax (ADC)", 7, binsTrack,xminTrack, xmaxTrack);
2744 xminTrack[0] =0.; xmaxTrack[0]=250;
2745 fHisQtot=new THnSparseS("Qtot (ADC)","Qtot (ADC)", 7, binsTrack,xminTrack, xmaxTrack);
2746 BinLogX(fHisDeltaY,3);
2747 BinLogX(fHisDeltaZ,3);
2748 BinLogX(fHisRMSY,3);
2749 BinLogX(fHisRMSZ,3);
2750 BinLogX(fHisQmax,3);
2751 BinLogX(fHisQtot,3);
2755 void AliTPCcalibTracks::AddHistos(AliTPCcalibTracks* calib){
2759 if (calib->fHisDeltaY) fHisDeltaY->Add(calib->fHisDeltaY);
2760 if (calib->fHisDeltaZ) fHisDeltaZ->Add(calib->fHisDeltaZ);
2761 if (calib->fHisRMSY) fHisRMSY->Add(calib->fHisRMSY);
2762 if (calib->fHisRMSZ) fHisRMSZ->Add(calib->fHisRMSZ);