TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
AliTPCcalibTracks * calibTracks = ( AliTPCcalibTracks *)array->FindObject("calibTracks");
+
+//USAGE of debug stream example
+ gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
+ gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
+ AliXRDPROOFtoolkit tool;
+ TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
+ chainres->Lookup();
*/
#include <TCollection.h>
#include <iostream>
#include <TLinearFitter.h>
+#include <TString.h>
//
// AliROOT includes
#include "TText.h"
#include "TPaveText.h"
#include "TSystem.h"
+#include "TStatToolkit.h"
+#include "TCut.h"
// amplitude
sprintf(chname,"Amp_Sector%d",i);
- his1 = new TH1F(chname,chname,250,0,500); // valgrind
+ his1 = new TH1F(chname,chname,100,0,32); // valgrind
his1->SetXTitle("Max Amplitude (ADC)");
fArrayAmp->AddAt(his1,i);
sprintf(chname,"Amp_Sector%d",i+36);
- his1 = new TH1F(chname,chname,200,0,600); // valgrind 3 13,408,208 bytes in 229 blocks are still reachable
+ his1 = new TH1F(chname,chname,100,0,32); // valgrind 3 13,408,208 bytes in 229 blocks are still reachable
his1->SetXTitle("Max Amplitude (ADC)");
fArrayAmp->AddAt(his1,i+36);
// driftlength
sprintf(chname, "driftlengt vs. charge, ROC %i", i);
- prof1 = new TProfile(chname, chname, 500, 0, 250);
+ prof1 = new TProfile(chname, chname, 25, 0, 250);
prof1->SetYTitle("Charge");
prof1->SetXTitle("Driftlength");
fArrayChargeVsDriftlength->AddAt(prof1,i);
sprintf(chname, "driftlengt vs. charge, ROC %i", i+36);
- prof1 = new TProfile(chname, chname, 500, 0, 250);
+ prof1 = new TProfile(chname, chname, 25, 0, 250);
prof1->SetYTitle("Charge");
prof1->SetXTitle("Driftlength");
fArrayChargeVsDriftlength->AddAt(prof1,i+36);
for (Int_t ipad = 0; ipad < 3; ipad++){
Int_t bin = GetBin(iq, ipad);
Float_t qmean = GetQ(bin);
- char name[200];
- sprintf(name,"ResolY Pad%d Qmiddle%f",ipad, qmean);
- his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1);
+ char hname[200];
+ sprintf(hname,"ResolY Pad%d Qmiddle%f",ipad, qmean);
+ his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
fArrayQDY->AddAt(his3D, bin);
- sprintf(name,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
- his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, -1,1);
+ sprintf(hname,"ResolZ Pad%d Qmiddle%f",ipad, qmean);
+ his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, -1,1);
fArrayQDZ->AddAt(his3D, bin);
- sprintf(name,"RMSY Pad%d Qmiddle%f",ipad, qmean);
- his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1);
+ sprintf(hname,"RMSY Pad%d Qmiddle%f",ipad, qmean);
+ his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
fArrayQRMSY->AddAt(his3D, bin);
- sprintf(name,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
- his3D = new TH3F(name, name, 20,10,250, 20, 0,1.5, 50, 0,1);
+ sprintf(hname,"RMSZ Pad%d Qmiddle%f",ipad, qmean);
+ his3D = new TH3F(hname, hname, 20,10,250, 20, 0,1.5, 100, 0,0.6);
fArrayQRMSZ->AddAt(his3D, bin);
}
}
}
-void AliTPCcalibTracks::AddInfo(TChain * chain, char* fileName){
- //
- // Add the neccessary information for processing to the chain
- // (cluster parametrization)
- //
- TFile clusterParamFile(fileName);
- AliTPCClusterParam *clusterParam = (AliTPCClusterParam *) clusterParamFile.Get("Param");
- chain->GetUserInfo()->AddLast((TObject*)clusterParam);
- cout << "Clusterparametrization added to the chain." << endl;
-}
void AliTPCcalibTracks::Process(AliTPCseed *track){
//
//if (TMath::Abs(track->GetZ())<10.) return kFALSE;
//if (TMath::Abs(track->GetTgl())>0.03) return kFALSE;
- if (GetDebugLevel() > 5) Info("AcceptTrack","Track has been accepted.");
+ if (GetDebugLevel() > 20) Info("AcceptTrack","Track has been accepted.");
return 0;
}
if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal"," ***** Start of FillResolutionHistoLocal *****");
const Int_t kDelta = 10; // delta rows to fit
const Float_t kMinRatio = 0.75; // minimal ratio
- const Float_t kCutChi2 = 6.; // cut chi2 - left right - kink removal
+ // const Float_t kCutChi2 = 6.; // cut chi2 - left right - kink removal
const Float_t kErrorFraction = 0.5; // use only clusters with small interpolation error - for error param
const Int_t kFirstLargePad = 127; // medium pads -> long pads
const Float_t kLargePadSize = 1.5; // factor between medium and long pads' area
paramZ0 -= paramZ1;
matrixY0 += matrixY1;
matrixZ0 += matrixZ1;
- Double_t chi2 = 0;
+ Double_t cchi2 = 0;
TMatrixD difY(2, 1, paramY0.GetMatrixArray());
TMatrixD difYT(1, 2, paramY0.GetMatrixArray());
matrixY0.Invert();
TMatrixD mulY(matrixY0, TMatrixD::kMult, difY);
TMatrixD chi2Y(difYT, TMatrixD::kMult, mulY);
- chi2 += chi2Y(0, 0);
+ cchi2 += chi2Y(0, 0);
TMatrixD difZ(2, 1, paramZ0.GetMatrixArray());
TMatrixD difZT(1, 2, paramZ0.GetMatrixArray());
matrixZ0.Invert();
TMatrixD mulZ(matrixZ0, TMatrixD::kMult, difZ);
TMatrixD chi2Z(difZT, TMatrixD::kMult, mulZ);
- chi2 += chi2Z(0, 0);
+ cchi2 += chi2Z(0, 0);
// REMOVE KINK - TO be fixed - proper chi2 calculation for curved track to be implemented
//if (chi2 * 0.25 > kCutChi2) fRejectedTracksHisto->Fill(8);
// if this chi2 is bigger than a given threshold, assume that the current cluster is
// a kink an goto next padrow
- TTreeSRedirector *cstream = GetDebugStreamer();
if (cstream){
(*cstream)<<"Cut8"<<
- "chi2="<<chi2<<
+ "chi2="<<cchi2<<
"\n";
}
}
TProfile *profAmpRow = (TProfile*)fArrayAmpRow->At(sector);
if ( irow >= kFirstLargePad) max /= kLargePadSize;
- profAmpRow->Fill( (Double_t)cluster0->GetRow(), max );
+ Double_t smax = TMath::Sqrt(max);
+ profAmpRow->Fill( (Double_t)cluster0->GetRow(), smax );
TH1F *hisAmp = (TH1F*)fArrayAmp->At(sector);
- hisAmp->Fill(max);
+ hisAmp->Fill(smax);
// remove the following two lines one day:
TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector);
- profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max );
+ profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize));
- profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max );
+ profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
fHclusterPerPadrow->Fill(irow); // fill histogram showing clusters per padrow
Int_t ipad = TMath::Nint(cluster0->GetPad());
Float_t sy = cluster0->GetSigmaY2();
Float_t sz = cluster0->GetSigmaZ2();
(*cstream)<<"Resol0"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
"padSize="<<padSize<<
"angley="<<angley<<
"anglez="<<anglez<<
if (useForResol && nclFound > 2 * kMinRatio * kDelta
&& irow % kDeltaWriteDebugStream == 0 && GetDebugLevel() > 4){
- if (GetDebugLevel() > 5) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
+ if (GetDebugLevel() > 20) Info("FillResolutionHistoLocal","Filling 'TPCSelectorDebug.root', irow = %i", irow);
FillResolutionHistoLocalDebugPart(track, cluster0, irow, angley, anglez, nclFound, kDelta);
} // if (useForResol && nclFound > 2 * kMinRatio * kDelta)
TTreeSRedirector *cstream = GetDebugStreamer();
if (cstream){
(*cstream)<<"ResolCl"<< // valgrind 3 40,000 bytes in 1 blocks are possibly
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
"Sector="<<sector<<
"Cl.="<<cluster0<<
"CSigmaY0="<<csigmaY0<< // cluster errorY
// tracklet dubug
(*cstream)<<"ResolTr"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
"padSize="<<padSize<<
"IPad="<<padSize<<
"Sector="<<sector<<
}
-void AliTPCcalibTracks::MakeReport(Int_t stat, char* pathName){
+void AliTPCcalibTracks::MakeReport(Int_t stat, const char* pathName){
//
// all functions are called, that produce pictures
// the histograms are written to the directory 'pathName'
}
-void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, char* pathName){
+void AliTPCcalibTracks::MakeAmpPlots(Int_t stat, const char* pathName){
//
// creates several plots:
// fArrayAmp.ps, fArrayAmpRow.ps and DeltaYZ.ps
}
-void AliTPCcalibTracks::MakeDeltaPlots(char* pathName){
+void AliTPCcalibTracks::MakeDeltaPlots(const char* pathName){
//
// creates several plots:
// DeltaYZ.ps: DeltaY and DeltaZ histogram with gaus fit
}
-void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(char* pathName){
+void AliTPCcalibTracks::MakeChargeVsDriftLengthPlotsOld(const char* pathName){
//
// creates charge vs. driftlength plots, one TProfile for each ROC
// is not correct like this, should be one TProfile for each sector and padsize
}
-void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(char* pathName){
+void AliTPCcalibTracks::MakeChargeVsDriftLengthPlots(const char* pathName){
//
// creates charge vs. driftlength plots, one TProfile for each ROC
// under development....
-void AliTPCcalibTracks::FitResolutionNew(char* pathName){
+void AliTPCcalibTracks::FitResolutionNew(const char* pathName){
//
// calculates different resulution fits in Y and Z direction
// the histograms are written to 'ResolutionYZ.ps'
}
-void AliTPCcalibTracks::FitRMSNew(char* pathName){
+void AliTPCcalibTracks::FitRMSNew(const char* pathName){
//
// calculates different resulution-rms fits in Y and Z direction
// the histograms are written to 'RMS_YZ.ps'
}
-void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, char* pathName){
+void AliTPCcalibTracks::MakeResPlotsQTree(Int_t minEntries, const char* pathName){
//
// Make tree with resolution parameters
// the result is written to 'resol.root' in directory 'pathname'
AliTPCcalibTracks *calibTracks = 0;
if (GetDebugLevel() > 1) cout << "start to iterate, filling lists" << endl;
Int_t counter = 0;
- while ( (calibTracks = (AliTPCcalibTracks*)listIterator->Next()) ){
+ while ( (calibTracks = dynamic_cast<AliTPCcalibTracks*> (listIterator->Next())) ){
// loop over all entries in the collectionList and get dataMembers into lists
if (!calibTracks) continue;
clusterCutHistoList->Add(calibTracks->GetfClusterCutHisto());
hclusterPerPadrowList->Add(calibTracks->GetfHclusterPerPadrow());
hclusterPerPadrowRawList->Add(calibTracks->GetfHclusterPerPadrowRaw());
- fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
- fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
+ //
+ if (fCalPadClusterPerPad && calibTracks->GetfCalPadClusterPerPad())
+ fCalPadClusterPerPad->Add(calibTracks->GetfCalPadClusterPerPad());
+ // fCalPadClusterPerPadRaw->Add(calibTracks->GetfCalPadClusterPerPadRaw());
counter++;
if (GetDebugLevel() > 5) cout << "filling lists, object " << counter << " added." << endl;
}
if (GetDebugLevel() > 0) cout << "merging fArrayAmps..." << endl;
// merge fArrayAmps
for (Int_t i = 0; i < fArrayAmp->GetEntriesFast(); i++ ) { // loop over data member, i<72
- TIterator *objListIterator = arrayAmpList->MakeIterator();
+ TIterator *cobjListIterator = arrayAmpList->MakeIterator();
histList = new TList;
- while (( objarray = (TObjArray*)objListIterator->Next() )) {
+ while (( objarray = (TObjArray*)cobjListIterator->Next() )) {
// loop over arrayAmpList, get TObjArray, get object at position i, cast it into TH1F
if (!objarray) continue;
hist = (TH1F*)(objarray->At(i));
}
((TH1F*)(fArrayAmp->At(i)))->Merge(histList);
delete histList;
- delete objListIterator;
+ delete cobjListIterator;
}
if (GetDebugLevel() > 0) cout << "merging fArrayQDY..." << endl;
}
+void AliTPCcalibTracks::MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints){
+ //
+ // Make position corrections
+ // for the moment Only using debug streamer
+ // chainres - debug tree
+ // param - parameters to be updated
+ // maxPoints - maximal number of points using for fit
+ // verbose - print info flag
+ //
+ // Current implementation - only using debug streamers
+ //
+
+ /*
+ //Defaults
+ Int_t maxPoints=100000;
+ */
+ /*
+ Usage:
+ //0. Load libraries
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libSTAT");
+ gSystem->Load("libTPCcalib");
+
+
+ //1. Load Parameters to be modified:
+ //e.g:
+
+ AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ AliCDBManager::Instance()->SetRun(0)
+ AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
+
+ //2. Load chain from debug streamers
+ //
+ //e.g
+ gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
+ gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
+ AliXRDPROOFtoolkit tool;
+ TChain * chainres = tool.MakeChain("tracks.txt","ResolCl",0,10200);
+ chainres->Lookup();
+ //3. Do fits and store results
+ //
+ AliTPCcalibTracks::MakeQPosNormAll(chainres,param,200000,0) ;
+ TFile f("paramout.root","recreate");
+ param->Write("clusterParam");
+ f.Close();
+ //4. Verification
+ TFile f2("paramout.root");
+ AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
+ param2->SetInstance(param2);
+ 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
+
+ */
+
+
+ TStatToolkit toolkit;
+ Double_t chi2;
+ TVectorD fitParamY0;
+ TVectorD fitParamY1;
+ TVectorD fitParamZ0;
+ TVectorD fitParamZ1;
+ TMatrixD covMatrix;
+ Int_t npoints;
+
+ chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
+ chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
+ chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
+ chainres->SetAlias("st","(sin(dt)-dt)");
+ //
+ chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
+ //
+ //
+ //
+ TCut cutA("1");
+ TString fstringY="";
+ //
+ fstringY+="(dp)++"; //1
+ fstringY+="(dp)*di++"; //2
+ fstringY+="(sp)++"; //3
+ fstringY+="(sp)*di++"; //4
+ TString fstringZ="";
+ fstringZ+="(dt)++"; //1
+ fstringZ+="(dt)*di++"; //2
+ fstringZ+="(st)++"; //3
+ fstringZ+="(st)*di++"; //4
+ //
+ // Z corrections
+ //
+ TString *strZ0 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
+ printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
+ param->fPosZcor[0] = (TVectorD*) fitParamZ0.Clone();
+ //
+ TString *strZ1 = toolkit.FitPlane(chainres,"(Cl.fZ-PZ0.fElements[0]):CSigmaZ0",fstringZ.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
+ //
+ printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
+ param->fPosZcor[1] = (TVectorD*) fitParamZ1.Clone();
+ param->fPosZcor[2] = (TVectorD*) fitParamZ1.Clone();
+ //
+ // Y corrections
+ //
+ TString *strY0 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector<36"+cutA, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
+ printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
+ param->fPosYcor[0] = (TVectorD*) fitParamY0.Clone();
+
+
+ TString *strY1 = toolkit.FitPlane(chainres,"(Cl.fY-PY0.fElements[0]):CSigmaY0",fstringY.Data(), "Cl.fDetector>36"+cutA, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
+ //
+ printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
+ param->fPosYcor[1] = (TVectorD*) fitParamY1.Clone();
+ param->fPosYcor[2] = (TVectorD*) fitParamY1.Clone();
+ //
+ //
+ //
+ chainres->SetAlias("fitZ0",strZ0->Data());
+ chainres->SetAlias("fitZ1",strZ1->Data());
+ chainres->SetAlias("fitY0",strY0->Data());
+ chainres->SetAlias("fitY1",strY1->Data());
+ // chainres->Draw("Cl.fZ-PZ0.fElements[0]","CSigmaY0<0.7&&CSigmaZ0<0.7"+cutA,"",10000);
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff