// the fitters after loading this object from file.
// (This will be gone for a new ROOT version > v5-17-05)
//
+//
+// In order to debug some numerical algorithm all data data which are used for
+// fitters can be stored in the debug streamers. In case of fitting roblems the
+// errors and tendencies can be checked.
+//
+// Debug Streamers:
+//
+//
+//
+//
+//
////////////////////////////////////////////////////////////////////////////
+/*
+ .x ~/UliStyle.C
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libSTAT");
+ gSystem->Load("libTPCcalib");
+
+ TFile fcalib("CalibObjects.root");
+ TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
+ AliTPCcalibTracksGain * gain = ( AliTPCcalibTracksGain *)array->FindObject("calibTracksGain");
+
+ //
+ // Angular and drift correction
+ //
+ AliTPCClusterParam *param = new AliTPCClusterParam;param->SetInstance(param);
+ gain->UpdateClusterParam(param);
+ TF2 fdrty("fdrty","AliTPCClusterParam::SQnorm(0,0,x,y,0)",0,1,0,1)
+
+ //
+ // Make visual Tree - compare with Kr calibration
+ //
+ AliTPCCalPad * gain010 = gain->CreateFitCalPad(0,kTRUE,0); gain010->SetName("CGain010");
+ AliTPCCalPad * gain110 = gain->CreateFitCalPad(1,kTRUE,0); gain110->SetName("CGain110");
+ AliTPCCalPad * gain210 = gain->CreateFitCalPad(2,kTRUE,0); gain210->SetName("CGain210");
+ TFile fkr("/u/miranov/GainMap.root");
+ AliTPCCalPad *gainKr = fkr.Get("GainMap"); fkr->SetName("KrGain");
+ //
+ AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
+ preprocesor->AddComponent(gain010);
+ preprocesor->AddComponent(gain110);
+ preprocesor->AddComponent(gain210);
+ preprocesor->AddComponent(gainKr);
+ preprocesor->DumpToFile("cosmicGain.root");
+ //
+ //
+ //
+ // Simple session using the debug streamers
+ //
+
+ gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
+ gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
+ AliXRDPROOFtoolkit tool;
+
+ TChain * chain0 = tool.MakeChain("gain.txt","dEdx",0,1000000);
+ TChain * chain1 = tool.MakeChain("gain.txt","Track",0,1000000);
+ TChain * chain2 = tool.MakeChain("gain.txt","TrackG",0,1000000);
+ chain0->Lookup();
+ chain1->Lookup();
+ chain2->Lookup();
+
+ chain2->SetAlias("k1","1/0.855");
+ chain2->SetAlias("k0","1/0.9928");
+ chain2->SetAlias("k2","1/1.152");
+
+*/
+
+
+
+
+
#include <TPDGCode.h>
#include <TStyle.h>
#include "TSystem.h"
#include <TFile.h>
#include <TCollection.h>
#include <TIterator.h>
+#include <TProfile.h>
+#include <TProfile2D.h>
+#include <TProof.h>
+#include <TStatToolkit.h>
//
// AliRoot includes
#include "AliTPCclusterMI.h"
#include "AliTPCcalibTracksCuts.h"
#include "AliTPCFitPad.h"
+#include "TStatToolkit.h"
+#include "TString.h"
+#include "TCut.h"
-// REMOVE ALL OF THIS
+//
#include <TTree.h>
#include "AliESDEvent.h"
AliTPCParamSR* AliTPCcalibTracksGain::fgTPCparam = new AliTPCParamSR();
AliTPCcalibTracksGain::AliTPCcalibTracksGain() :
- TNamed(),
- fDebugStream(0), //! debug stream for debugging
+ AliTPCcalibBase(),
fCuts(0), // cuts that are used for sieving the tracks used for calibration
+ fGainMap(0), // gain map to be applied
//
// Simple Array of histograms
//
fSimpleFitter(0), // simple fitter for short pads
fSqrtFitter(0), // sqrt fitter for medium pads
fLogFitter(0), // log fitter for long pads
+
fFitter0M(0), // fitting of the atenuation, angular correction, and mean chamber gain
fFitter1M(0), // fitting of the atenuation, angular correction, and mean chamber gain
fFitter2M(0), // fitting of the atenuation, angular correction, and mean chamber gain
fFitter0T(0), // fitting of the atenuation, angular correction, and mean chamber gain
fFitter1T(0), // fitting of the atenuation, angular correction, and mean chamber gain
fFitter2T(0), // fitting of the atenuation, angular correction, and mean chamber gain
+ //
+ fDFitter0M(0), // fitting of the atenuation, angular correction
+ fDFitter1M(0), // fitting of the atenuation, angular correction
+ fDFitter2M(0), // fitting of the atenuation, angular correction
+ fDFitter0T(0), // fitting of the atenuation, angular correction
+ fDFitter1T(0), // fitting of the atenuation, angular correction
+ fDFitter2T(0), // fitting of the atenuation, angular correction
+ //
fSingleSectorFitter(0), // just for debugging
//
// Counters
fTotalTracks(0), // just for debugging
fAcceptedTracks(0), // just for debugging
fDebugCalPadRaw(0), // just for debugging
- fDebugCalPadCorr(0), // just for debugging
- fPrevIter(0) // the calibration object in its previous iteration (will not be owned by the new object, don't forget to delete it!)
+ fDebugCalPadCorr(0) // just for debugging
{
//
}
AliTPCcalibTracksGain::AliTPCcalibTracksGain(const AliTPCcalibTracksGain& obj) :
- TNamed(obj),
- fDebugStream(0), //! debug stream for debugging
+ AliTPCcalibBase(obj),
fCuts(obj.fCuts), // cuts that are used for sieving the tracks used for calibration
+ fGainMap(new AliTPCCalPad(*(obj.fGainMap))), // gain map to be applied
fArrayQM(0), // Qmax normalized
fArrayQT(0), // Qtot normalized
//
fFitter0T(obj.fFitter0T),
fFitter1T(obj.fFitter1T),
fFitter2T(obj.fFitter2T),
+ //
+ fDFitter0M(obj.fDFitter0M),
+ fDFitter1M(obj.fDFitter1M),
+ fDFitter2M(obj.fDFitter2M),
+ fDFitter0T(obj.fDFitter0T),
+ fDFitter1T(obj.fDFitter1T),
+ fDFitter2T(obj.fDFitter2T),
fSingleSectorFitter(obj.fSingleSectorFitter), // just for debugging
//
// Counters
fTotalTracks(obj.fTotalTracks), // just for debugging
fAcceptedTracks(obj.fAcceptedTracks), // just for debugging
fDebugCalPadRaw(obj.fDebugCalPadRaw), // just for debugging
- fDebugCalPadCorr(obj.fDebugCalPadCorr), // just for debugging
- fPrevIter(obj.fPrevIter) // the calibration object in its previous iteration (will not be owned by the new object, don't forget to delete it!)
+ fDebugCalPadCorr(obj.fDebugCalPadCorr) // just for debugging
{
//
fSqrtFitter = new AliTPCFitPad(*(rhs.fSqrtFitter));
fLogFitter = new AliTPCFitPad(*(rhs.fLogFitter));
fSingleSectorFitter = new AliTPCFitPad(*(rhs.fSingleSectorFitter));
- fPrevIter = new AliTPCcalibTracksGain(*(rhs.fPrevIter));
fCuts = new AliTPCcalibTracksCuts(*(rhs.fCuts));
+ fGainMap = new AliTPCCalPad(*(rhs.fGainMap));
}
return *this;
}
-AliTPCcalibTracksGain::AliTPCcalibTracksGain(const char* name, const char* title, AliTPCcalibTracksCuts* cuts, TNamed* /*debugStreamPrefix*/, AliTPCcalibTracksGain* prevIter) :
- TNamed(name, title),
- fDebugStream(0), //! debug stream for debugging
+AliTPCcalibTracksGain::AliTPCcalibTracksGain(const char* name, const char* title, AliTPCcalibTracksCuts* cuts) :
+ AliTPCcalibBase(),
fCuts(0), // cuts that are used for sieving the tracks used for calibration
+ fGainMap(0), // gain map to be applied
fArrayQM(0), // Qmax normalized
fArrayQT(0), // Qtot normalized
//
fFitter0T(0), // fitting of the atenuation, angular correction, and mean chamber gain
fFitter1T(0), // fitting of the atenuation, angular correction, and mean chamber gain
fFitter2T(0), // fitting of the atenuation, angular correction, and mean chamber gain
+ //
+ fDFitter0M(0), // fitting of the atenuation, angular correction
+ fDFitter1M(0), // fitting of the atenuation, angular correction
+ fDFitter2M(0), // fitting of the atenuation, angular correction
+ fDFitter0T(0), // fitting of the atenuation, angular correction
+ fDFitter1T(0), // fitting of the atenuation, angular correction
+ fDFitter2T(0), // fitting of the atenuation, angular correction
fSingleSectorFitter(0), // just for debugging
//
// Counters
fTotalTracks(0), // just for debugging
fAcceptedTracks(0), // just for debugging
fDebugCalPadRaw(0), // just for debugging
- fDebugCalPadCorr(0), // just for debugging
- fPrevIter(0) // the calibration object in its previous iteration (will not be owned by the new object, don't forget to delete it!)
+ fDebugCalPadCorr(0) // just for debugging
{
//
// Constructor.
//
- G__SetCatchException(0);
+ this->SetNameTitle(name, title);
fCuts = cuts;
- fPrevIter = prevIter;
//
// Fitter initialization
//
fFitter1T = new TLinearFitter(45,"hyp44");
fFitter2T = new TLinearFitter(45,"hyp44");
//
+ fDFitter0M = new TLinearFitter(10,"hyp9");
+ fDFitter1M = new TLinearFitter(10,"hyp9");
+ fDFitter2M = new TLinearFitter(10,"hyp9");
+ fDFitter0T = new TLinearFitter(10,"hyp9");
+ fDFitter1T = new TLinearFitter(10,"hyp9");
+ fDFitter2T = new TLinearFitter(10,"hyp9");
+ //
+ //
+ fFitter0M->StoreData(kFALSE);
+ fFitter1M->StoreData(kFALSE);
+ fFitter2M->StoreData(kFALSE);
+ fFitter0T->StoreData(kFALSE);
+ fFitter1T->StoreData(kFALSE);
+ fFitter2T->StoreData(kFALSE);
+ //
+ fDFitter0M->StoreData(kFALSE);
+ fDFitter1M->StoreData(kFALSE);
+ fDFitter2M->StoreData(kFALSE);
+ fDFitter0T->StoreData(kFALSE);
+ fDFitter1T->StoreData(kFALSE);
+ fDFitter2T->StoreData(kFALSE);
+ //
//
// Add profile histograms -JUST for visualization - Not used for real calibration
//
if (fLogFitter) delete fLogFitter;
if (fSingleSectorFitter) delete fSingleSectorFitter;
-
- if (fDebugStream) {
- //fDebugStream->GetFile()->Close();
- printf("Deleting debug stream object\n");
- delete fDebugStream;
- }
-
-
if (fDebugCalPadRaw) delete fDebugCalPadRaw;
if (fDebugCalPadCorr) delete fDebugCalPadCorr;
}
//
Evaluate();
- if (fDebugStream) {
- delete fDebugStream;
- fDebugStream = 0;
- }
+ AliTPCcalibBase::Terminate();
}
-void AliTPCcalibTracksGain::AddInfo(TChain* chain, char* debugStreamPrefix, char* prevIterFileName) {
- //
- // Add some parameters to the chain.
- // debugStreamPrefix: If specified, contains the location (either normal or xrootd directory)
- // where the debug stream is moved (normal directory) or copied to (xrootd).
- // prevIterFileName: If specified, contains an AliTPCcalibTracksGain object from a previous run
- // for doing an iterative calibration procedure (right now unused).
- // Note: The parameters are *not* added to this class, you need to do it later by retrieving
- // the parameters from the chain and passing them to the constructor!
- //
-
- if (debugStreamPrefix) {
- TNamed* objDebugStreamPrefix = new TNamed("debugStreamPrefix", debugStreamPrefix);
- chain->GetUserInfo()->AddLast((TObject*)objDebugStreamPrefix);
- }
-
- if (prevIterFileName) {
- TFile paramFile(prevIterFileName);
- if (paramFile.IsZombie()) {
- printf("File %s not found. Continuing without previous iteration.\n", prevIterFileName);
- return;
- }
-
- AliTPCcalibTracksGain *prevIter = (AliTPCcalibTracksGain*)paramFile.Get("calibTracksGain");
- if (prevIter) {
- chain->GetUserInfo()->AddLast((TObject*)prevIter);
- } else
- printf("No calibTracksGain object found. Continuing without previous iteration.\n");
- }
-}
-Bool_t AliTPCcalibTracksGain::AcceptTrack(AliTPCseed* track) {
- //
- // Decides whether to accept a track or not depending on track parameters and cuts
- // contained as AliTPCcalibTracksCuts object fCuts.
- // Tracks are discarded if the number of clusters is too low or the transverse
- // momentum is too low.
- // The corresponding cut values are specified in the fCuts member.
- //
-
- if (track->GetNumberOfClusters() < fCuts->GetMinClusters()) return kFALSE;
- //if ((TMath::Abs(track->GetY() / track->GetX()) > fCuts->GetEdgeYXCutNoise())
- // && (TMath::Abs(track->GetTgl()) < fCuts->GetEdgeThetaCutNoise())) return kFALSE;
- //if (track->GetNumberOfClusters() / (track->GetNFoundable()+1.) < fCuts->GetMinRatio()) return kFALSE;
- if (TMath::Abs(track->GetSigned1Pt()) > fCuts->GetMax1pt()) return kFALSE;
-
- //if (track->GetPt() < 50.) return kFALSE;
- return kTRUE;
-}
void AliTPCcalibTracksGain::Process(AliTPCseed* seed) {
//
// is added.
//
+
fTotalTracks++;
- if (!AcceptTrack(seed)) return;
+ if (!fCuts->AcceptTrack(seed)) return;
+ //
+ // reinint on proof
+ // if (gProof){
+ static Bool_t doinit= kTRUE;
+ if (doinit){
+ fSimpleFitter = new AliTPCFitPad(8, "hyp7", "");
+ fSqrtFitter = new AliTPCFitPad(8, "hyp7", "");
+ fLogFitter = new AliTPCFitPad(8, "hyp7", "");
+ fSingleSectorFitter = new AliTPCFitPad(8, "hyp7", "");
+ //
+ fFitter0M = new TLinearFitter(45,"hyp44");
+ fFitter1M = new TLinearFitter(45,"hyp44");
+ fFitter2M = new TLinearFitter(45,"hyp44");
+ fFitter0T = new TLinearFitter(45,"hyp44");
+ fFitter1T = new TLinearFitter(45,"hyp44");
+ fFitter2T = new TLinearFitter(45,"hyp44");
+ //
+ fDFitter0M = new TLinearFitter(10,"hyp9");
+ fDFitter1M = new TLinearFitter(10,"hyp9");
+ fDFitter2M = new TLinearFitter(10,"hyp9");
+ fDFitter0T = new TLinearFitter(10,"hyp9");
+ fDFitter1T = new TLinearFitter(10,"hyp9");
+ fDFitter2T = new TLinearFitter(10,"hyp9");
+ doinit=kFALSE;
+ }
+ //}
+
fAcceptedTracks++;
AddTrack(seed);
}
while ((cal = (AliTPCcalibTracksGain*)iter->Next())) {
if (!cal->InheritsFrom(AliTPCcalibTracksGain::Class())) {
- Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
+ //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
return -1;
}
return 0;
}
+Float_t AliTPCcalibTracksGain::GetGain(AliTPCclusterMI* cl){
+ //
+ // Return local gain at cluster position
+ //
+ Float_t factor = 1;
+ if(!fGainMap) return factor;
+
+ AliTPCCalROC * roc = fGainMap->GetCalROC(cl->GetDetector());
+ Int_t irow = cl->GetRow();
+ if (roc){
+ if (irow < 63) { // IROC
+ factor = roc->GetValue(irow, TMath::Nint(cl->GetPad()));
+ } else { // OROC
+ factor = roc->GetValue(irow - 63, TMath::Nint(cl->GetPad()));
+ }
+ }
+ if (factor<0.1) factor=0.1;
+ return factor;
+}
+
+
+Float_t AliTPCcalibTracksGain::GetMaxNorm(AliTPCclusterMI * cl){
+ //
+ // Get normalized amplituded if the gain map provided
+ //
+ return cl->GetMax()/GetGain(cl);
+}
+
+
+Float_t AliTPCcalibTracksGain::GetQNorm(AliTPCclusterMI * cl){
+ //
+ // Get normalized amplituded if the gain map provided
+ //
+ return cl->GetQ()/GetGain(cl);
+}
+
+
+
void AliTPCcalibTracksGain::Add(AliTPCcalibTracksGain* cal) {
//
// Adds another AliTPCcalibTracksGain object to this object.
//
- fSimpleFitter->Add(cal->fSimpleFitter);
- fSqrtFitter->Add(cal->fSqrtFitter);
- fLogFitter->Add(cal->fLogFitter);
- fSingleSectorFitter->Add(cal->fSingleSectorFitter);
- //
- //
- //
- fFitter0M->Add(cal->fFitter0M);
- fFitter1M->Add(cal->fFitter1M);
- fFitter2M->Add(cal->fFitter2M);
- fFitter0T->Add(cal->fFitter0T);
- fFitter1T->Add(cal->fFitter1T);
- fFitter2T->Add(cal->fFitter2T);
+ fSimpleFitter->Add(cal->fSimpleFitter);
+ fSqrtFitter->Add(cal->fSqrtFitter);
+ fLogFitter->Add(cal->fLogFitter);
+ fSingleSectorFitter->Add(cal->fSingleSectorFitter);
+ //
+ //
+ //
+ if (cal->fFitter0M->GetNpoints()>0) fFitter0M->Add(cal->fFitter0M);
+ if (cal->fFitter1M->GetNpoints()>0) fFitter1M->Add(cal->fFitter1M);
+ if (cal->fFitter2M->GetNpoints()>0) fFitter2M->Add(cal->fFitter2M);
+ if (cal->fFitter0T->GetNpoints()>0) fFitter0T->Add(cal->fFitter0T);
+ if (cal->fFitter1T->GetNpoints()>0) fFitter1T->Add(cal->fFitter1T);
+ if (cal->fFitter2T->GetNpoints()>0) fFitter2T->Add(cal->fFitter2T);
+ //
+ if (cal->fDFitter0M->GetNpoints()>0) fDFitter0M->Add(cal->fDFitter0M);
+ if (cal->fDFitter1M->GetNpoints()>0) fDFitter1M->Add(cal->fDFitter1M);
+ if (cal->fDFitter2M->GetNpoints()>0) fDFitter2M->Add(cal->fDFitter2M);
+ if (cal->fDFitter0T->GetNpoints()>0) fDFitter0T->Add(cal->fDFitter0T);
+ if (cal->fDFitter1T->GetNpoints()>0) fDFitter1T->Add(cal->fDFitter1T);
+ if (cal->fDFitter2T->GetNpoints()>0) fDFitter2T->Add(cal->fDFitter2T);
//
//
// histograms
// See AddCluster(...) for more detail.
//
- if (!fDebugStream) fDebugStream = new TTreeSRedirector(fgkDebugStreamFileName);
DumpTrack(seed);
//
// simple histograming part
void AliTPCcalibTracksGain::AddCluster(AliTPCclusterMI* cluster){
//
// Adding cluster information to the simple histograms
- // No correction, fittings are applied
- //
+ // No correction, fittings are applied
+
+
Float_t kThreshold=5;
if (cluster->GetX()<=0) return;
if (cluster->GetQ()<=kThreshold) return;
//
-
Int_t sector = cluster->GetDetector();
TH1F * his;
his = GetQT(sector);
- if (his) his->Fill(cluster->GetQ());
+ if (his) his->Fill(GetQNorm(cluster));
his = GetQT(-1);
- if (his) his->Fill(cluster->GetQ());
+ if (his) his->Fill(GetQNorm(cluster));
his = GetQM(sector);
- if (his) his->Fill(cluster->GetMax());
+ if (his) his->Fill(GetMaxNorm(cluster));
his = GetQM(-1);
- if (his) his->Fill(cluster->GetMax());
+ if (his) his->Fill(GetMaxNorm(cluster));
//
sector = sector%36;
TProfile *prof;
prof = GetProfileQT(sector);
- if (prof) prof->Fill(cluster->GetX(),cluster->GetQ());
+ if (prof) prof->Fill(cluster->GetX(),GetQNorm(cluster));
prof = GetProfileQT(-1);
- if (prof) prof->Fill(cluster->GetX(),cluster->GetQ());
+ if (prof) prof->Fill(cluster->GetX(),GetQNorm(cluster));
prof = GetProfileQM(sector);
- if (prof) prof->Fill(cluster->GetX(),cluster->GetMax());
+ if (prof) prof->Fill(cluster->GetX(),GetMaxNorm(cluster));
prof = GetProfileQM(-1);
- if (prof) prof->Fill(cluster->GetX(),cluster->GetMax());
+ if (prof) prof->Fill(cluster->GetX(),GetMaxNorm(cluster));
//
Float_t phi = cluster->GetY()/cluster->GetX();
TProfile2D *prof2;
prof2 = GetProfileQT2D(sector);
- if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetQ());
+ if (prof2) prof2->Fill(cluster->GetX(),phi,GetQNorm(cluster));
prof2 = GetProfileQT2D(-1);
- if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetQ());
+ if (prof2) prof2->Fill(cluster->GetX(),phi,GetQNorm(cluster));
prof2 = GetProfileQM2D(sector);
- if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetMax());
+ if (prof2) prof2->Fill(cluster->GetX(),phi,GetMaxNorm(cluster));
prof2 = GetProfileQM2D(-1);
- if (prof2) prof2->Fill(cluster->GetX(),phi,cluster->GetMax());
+ if (prof2) prof2->Fill(cluster->GetX(),phi,GetMaxNorm(cluster));
//
}
// gets the square root of the charge, and the log fitter gets fgkM*(1+q/fgkM), where q is the original charge
// and fgkM==25.
//
+ Float_t kedge = 3;
+ Float_t kfraction = 0.7;
+ Int_t kinorm = 2;
+
+
+ // Where to put selection on threshold?
+ // Defined by the Q/dEdxT variable - see debug streamer:
+ //
+ // Debug stream variables: (Where tu cut ?)
+ // chain0->Draw("Cl.fQ/dedxQ.fElements[1]>>his(100,0,3)","fraction2<0.6&&dedge>3","",1000000);
+ // mean 1 sigma 0.25
+ // chain0->Draw("Cl.fMax/dedxM.fElements[1]>>his(100,0,3)","fraction2<0.6&&dedge>3","",1000000)
+ // mean 1 sigma 0.25
+ // chain0->Draw("Cl.fQ/dedxQ.fElements[2]>>his(100,0,3)","fraction2<0.7&&dedge>3","",1000000)
+ // mean 1 sigma 0.29
+ // chain0->Draw("Cl.fMax/dedxM.fElements[2]>>his(100,0,3)","fraction2<0.7&&dedge>3","",1000000)
+ // mean 1 sigma 0.27
+ // chain0->Draw("Cl.fQ/dedxQ.fElements[3]>>his(100,0,3)","fraction2<0.8&&dedge>3","",1000000)
+ // mean 1 sigma 0.32
+ //
+ // chain0->Draw("Cl.fQ/dedxQ.fElements[4]>>his(100,0,3)","fraction2<0.9&&dedge>3","",1000000)
+ // mean 1 sigma 0.4
+
+ // Fraction choosen 0.7
if (!cluster) {
Error("AddCluster", "Cluster not valid.");
return;
}
- if (dedge < 3.) return;
- if (fraction2 > 0.7) return;
+ if (dedge < kedge) return;
+ if (fraction2 > kfraction) return;
//Int_t padType = GetPadType(cluster->GetX());
Double_t xx[7];
// Update fitters
//
Int_t segment = cluster->GetDetector() % 36;
- Double_t q = fgkUseTotalCharge ? ((Double_t)(cluster->GetQ())) : ((Double_t)(cluster->GetMax())); // note: no normalization to pad size!
+ Double_t q = fgkUseTotalCharge ? ((Double_t)(GetQNorm(cluster))) : ((Double_t)(GetMaxNorm(cluster))); // note: no normalization to pad size!
// just for debugging
Int_t row = 0;
fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadRaw->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
// correct charge by normalising to mean charge per track
- q /= dedxQ[2];
+ q /= dedxQ[kinorm];
// just for debugging
fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->SetValue(row, pad, q + fDebugCalPadCorr->GetCalROC(cluster->GetDetector())->GetValue(row, pad));
Double_t sqrtQ = TMath::Sqrt(q);
Double_t logQ = fgkM * TMath::Log(1 + q / fgkM);
- fSimpleFitter->GetFitter(segment, padType)->AddPoint(xx, q);
- fSqrtFitter->GetFitter(segment, padType)->AddPoint(xx, sqrtQ);
- fLogFitter->GetFitter(segment, padType)->AddPoint(xx, logQ);
- fSingleSectorFitter->GetFitter(0, padType)->AddPoint(xx, q);
+ TLinearFitter * fitter =0;
+ //
+ fitter = fSimpleFitter->GetFitter(segment, padType);
+ fitter->AddPoint(xx, q);
+ //
+ fitter = fSqrtFitter->GetFitter(segment, padType);
+ fitter->AddPoint(xx, sqrtQ);
+ //
+ fitter = fLogFitter->GetFitter(segment, padType);
+ fitter->AddPoint(xx, logQ);
+ //
+ fitter=fSingleSectorFitter->GetFitter(0, padType);
+ fitter->AddPoint(xx, q);
// this will be gone for the a new ROOT version > v5-17-05
if (padType == kShortPads)
fFitter0T->Eval();
fFitter1T->Eval();
fFitter2T->Eval();
+ //
+ fDFitter0M->Eval();
+ fDFitter1M->Eval();
+ fDFitter2M->Eval();
+ fDFitter0T->Eval();
+ fDFitter1T->Eval();
+ fDFitter2T->Eval();
}
AliTPCCalPad* AliTPCcalibTracksGain::CreateFitCalPad(UInt_t fitType, Bool_t undoTransformation, Bool_t normalizeToPadSize) {
return roc;
}
-void AliTPCcalibTracksGain::GetParameters(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitParam) {
+Bool_t AliTPCcalibTracksGain::GetParameters(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitParam) {
//
// Puts the fit parameters for the specified segment (IROC & OROC), padType and fitType
// into the fitParam TVectorD (which should contain 8 elements).
// padType is one of kShortPads, kMediumPads, kLongPads. fitType is one of kSimpleFitter, kSqrtFitter, kLogFitter.
// Note: The fitter has to be evaluated first!
//
-
- GetFitter(segment, padType, fitType)->GetParameters(fitParam);
+ TLinearFitter * fitter = GetFitter(segment, padType, fitType);
+ if (fitter){
+ fitter->Eval();
+ fitter->GetParameters(fitParam);
+ return kTRUE;
+ }else{
+ Error("AliTPCcalibTracksGain::GetParameters",
+ Form("Fitter%d_%d_%d not availble", segment, padType, fitType));
+ return kFALSE;
+ }
+ return kFALSE;
}
void AliTPCcalibTracksGain::GetErrors(UInt_t segment, UInt_t padType, UInt_t fitType, TVectorD &fitError) {
AddTracklet(sector[ipad],ipad, dedxQ[ipad], dedxM[ipad], parY[ipad], parZ[ipad], meanPos[ipad] );
if (isOK) count++;
}
-
- (*fDebugStream) << "Track" <<
- "Track.=" << track << // track information
- "\n";
- //
- //
- //
- if (count>1){
- (*fDebugStream) << "TrackG" <<
+
+ TTreeSRedirector * cstream = GetDebugStreamer();
+ if (cstream){
+ (*cstream) << "Track" <<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
"Track.=" << track << // track information
- "ncount="<<count<<
- // info for pad type 0
- "sector0="<<sector[0]<<
- "npoints0="<<npoints[0]<<
- "dedxM0.="<<&dedxM[0]<<
- "dedxQ0.="<<&dedxQ[0]<<
- "parY0.="<<&parY[0]<<
- "parZ0.="<<&parZ[0]<<
- "meanPos0.="<<&meanPos[0]<<
- //
- // info for pad type 1
- "sector1="<<sector[1]<<
- "npoints1="<<npoints[1]<<
- "dedxM1.="<<&dedxM[1]<<
- "dedxQ1.="<<&dedxQ[1]<<
- "parY1.="<<&parY[1]<<
- "parZ1.="<<&parZ[1]<<
- "meanPos1.="<<&meanPos[1]<<
- //
- // info for pad type 2
- "sector2="<<sector[2]<<
- "npoints2="<<npoints[2]<<
- "dedxM2.="<<&dedxM[2]<<
- "dedxQ2.="<<&dedxQ[2]<<
- "parY2.="<<&parY[2]<<
- "parZ2.="<<&parZ[2]<<
- "meanPos2.="<<&meanPos[2]<<
- //
"\n";
+ //
+ //
+ //
+ if ( GetStreamLevel()>1 && count>1){
+ (*cstream) << "TrackG" <<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ "Track.=" << track << // track information
+ "ncount="<<count<<
+ // info for pad type 0
+ "sector0="<<sector[0]<<
+ "npoints0="<<npoints[0]<<
+ "dedxM0.="<<&dedxM[0]<<
+ "dedxQ0.="<<&dedxQ[0]<<
+ "parY0.="<<&parY[0]<<
+ "parZ0.="<<&parZ[0]<<
+ "meanPos0.="<<&meanPos[0]<<
+ //
+ // info for pad type 1
+ "sector1="<<sector[1]<<
+ "npoints1="<<npoints[1]<<
+ "dedxM1.="<<&dedxM[1]<<
+ "dedxQ1.="<<&dedxQ[1]<<
+ "parY1.="<<&parY[1]<<
+ "parZ1.="<<&parZ[1]<<
+ "meanPos1.="<<&meanPos[1]<<
+ //
+ // info for pad type 2
+ "sector2="<<sector[2]<<
+ "npoints2="<<npoints[2]<<
+ "dedxM2.="<<&dedxM[2]<<
+ "dedxQ2.="<<&dedxQ[2]<<
+ "parY2.="<<&parY[2]<<
+ "parZ2.="<<&parZ[2]<<
+ "meanPos2.="<<&meanPos[2]<<
+ //
+ "\n";
+ }
}
-
}
Bool_t AliTPCcalibTracksGain::GetDedx(AliTPCseed* track, Int_t padType, Int_t* /*rows*/,
static TLinearFitter fitY(2, "pol1");
static TLinearFitter fitZ(2, "pol1");
+ fitY.StoreData(kFALSE);
+ fitZ.StoreData(kFALSE);
//
//
Int_t firstRow = 0, lastRow = 0;
Int_t detector = cluster->GetDetector() ;
if (lastSector == -1) lastSector = detector;
if (lastSector != detector) continue;
- amplitudeQ[nclusters] = cluster->GetQ();
- amplitudeM[nclusters] = cluster->GetMax();
+ amplitudeQ[nclusters] = GetQNorm(cluster);
+ amplitudeM[nclusters] = GetMaxNorm(cluster);
rowIn[nclusters] = iCluster;
nclusters++;
Double_t dx = cluster->GetX() - xcenter;
dedxQ[i] /= ndedx[i];
dedxM[i] /= ndedx[i];
}
-
+ TTreeSRedirector * cstream = GetDebugStreamer();
inonEdge = 0;
Float_t momenta = track->GetP();
Float_t mdedx = track->GetdEdx();
Float_t fraction2 = Float_t(inonEdge) / Float_t(nclustersNE);
AddCluster(cluster, momenta, mdedx, padType, xcenter, dedxQ, dedxM, fraction, fraction2, dedge, parY, parZ, meanPos);
-
- (*fDebugStream) << "dEdx" <<
- "Cl.=" << cluster << // cluster of interest
- "P=" << momenta << // track momenta
- "dedx=" << mdedx << // mean dedx - corrected for angle
- "IPad=" << padType << // pad type 0..2
- "xc=" << xcenter << // x center of chamber
- "dedxQ.=" << &dedxQ << // dedxQ - total charge
- "dedxM.=" << &dedxM << // dedxM - maximal charge
- "fraction=" << fraction << // fraction - order in statistic (0,1)
- "fraction2=" << fraction2 << // fraction - order in statistic (0,1)
- "dedge=" << dedge << // distance to the edge
- "parY.=" << &parY << // line fit
- "parZ.=" << &parZ << // line fit
- "meanPos.=" << &meanPos << // mean position (dx, dx^2, y,y^2, z, z^2)
- "\n";
+ Float_t gain = GetGain(cluster);
+ if (cstream) (*cstream) << "dEdx" <<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+
+ "Cl.=" << cluster << // cluster of interest
+ "gain="<<gain<< // gain at cluster position
+ "P=" << momenta << // track momenta
+ "dedx=" << mdedx << // mean dedx - corrected for angle
+ "IPad=" << padType << // pad type 0..2
+ "xc=" << xcenter << // x center of chamber
+ "dedxQ.=" << &dedxQ << // dedxQ - total charge
+ "dedxM.=" << &dedxM << // dedxM - maximal charge
+ "fraction=" << fraction << // fraction - order in statistic (0,1)
+ "fraction2=" << fraction2 << // fraction - order in statistic (0,1)
+ "dedge=" << dedge << // distance to the edge
+ "parY.=" << &parY << // line fit
+ "parZ.=" << &parZ << // line fit
+ "meanPos.=" << &meanPos << // mean position (dx, dx^2, y,y^2, z, z^2)
+ "\n";
}
- (*fDebugStream) << "dEdxT" <<
+ if (cstream) (*cstream) << "dEdxT" <<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
"P=" << momenta << // track momenta
"npoints="<<inonEdge<< // number of points
"sector="<<lastSector<< // sector number
if (padType==1) fitterT = fFitter1T;
if (padType==2) fitterT = fFitter2T;
fitterT->AddPoint(xxx,dedxQ[1]);
+ //
+ TLinearFitter *dfitterM = fDFitter0M;
+ if (padType==1) dfitterM=fDFitter1M;
+ if (padType==2) dfitterM=fDFitter2M;
+ dfitterM->AddPoint(xxx,dedxM[1]);
+ //
+ TLinearFitter *dfitterT = fDFitter0T;
+ if (padType==1) dfitterT = fDFitter1T;
+ if (padType==2) dfitterT = fDFitter2T;
+ dfitterT->AddPoint(xxx,dedxQ[1]);
}
// void SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm);
TVectorD vec;
+
//
- fFitter0T->GetParameters(vec);
+ fDFitter0T->Eval();
+ fDFitter1T->Eval();
+ fDFitter2T->Eval();
+ fDFitter0M->Eval();
+ fDFitter1M->Eval();
+ fDFitter2M->Eval();
+ fDFitter0T->GetParameters(vec);
clparam->SetQnorm(0,0,&vec);
- fFitter1T->GetParameters(vec);
+ fDFitter1T->GetParameters(vec);
clparam->SetQnorm(1,0,&vec);
- fFitter2T->GetParameters(vec);
+ fDFitter2T->GetParameters(vec);
clparam->SetQnorm(2,0,&vec);
//
- fFitter0M->GetParameters(vec);
+ fDFitter0M->GetParameters(vec);
clparam->SetQnorm(0,1,&vec);
- fFitter1M->GetParameters(vec);
+ fDFitter1M->GetParameters(vec);
clparam->SetQnorm(1,1,&vec);
- fFitter2M->GetParameters(vec);
+ fDFitter2M->GetParameters(vec);
clparam->SetQnorm(2,1,&vec);
//
}
+
+
+void AliTPCcalibTracksGain::Analyze(){
+
+ Evaluate();
+
+}
+
+
+
+TVectorD * AliTPCcalibTracksGain::MakeQPosNorm(TTree * chain0, Int_t ipad, Bool_t isMax, Int_t maxPoints, Int_t verbose){
+ //
+ // Input parameters
+ // chain0 - the tree with information -Debug stream
+ // ipad - 0 IROC
+ // - 1 OROC medium
+ // - 2 OROC LONG
+ // isMax - kFALSE - total charge param
+ // kTRUE - Max charge param
+ //
+ // maxPoints - number of points for fit
+ //
+ // verbose -
+ //
+ /* e.g
+ ipad=0
+ isMax=kTRUE;
+ maxPoints=1000000;
+ */
+ // Make Q normalization as function of following parameters
+ // 1 - dp - relative pad position
+ // 2 - dt - relative time position
+ // 3 - di - drift length (norm to 1);
+ // 4 - dq0 - Tot/Max charge
+ // 5 - dq1 - Max/Tot charge
+ // 6 - sy - sigma y - shape
+ // 7 - sz - sigma z - shape
+ //
+ // Coeficient of Taylor expansion fitted
+ // Fit parameters returned as TVectorD
+ // Fit parameters to be used in corresponding correction function
+ // in AliTPCclusterParam
+ //
+ //
+ TStatToolkit toolkit;
+ Double_t chi2;
+ TVectorD fitParam;
+ TMatrixD covMatrix;
+ Int_t npoints;
+ TCut cutA("dedge>3&&fraction2<0.7");
+ chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
+ chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
+ chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
+ chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
+ chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
+ chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
+ chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
+ //
+ TString fstring="";
+ fstring+="dp++"; //1
+ fstring+="dt++"; //2
+ fstring+="dp*dp++"; //3
+ fstring+="dt*dt++"; //4
+ fstring+="dt*dt*dt++"; //5
+ fstring+="dp*dt++"; //6
+ fstring+="dp*dt*dt++"; //7
+ fstring+="(dq0)++"; //8
+ fstring+="(dq1)++"; //9
+ //
+ //
+ fstring+="dp*dp*(di)++"; //10
+ fstring+="dt*dt*(di)++"; //11
+ fstring+="dp*dp*sy++"; //12
+ fstring+="dt*sz++"; //13
+ fstring+="dt*dt*sz++"; //14
+ fstring+="dt*dt*dt*sz++"; //15
+ //
+ fstring+="dp*dp*1*sy*sz++"; //16
+ fstring+="dt*sy*sz++"; //17
+ fstring+="dt*dt*sy*sz++"; //18
+ fstring+="dt*dt*dt*sy*sz++"; //19
+ //
+ fstring+="dp*dp*(dq0)++"; //20
+ fstring+="dt*1*(dq0)++"; //21
+ fstring+="dt*dt*(dq0)++"; //22
+ fstring+="dt*dt*dt*(dq0)++"; //23
+ //
+ fstring+="dp*dp*(dq1)++"; //24
+ fstring+="dt*(dq1)++"; //25
+ fstring+="dt*dt*(dq1)++"; //26
+ fstring+="dt*dt*dt*(dq1)++"; //27
+
+ TString var;
+ if (isMax) var = "Cl.fMax/gain/dedxM.fElements[2]";
+ if (!isMax) var = "Cl.fQ/gain/dedxQ.fElements[2]";
+ TString cutP="IPad==";
+ cutP+=ipad;
+ //
+ TString *strq0 = toolkit.FitPlane(chain0,var.Data(),fstring.Data(), cutP.Data()+cutA, chi2,npoints,fitParam,covMatrix,-1,0,maxPoints);
+ //
+ //
+ if (verbose){
+ printf("Chi2/npoints = %f",TMath::Sqrt(chi2/npoints));
+ printf("\nFit function\n:%s\n",strq0->Data());
+ }
+ TVectorD *vec = new TVectorD(fitParam);
+ return vec;
+}
+
+void AliTPCcalibTracksGain::MakeQPosNormAll(TTree * chain, AliTPCClusterParam * param, Int_t maxPoints, Int_t verbose){
+ //
+ // Fill the content of the of the AliTPCclusterParam
+ // with fitted values of corrections
+ //
+ param->fPosQTnorm[0] = MakeQPosNorm(chain,0,kTRUE,maxPoints,verbose);
+ param->fPosQTnorm[1] = MakeQPosNorm(chain,1,kTRUE,maxPoints,verbose);
+ param->fPosQTnorm[2] = MakeQPosNorm(chain,1,kTRUE,maxPoints,verbose);
+ //
+ param->fPosQMnorm[0] = MakeQPosNorm(chain,0,kFALSE,maxPoints,verbose);
+ param->fPosQMnorm[1] = MakeQPosNorm(chain,1,kFALSE,maxPoints,verbose);
+ param->fPosQMnorm[2] = MakeQPosNorm(chain,2,kFALSE,maxPoints,verbose);
+}
+
+
+
+/*
+
+ Position correction fit:
+ //
+TStatToolkit toolkit;
+Double_t chi2;
+TVectorD fitParam;
+TMatrixD covMatrix;
+Int_t npoints;
+//
+TCut cutA("dedge>3&&fraction2<0.7");
+chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
+chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
+chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
+chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
+chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
+chain0->SetAlias("sy","(0.2/sqrt(0.01^2+Cl.fSigmaY2))");
+chain0->SetAlias("sz","(0.2/sqrt(0.01^2+Cl.fSigmaZ2))");
+
+TString fstring="";
+
+fstring+="dp++"; //1
+fstring+="dt++"; //2
+fstring+="dp*dp++"; //3
+fstring+="dt*dt++"; //4
+fstring+="dt*dt*dt++"; //5
+fstring+="dp*dt++"; //6
+fstring+="dp*dt*dt++"; //7
+fstring+="(dq0)++"; //8
+fstring+="(dq1)++"; //9
+//
+//
+fstring+="dp*dp*(di)++"; //10
+fstring+="dt*dt*(di)++"; //11
+fstring+="dp*dp*sy++"; //12
+fstring+="dt*sz++"; //13
+fstring+="dt*dt*sz++"; //14
+fstring+="dt*dt*dt*sz++"; //15
+//
+fstring+="dp*dp*1*sy*sz++"; //16
+fstring+="dt*sy*sz++"; //17
+fstring+="dt*dt*sy*sz++"; //18
+fstring+="dt*dt*dt*sy*sz++"; //19
+//
+fstring+="dp*dp*(dq0)++"; //20
+fstring+="dt*1*(dq0)++"; //21
+fstring+="dt*dt*(dq0)++"; //22
+fstring+="dt*dt*dt*(dq0)++"; //23
+
+fstring+="dp*dp*(dq1)++"; //24
+fstring+="dt*(dq1)++"; //25
+fstring+="dt*dt*(dq1)++"; //26
+fstring+="dt*dt*dt*(dq1)++"; //27
+
+
+ TString *strq0 = toolkit.FitPlane(chain0,"Cl.fMax/gain/dedxM.fElements[2]",fstring->Data(), "IPad==0"+cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
+ TString *strqt0 = toolkit.FitPlane(chain0,"Cl.fQ/gain/dedxQ.fElements[2]",fstring->Data(), "IPad==0"+cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
+
+ chain0->SetAlias("qcorM0",strq0->Data());
+ chain0->SetAlias("qcorT0",strqt0->Data());
+//chain0->SetAlias("mmqcorM0","min(max(qcorM0,0.75),1.15)");
+ chain0->Draw("(Cl.fMax/gain/dedxM.fElements[2]):min(max(qcorM0,0.75),1.15)","IPad==0"+cutA,"prof",100000)
+
+ fraction05 -
+ sigma 0.2419
+ sigma fit 0.2302
+ sigma fit with shape 0.2257
+ fraction 07
+ qtot sigma 0.322
+ qmax sigma 0.292
+ qmax sigma fit 0.2702
+ qmax sigma fit+ratio 0.2638
+
+*/
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff