AliTPCcalibCosmic * cosmic = ( AliTPCcalibCosmic *)array->FindObject("cosmicTPC");
-
+ //
+ //
+ gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
+ gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
+ AliXRDPROOFtoolkit tool;
+ TChain * chainCosmic = tool.MakeChainRandom("cosmic.txt","Track0",0,10000);
+ TChain * chainBudget = tool.MakeChainRandom("cosmic.txt","material",0,1000);
+ TCut cutptV="abs(1/pt0V-1/pt1V)<0.1";
+ TCut cutptI="abs(1/pt0In-1/pt1In)<0.5";
+ TCut cutncl="nclmin>120";
+ TCut cutDz="abs(p0.fP[1])<50";
+ TCut cutDr="abs(p0.fP[0])<50";
+ //
+ chainBudget->Draw(">>listB",cutptV+cutptI+cutncl+cutDr+cutDz,"entryList");
+ TEntryList *elistB = (TEntryList*)gDirectory->Get("listB");
+ chainBudget->SetEntryList(elistB);
+
+ chainBudget->SetAlias("dptrel","(pt0V-pt1V)/((pt0V+pt1V)*0.5)");
+ chainBudget->SetAlias("dptInrel","(pt0In-pt1In)/((pt0In+pt1In)*0.5)");
+ chainBudget->SetAlias("ptcorr","(pt0In-pt0V)/(pt0V)+(pt1V-pt1In)/(pt1In)");
*/
#include "TCanvas.h"
#include "TFile.h"
#include "TF1.h"
+#include "THnSparse.h"
+#include "TDatabasePDG.h"
#include "AliTPCclusterMI.h"
#include "AliTPCseed.h"
#include "AliAnalysisManager.h"
#include "AliTracker.h"
-#include "AliMagFMaps.h"
+#include "AliMagF.h"
#include "AliTPCCalROC.h"
#include "AliLog.h"
#include "AliTPCcalibCosmic.h"
-
#include "TTreeStream.h"
#include "AliTPCTracklet.h"
+//#include "AliESDcosmic.h"
+
ClassImp(AliTPCcalibCosmic)
AliTPCcalibCosmic::AliTPCcalibCosmic()
:AliTPCcalibBase(),
- fGainMap(0),
fHistNTracks(0),
fClusters(0),
fModules(0),
fDeDxMIP(0),
fMIPvalue(1),
fCutMaxD(5), // maximal distance in rfi ditection
+ fCutMaxDz(40), // maximal distance in z ditection
fCutTheta(0.03), // maximal distan theta
fCutMinDir(-0.99) // direction vector products
{
AliInfo("Default Constructor");
+ for (Int_t ihis=0; ihis<6;ihis++){
+ fHistoDelta[ihis]=0;
+ fHistoPull[ihis]=0;
+ }
+ for (Int_t ihis=0; ihis<4;ihis++){
+ fHistodEdxMax[ihis] =0;
+ fHistodEdxTot[ihis] =0;
+ }
}
AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
:AliTPCcalibBase(),
- fGainMap(0),
fHistNTracks(0),
fClusters(0),
fModules(0),
fDeDx(0),
fDeDxMIP(0),
fMIPvalue(1),
- fCutMaxD(5), // maximal distance in rfi ditection
+ fCutMaxD(5), // maximal distance in rfi ditection
+ fCutMaxDz(40), // maximal distance in z ditection
fCutTheta(0.03), // maximal distan theta
fCutMinDir(-0.99) // direction vector products
{
SetName(name);
SetTitle(title);
- AliMagFMaps * field = new AliMagFMaps("dummy1", "dummy2",0,5,0);
- AliTracker::SetFieldMap(field, kTRUE);
fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
BinLogX(fDeDx);
fDeDxMIP = new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
-
+ Init();
AliInfo("Non Default Constructor");
//
}
//
//
//
+ for (Int_t ihis=0; ihis<6;ihis++){
+ delete fHistoDelta[ihis];
+ delete fHistoPull[ihis];
+ }
+ for (Int_t ihis=0; ihis<4;ihis++){
+ delete fHistodEdxTot[ihis];
+ delete fHistodEdxMax[ihis];
+ }
+
+ delete fHistNTracks; // histogram showing number of ESD tracks per event
+ delete fClusters; // histogram showing the number of clusters per track
+ delete fModules; // 2d histogram of tracks which are propagated to the ACORDE scintillator array
+ delete fHistPt; // Pt histogram of reconstructed tracks
+ delete fDeDx; // dEdx spectrum showing the different particle types
+ delete fDeDxMIP; // TPC signal close to the MIP region of muons 0.4 < p < 0.45 GeV
}
+void AliTPCcalibCosmic::Init(){
+ //
+ // init component
+ // Make performance histograms
+ //
+
+ // tracking performance bins
+ // 0 - delta of interest
+ // 1 - min (track0, track1) number of clusters
+ // 2 - R - vertex radius
+ // 3 - P1 - mean z
+ // 4 - P2 - snp(phi) at inner wall of TPC
+ // 5 - P3 - tan(theta) at inner wall of TPC
+ // 6 - P4 - 1/pt mean
+ // 7 - pt - pt mean
+ // 8 - alpha
+
+ Double_t xminTrack[9], xmaxTrack[9];
+ Int_t binsTrack[9];
+ TString axisName[9];
+ //
+ binsTrack[0] =100;
+ axisName[0] ="#Delta";
+ //
+ binsTrack[1] =8;
+ xminTrack[1] =80; xmaxTrack[1]=160;
+ axisName[1] ="N_{cl}";
+ //
+ binsTrack[2] =10;
+ xminTrack[2] =0; xmaxTrack[2]=90; //
+ axisName[2] ="dca_{r} (cm)";
+ //
+ binsTrack[3] =25;
+ xminTrack[3] =-250; xmaxTrack[3]=250; //
+ axisName[3] ="z (cm)";
+ //
+ binsTrack[4] =10;
+ xminTrack[4] =-0.8; xmaxTrack[4]=0.8; //
+ axisName[4] ="sin(#phi)";
+ //
+ binsTrack[5] =10;
+ xminTrack[5] =-1; xmaxTrack[5]=1; //
+ axisName[5] ="tan(#theta)";
+ //
+ binsTrack[6] =40;
+ xminTrack[6] =-2; xmaxTrack[6]=2; //
+ axisName[6] ="1/pt (1/GeV)";
+ //
+ binsTrack[7] =40;
+ xminTrack[7] =0.2; xmaxTrack[7]=50; //
+ axisName[7] ="pt (GeV)";
+ //
+ binsTrack[8] =18;
+ xminTrack[8] =0; xmaxTrack[8]=TMath::Pi(); //
+ axisName[8] ="alpha";
+ //
+ // delta y
+ xminTrack[0] =-1; xmaxTrack[0]=1; //
+ fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0] =-5; xmaxTrack[0]=5; //
+ fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
+ //
+ // delta z
+ xminTrack[0] =-1; xmaxTrack[0]=1; //
+ fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0] =-5; xmaxTrack[0]=5; //
+ fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
+ //
+ // delta P2
+ xminTrack[0] =-10; xmaxTrack[0]=10; //
+ fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", 9, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0] =-5; xmaxTrack[0]=5; //
+ fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
+ //
+ // delta P3
+ xminTrack[0] =-10; xmaxTrack[0]=10; //
+ fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", 9, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0] =-5; xmaxTrack[0]=5; //
+ fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
+ //
+ // delta P4
+ xminTrack[0] =-0.2; xmaxTrack[0]=0.2; //
+ fHistoDelta[4] = new THnSparseS("#Delta_{1/pt} (1/GeV)","#Delta_{1/pt} (1/GeV)", 9, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0] =-5; xmaxTrack[0]=5; //
+ fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
+
+ //
+ // delta Pt
+ xminTrack[0] =-0.5; xmaxTrack[0]=0.5; //
+ fHistoDelta[5] = new THnSparseS("#Delta_{pt}/p_{t}","#Delta_{pt}/p_{t}", 9, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0] =-5; xmaxTrack[0]=5; //
+ fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
+ //
+
+ for (Int_t idedx=0;idedx<4;idedx++){
+ xminTrack[0] =0.5; xmaxTrack[0]=1.5; //
+ binsTrack[1] =40;
+ xminTrack[1] =10; xmaxTrack[1]=160;
+
+ fHistodEdxMax[idedx] = new THnSparseS(Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
+ Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
+ 9, binsTrack,xminTrack, xmaxTrack);
+ fHistodEdxTot[idedx] = new THnSparseS(Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
+ Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
+ 9, binsTrack,xminTrack, xmaxTrack);
+ }
+
+
+
+ for (Int_t ivar=0;ivar<6;ivar++){
+ for (Int_t ivar2=0;ivar2<9;ivar2++){
+ fHistoDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
+ fHistoDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
+ fHistoPull[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
+ fHistoPull[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
+ BinLogX(fHistoDelta[ivar],7);
+ BinLogX(fHistoPull[ivar],7);
+ if (ivar<4){
+ fHistodEdxMax[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
+ fHistodEdxMax[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
+ fHistodEdxTot[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
+ fHistodEdxTot[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
+ BinLogX(fHistodEdxMax[ivar],7);
+ BinLogX(fHistodEdxTot[ivar],7);
+ }
+ }
+ }
+}
+
+
+void AliTPCcalibCosmic::Add(const AliTPCcalibCosmic* cosmic){
+ //
+ //
+ //
+ for (Int_t ivar=0; ivar<6;ivar++){
+ if (fHistoDelta[ivar] && cosmic->fHistoDelta[ivar]){
+ fHistoDelta[ivar]->Add(cosmic->fHistoDelta[ivar]);
+ }
+ if (fHistoPull[ivar] && cosmic->fHistoPull[ivar]){
+ fHistoPull[ivar]->Add(cosmic->fHistoPull[ivar]);
+ }
+ }
+ for (Int_t ivar=0; ivar<4;ivar++){
+ if (fHistodEdxMax[ivar] && cosmic->fHistodEdxMax[ivar]){
+ fHistodEdxMax[ivar]->Add(cosmic->fHistodEdxMax[ivar]);
+ }
+ if (fHistodEdxTot[ivar] && cosmic->fHistodEdxTot[ivar]){
+ fHistodEdxTot[ivar]->Add(cosmic->fHistodEdxTot[ivar]);
+ }
+ }
+}
+
Printf("ERROR: ESD not available");
return;
}
- AliESDfriend *
ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
- if (!ESDfriend) {
- Printf("ERROR: ESDfriend not available");
+ AliESDfriend *
esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
+ if (!esdFriend) {
+ Printf("ERROR: esdFriend not available");
return;
}
Int_t ntracks=event->GetNumberOfTracks();
fHistNTracks->Fill(ntracks);
if (ntracks==0) return;
-
+ // AliESDcosmic cosmicESD;
+// TTreeSRedirector * cstream = GetDebugStreamer();
+// cosmicESD.SetDebugStreamer(cstream);
+// cosmicESD.ProcessEvent(event);
+// if (cstream) cosmicESD.DumpToTree();
+
+
}
+void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, const AliExternalTrackParam *par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam */*inner1*/, AliTPCseed *seed0, AliTPCseed *seed1, const AliExternalTrackParam *param0Combined ){
+ //
+ // par0,par1 - parameter of tracks at DCA 0
+ // inner0,inner1 - parameter of tracks at the TPC entrance
+ // seed0, seed1 - detailed track information
+ // param0Combined - Use combined track parameters for binning
+ //
+ Int_t kMinCldEdx =20;
+ Int_t ncl0 = seed0->GetNumberOfClusters();
+ Int_t ncl1 = seed1->GetNumberOfClusters();
+
+ const Double_t kpullCut = 10;
+ Double_t x[9];
+ Double_t xyz0[3];
+ Double_t xyz1[3];
+ par0->GetXYZ(xyz0);
+ par1->GetXYZ(xyz1);
+ Double_t radius0 = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]);
+ Double_t radius1 = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
+ inner0->GetXYZ(xyz0);
+ Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
+ // bin parameters
+ x[1] = TMath::Min(ncl0,ncl1);
+ x[2] = (radius0+radius1)*0.5;
+ x[3] = param0Combined->GetZ();
+ x[4] = inner0->GetSnp();
+ x[5] = param0Combined->GetTgl();
+ x[6] = param0Combined->GetSigned1Pt();
+ x[7] = param0Combined->Pt();
+ x[8] = alpha;
+ // deltas
+ Double_t delta[6];
+ Double_t sigma[6];
+ delta[0] = (par0->GetY()+par1->GetY());
+ delta[1] = (par0->GetZ()-par1->GetZ());
+ delta[2] = (par0->GetAlpha()-par1->GetAlpha()-TMath::Pi());
+ delta[3] = (par0->GetTgl()+par1->GetTgl());
+ delta[4] = (par0->GetParameter()[4]+par1->GetParameter()[4]);
+ delta[5] = (par0->Pt()-par1->Pt())/((par0->Pt()+par1->Pt())*0.5);
+ //
+ sigma[0] = TMath::Sqrt(par0->GetSigmaY2()+par1->GetSigmaY2());
+ sigma[1] = TMath::Sqrt(par0->GetSigmaZ2()+par1->GetSigmaZ2());
+ sigma[2] = TMath::Sqrt(par0->GetSigmaSnp2()+par1->GetSigmaSnp2());
+ sigma[3] = TMath::Sqrt(par0->GetSigmaTgl2()+par1->GetSigmaTgl2());
+ sigma[4] = TMath::Sqrt(par0->GetSigma1Pt2()+par1->GetSigma1Pt2());
+ sigma[5] = sigma[4]*((par0->Pt()+par1->Pt())*0.5);
+ //
+ Bool_t isOK = kTRUE;
+ for (Int_t ivar=0;ivar<6;ivar++){
+ if (sigma[ivar]==0) isOK=kFALSE;
+ x[0]= delta[ivar]/sigma[ivar];
+ if (TMath::Abs(x[0])>kpullCut) isOK = kFALSE;
+ }
+ //
+
+ if (isOK) for (Int_t ivar=0;ivar<6;ivar++){
+ x[0]= delta[ivar]/TMath::Sqrt(2);
+ if (ivar==2 || ivar ==3) x[0]*=1000;
+ fHistoDelta[ivar]->Fill(x);
+ if (sigma[ivar]>0){
+ x[0]= delta[ivar]/sigma[ivar];
+ fHistoPull[ivar]->Fill(x);
+ }
+ }
+
+ //
+ // Fill dedx performance
+ //
+ for (Int_t ipad=0; ipad<4;ipad++){
+ //
+ //
+ //
+ Int_t row0=0;
+ Int_t row1=160;
+ if (ipad==0) row1=63;
+ if (ipad==1) {row0=63; row1=63+64;}
+ if (ipad==2) {row0=128;}
+ Int_t nclUp = TMath::Nint(seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
+ Int_t nclDown = TMath::Nint(seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
+ Int_t minCl = TMath::Min(nclUp,nclDown);
+ if (minCl<kMinCldEdx) continue;
+ x[1] = minCl;
+ //
+ Float_t dEdxTotUp = seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
+ Float_t dEdxTotDown = seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
+ Float_t dEdxMaxUp = seed0->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
+ Float_t dEdxMaxDown = seed1->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
+ //
+ if (dEdxTotDown<=0) continue;
+ if (dEdxMaxDown<=0) continue;
+ x[0]=dEdxTotUp/dEdxTotDown;
+ fHistodEdxTot[ipad]->Fill(x);
+ x[0]=dEdxMaxUp/dEdxMaxDown;
+ fHistodEdxMax[ipad]->Fill(x);
+ }
+
+
+
+}
+
+void AliTPCcalibCosmic::MaterialBudgetDump(AliExternalTrackParam *const par0, AliExternalTrackParam *const par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam *inner1, AliTPCseed *const seed0, AliTPCseed *const seed1, AliExternalTrackParam *const param0Combined, AliExternalTrackParam *const param1Combined){
+ //
+ // matrial budget AOD dump
+ //
+ // par0,par1 - parameter of tracks at DCA 0
+ // inner0,inner1 - parameter of tracks at the TPC entrance
+ // seed0, seed1 - detailed track information
+ // param0Combined - Use combined track parameters for binning
+ // param1Combined -
+ Double_t p0In = inner0->GetP();
+ Double_t p1In = inner1->GetP();
+ Double_t p0V = par0->GetP();
+ Double_t p1V = par1->GetP();
+ //
+ Double_t pt0In = inner0->Pt();
+ Double_t pt1In = inner1->Pt();
+ Double_t pt0V = par0->Pt();
+ Double_t pt1V = par1->Pt();
+ Int_t ncl0 = seed0->GetNumberOfClusters();
+ Int_t ncl1 = seed1->GetNumberOfClusters();
+ Int_t nclmin=TMath::Min(ncl0,ncl1);
+ Double_t sign = (param0Combined->GetSigned1Pt()>0) ? 1:-1.;
+ //
+ TTreeSRedirector * pcstream = GetDebugStreamer();
+ if (pcstream){
+ (*pcstream)<<"material"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "triggerClass="<<&fTriggerClass<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ "sign="<<sign<< // sign of the track
+ //
+ "ncl0="<<ncl0<<
+ "ncl1="<<ncl1<<
+ "nclmin="<<nclmin<<
+ //
+ "p0In="<<p0In<<
+ "p1In="<<p1In<<
+ "p0V="<<p0V<<
+ "p1V="<<p1V<<
+ "pt0In="<<pt0In<<
+ "pt1In="<<pt1In<<
+ "pt0V="<<pt0V<<
+ "pt1V="<<pt1V<<
+ "p0.="<<par0<<
+ "p1.="<<par1<<
+ "up0.="<<param0Combined<<
+ "up1.="<<param1Combined<<
+ "\n";
+ }
+
+}
void AliTPCcalibCosmic::Analyze() {
// Track1 is choosen in lower TPC part
//
if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
- AliESDfriend *
ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
+ AliESDfriend *
esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
Int_t ntracks=event->GetNumberOfTracks();
TObjArray tpcSeeds(ntracks);
if (ntracks==0) return;
const AliExternalTrackParam * trackOut = track->GetOuterParam();
if (!trackIn) continue;
if (!trackOut) continue;
-
- AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
+ if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // filter laser
+
+
+ AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
TObject *calibObject;
AliTPCseed *seed = 0;
for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
- fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
+ fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159));
//
- if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
+ if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159));
//
- if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap)>300) {
- TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
- if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
- if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
- }
+ // if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159)>300) {
+// //TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
+// //if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
+// // if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
+// }
}
AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
if (! seed0) continue;
if (! seed1) continue;
- Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
- Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
+ Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159);
+ Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159);
//
- Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
- Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
+ Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63);
+ Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63);
//
- Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
- Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
+ Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159);
+ Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159);
//
Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
Float_t d0 = track0->GetLinearD(0,0);
//
// Propagate using Magnetic field and correct fo material budget
//
- AliTracker::PropagateTrackTo(¶m0,dmax+1,0.0005,3,kTRUE);
- AliTracker::PropagateTrackTo(¶m1,dmax+1,0.0005,3,kTRUE);
+ Double_t sign0=-1;
+ Double_t sign1=1;
+ Double_t maxsnp=0.90;
+ AliTracker::PropagateTrackToBxByBz(¶m0,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign0);
+ AliTracker::PropagateTrackToBxByBz(¶m1,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign1);
//
// Propagate rest to the 0,0 DCA - z should be ignored
//
//
param0.GetDZ(0,0,0,bz,dvertex0);
param1.GetDZ(0,0,0,bz,dvertex1);
+ if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
//
Double_t xyz0[3];//,pxyz0[3];
Double_t xyz1[3];//,pxyz1[3];
Bool_t isPair = IsPair(¶m0,¶m1);
//
if (isPair) FillAcordeHist(track0);
+ //
+ // combined track params
+ //
+ AliExternalTrackParam *par0U=MakeCombinedTrack(¶m0,¶m1);
+ AliExternalTrackParam *par1U=MakeCombinedTrack(¶m1,¶m0);
+
+
//
if (fStreamLevel>0){
TTreeSRedirector * cstream = GetDebugStreamer();
Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
+ //
+ //
+ //
+ FillHistoPerformance(¶m0, ¶m1, ip0, ip1, seed0, seed1,par0U);
+ MaterialBudgetDump(¶m0, ¶m1, ip0, ip1, seed0, seed1,par0U,par1U);
if (cstream) {
(*cstream) << "Track0" <<
"run="<<fRun<< // run number
"event="<<fEvent<< // event number
"time="<<fTime<< // time stamp of event
"trigger="<<fTrigger<< // trigger
+ "triggerClass="<<&fTriggerClass<< // trigger
"mag="<<fMagF<< // magnetic field
"dir="<<dir<< // direction
"OK="<<isPair<< // will be accepted
"Ip1.="<<ip1<< // inner param - lower
"Op0.="<<op0<< // outer param - upper
"Op1.="<<op1<< // outer param - lower
+ "Up0.="<<par0U<< // combined track 0
+ "Up1.="<<par1U<< // combined track 1
//
"v00="<<dvertex0[0]<< // distance using kalman
"v01="<<dvertex0[1]<< //
"dedx1O="<<dedx1O<< // dedx1 - outer ROC
"\n";
}
- }
+ }
+ delete par0U;
+ delete par1U;
}
}
}
// Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
- const Double_t AcordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
+ const Double_t acordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
Double_t r[3];
Double_t d[3];
upperTrack->GetDirection(d);
Double_t x,z;
- z = r[2] + (d[2]/d[1])*(AcordePlane - r[1]);
- x = r[0] + (d[0]/d[1])*(AcordePlane - r[1]);
+ z = r[2] + (d[2]/d[1])*(acordePlane - r[1]);
+ x = r[0] + (d[0]/d[1])*(acordePlane - r[1]);
if (x > roof) {
- x = r[0] + (d[0]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
- z = r[2] + (d[2]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
+ x = r[0] + (d[0]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
+ z = r[2] + (d[2]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
}
if (x < -roof) {
- x = r[0] + (d[0]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
- z = r[2] + (d[2]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
+ x = r[0] + (d[0]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
+ z = r[2] + (d[2]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
}
fModules->Fill(z, x);
-Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
+Long64_t AliTPCcalibCosmic::Merge(TCollection *const li) {
TIterator* iter = li->MakeIterator();
AliTPCcalibCosmic* cal = 0;
while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
if (!cal->InheritsFrom(AliTPCcalibCosmic::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;
}
fHistPt->Add(cal->GetHistPt());
fDeDx->Add(cal->GetHistDeDx());
fDeDxMIP->Add(cal->GetHistMIP());
-
+ Add(cal);
}
-
return 0;
}
const Double_t *p0 = tr0->GetParameter();
const Double_t *p1 = tr1->GetParameter();
if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
+ if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
+
Double_t d0[3], d1[3];
tr0->GetDirection(d0);
tr1->GetDirection(d1);
funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
hist->Fit(funcDoubleGaus);
- Double_t MIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
+ Double_t aMIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
delete funcDoubleGaus;
- return MIPvalue;
+ return aMIPvalue;
}
-void AliTPCcalibCosmic::CalculateBetheParams(TH2F *hist, Double_t * initialParam) {
-
+void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) {
+ //
+ // Not implemented yet
+ //
return;
}
-void AliTPCcalibCosmic::BinLogX(TH1 *h) {
+void AliTPCcalibCosmic::BinLogX(THnSparse *const h, Int_t axisDim) {
+
+ // Method for the correct logarithmic binning of histograms
+
+ TAxis *axis = h->GetAxis(axisDim);
+ int bins = axis->GetNbins();
+
+ Double_t from = axis->GetXmin();
+ Double_t to = axis->GetXmax();
+ Double_t *newBins = new Double_t[bins + 1];
+
+ newBins[0] = from;
+ Double_t factor = pow(to/from, 1./bins);
+
+ for (int i = 1; i <= bins; i++) {
+ newBins[i] = factor * newBins[i-1];
+ }
+ axis->Set(bins, newBins);
+ delete newBins;
+
+}
+
+
+void AliTPCcalibCosmic::BinLogX(TH1 *const h) {
// Method for the correct logarithmic binning of histograms
Double_t from = axis->GetXmin();
Double_t to = axis->GetXmax();
- Double_t *new_bins = new Double_t[bins + 1];
+ Double_t *newBins = new Double_t[bins + 1];
- new_bins[0] = from;
+ newBins[0] = from;
Double_t factor = pow(to/from, 1./bins);
for (int i = 1; i <= bins; i++) {
- new_bins[i] = factor * new_bins[i-1];
+ newBins[i] = factor * newBins[i-1];
}
- axis->Set(bins, new_bins);
- delete new_bins;
+ axis->Set(bins, newBins);
+ delete newBins;
}
-AliExternalTrackParam *AliTPCcalibCosmic::Invert(AliExternalTrackParam *input)
-{
+
+AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
//
- // Invert paramerameter - not covariance yet
+ //
+ //
+ AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
+ par1R->Rotate(track0->GetAlpha());
+ par1R->PropagateTo(track0->GetX(),AliTracker::GetBz());
//
- AliExternalTrackParam *output = new AliExternalTrackParam(*input);
- Double_t * param = (Double_t*)output->GetParameter();
- param[0]*=-1;
- param[3]*=-1;
- param[4]*=-1;
//
- return output;
+ Double_t * param = (Double_t*)par1R->GetParameter();
+ Double_t * covar = (Double_t*)par1R->GetCovariance();
+
+ param[0]*=1; //OK
+ param[1]*=1; //OK
+ param[2]*=1; //?
+ param[3]*=-1; //OK
+ param[4]*=-1; //OK
+ //
+ covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
+ //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
+ covar[13]*=-1.;
+ return par1R;
}
-/*
-
-void Init(){
-
-.x ~/UliStyle.C
-.x ~/rootlogon.C
-gSystem->Load("libSTAT.so");
-gSystem->Load("libANALYSIS");
-gSystem->Load("libTPCcalib");
-gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
-
-gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
-AliXRDPROOFtoolkit tool;
-TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
-chain->Lookup();
-
-TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01"); // OK
-TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2"); // OK
-TCut cutP1("cutP1","abs(Tr0.fP[1]-Tr1.fP[1])<10"); // OK
-TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
-TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
-TCut cutA=cutT+cutD+cutPt+cutN+cutP1;
-
-TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
-
-//
-chain->Draw(">>listEL",cutA,"entryList");
-TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
-chain->SetEntryList(elist);
-//
-chain->Draw(">>listV40Z100","abs(d0)<40&&abs(v01)<100","entryList");
-TEntryList *elistV40Z100 = (TEntryList*)gDirectory->Get("listV40Z100");
-chain->SetEntryList(elistV40Z100);
-
-//
-// Aliases
-//
-chain->SetAlias("side","(-1+(Tr0.fP[1]>0)*2)");
-chain->SetAlias("hpt","abs(Tr0.fP[4])<0.2");
-chain->SetAlias("signy","(-1+(Tr0.fP[0]>0)*2)");
-
-chain->SetAlias("dy","Tr0.fP[0]+Tr1.fP[0]");
-chain->SetAlias("dz","Tr0.fP[1]-Tr1.fP[1]");
-chain->SetAlias("d1pt","Tr0.fP[4]+Tr1.fP[4]");
-chain->SetAlias("dtheta","(Tr0.fP[3]+Tr1.fP[3])");
-chain->SetAlias("dphi","(Tr0.fAlpha-Tr1.fAlpha-pi)");
-
-chain->SetAlias("mtheta","(Tr0.fP[3]-Tr1.fP[3])*0.5")
-chain->SetAlias("sa","(sin(Tr0.fAlpha+0.))");
-chain->SetAlias("ca","(cos(Tr0.fAlpha+0.))");
-
-
-
-chain->Draw("dy:sqrt(abs(Tr0.fP[4]))>>hisdyA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
-hisdyA->FitSlicesY();
-hisdyA_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdyA_2->SetYTitle("#sigma_{y}(cm)");
-hisdyA_2->SetTitle("Cosmic - Y matching");
-hisdyA_2->SetMaximum(0.5);
-
-
-chain->Draw("dy:sqrt(abs(Tr0.fP[4]))>>hisdyC(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
-hisdyC->FitSlicesY();
-hisdyC_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdyC_2->SetYTitle("#sigma_{y}(cm)");
-hisdyC_2->SetTitle("Cosmic - Y matching");
-hisdyC_2->SetMaximum(1);
-hisdyC_2->SetMinimum(0);
-hisdyC_2->SetMarkerStyle(22);
-hisdyA_2->SetMarkerStyle(21);
-hisdyC_2->SetMarkerSize(1.5);
-hisdzA_2->SetMarkerSize(1.5);
-hisdyC_2->Draw();
-hisdyA_2->Draw("same");
-gPad->SaveAs("~/Calibration/Cosmic/pic/ymatching.gif")
-
-chain->Draw("dz:sqrt(abs(Tr0.fP[4]))>>hisdzA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
-hisdzA->FitSlicesY();
-hisdzA_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdzA_2->SetYTitle("#sigma_{z}(cm)");
-hisdzA_2->SetTitle("Cosmic - Z matching - A side ");
-hisdzA_2->SetMaximum(0.5);
-
-chain->Draw("dz:sqrt(abs(Tr0.fP[4]))>>hisdzC(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
-hisdzC->FitSlicesY();
-hisdzC_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdzC_2->SetYTitle("#sigma_{z}(cm)");
-hisdzC_2->SetTitle("Cosmic - Z matching");
-hisdzC_2->SetMaximum(0.5);
-hisdzC_2->SetMarkerStyle(22);
-hisdzA_2->SetMarkerStyle(21);
-hisdzC_2->SetMarkerSize(1.5);
-hisdzA_2->SetMarkerSize(1.5);
-
-hisdzC_2->Draw();
-hisdzA_2->Draw("same");
-
-
-//
-// PICTURE 1/pt
-//
-chain->Draw("d1pt:sqrt(abs(Tr0.fP[4]))>>hisd1ptA(5,0,1,30,-0.1,0.1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
-hisd1ptA->FitSlicesY();
-hisd1ptA_2->SetXTitle("#sqrt{1/p_{t}}");
-hisd1ptA_2->SetYTitle("#sigma_{z}(cm)");
-hisd1ptA_2->SetTitle("Cosmic - Z matching - A side ");
-hisd1ptA_2->SetMaximum(0.5);
-
-chain->Draw("d1pt:sqrt(abs(Tr0.fP[4]))>>hisd1ptC(5,0,1,30,-0.1,0.1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
-hisd1ptC->FitSlicesY();
-hisd1ptC_2->SetXTitle("#sqrt{1/p_{t}}");
-hisd1ptC_2->SetYTitle("#sigma_{1/pt}(1/GeV)");
-hisd1ptC_2->SetTitle("Cosmic - 1/pt matching");
-hisd1ptC_2->SetMaximum(0.05);
-hisd1ptC_2->SetMarkerStyle(22);
-hisd1ptA_2->SetMarkerStyle(21);
-hisd1ptC_2->SetMarkerSize(1.5);
-hisd1ptA_2->SetMarkerSize(1.5);
-
-hisd1ptC_2->Draw();
-hisd1ptA_2->Draw("same");
-gPad->SaveAs("~/Calibration/Cosmic/pic/1ptmatching.gif")
-
-//
-// Theta
-//
-chain->Draw("dtheta:sqrt(abs(Tr0.fP[4]))>>hisdthetaA(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
-hisdthetaA->FitSlicesY();
-hisdthetaA_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdthetaA_2->SetYTitle("#sigma_{#theta}(cm)");
-hisdthetaA_2->SetTitle("Cosmic - Z matching - A side ");
-hisdthetaA_2->SetMaximum(0.5);
-
-chain->Draw("dtheta:sqrt(abs(Tr0.fP[4]))>>hisdthetaC(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
-hisdthetaC->FitSlicesY();
-hisdthetaC_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdthetaC_2->SetYTitle("#sigma_{#theta}(rad)");
-hisdthetaC_2->SetTitle("Cosmic - Theta matching");
-hisdthetaC_2->SetMaximum(0.01);
-hisdthetaC_2->SetMinimum(0.0);
-hisdthetaC_2->SetMarkerStyle(22);
-hisdthetaA_2->SetMarkerStyle(21);
-hisdthetaC_2->SetMarkerSize(1.5);
-hisdthetaA_2->SetMarkerSize(1.5);
-
-hisdthetaC_2->Draw();
-hisdthetaA_2->Draw("same");
-gPad->SaveAs("~/Calibration/Cosmic/pic/thetamatching.gif")
-//
-// Phi
-//
-chain->Draw("dphi:sqrt(abs(Tr0.fP[4]))>>hisdphiA(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
-hisdphiA->FitSlicesY();
-hisdphiA_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdphiA_2->SetYTitle("#sigma_{#phi}(rad)");
-hisdphiA_2->SetTitle("Cosmic - Z matching - A side ");
-hisdphiA_2->SetMaximum(0.5);
-
-chain->Draw("dphi:sqrt(abs(Tr0.fP[4]))>>hisdphiC(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
-hisdphiC->FitSlicesY();
-hisdphiC_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdphiC_2->SetYTitle("#sigma_{#phi}(rad)");
-hisdphiC_2->SetTitle("Cosmic - Phi matching");
-hisdphiC_2->SetMaximum(0.01);
-hisdphiC_2->SetMinimum(0.0);
-hisdphiC_2->SetMarkerStyle(22);
-hisdphiA_2->SetMarkerStyle(21);
-hisdphiC_2->SetMarkerSize(1.5);
-hisdphiA_2->SetMarkerSize(1.5);
-
-hisdphiC_2->Draw();
-hisdphiA_2->Draw("same");
-gPad->SaveAs("~/Calibration/Cosmic/pic/phimatching.gif")
-
-
-
+AliExternalTrackParam *AliTPCcalibCosmic::MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
+ //
+ // Make combined track
+ //
+ //
+ AliExternalTrackParam * par1T = MakeTrack(track0,track1);
+ AliExternalTrackParam * par0U = new AliExternalTrackParam(*track0);
+ //
+ UpdateTrack(*par0U,*par1T);
+ delete par1T;
+ return par0U;
}
-*/
-
-
-/*
-void MatchTheta(){
-
-TStatToolkit toolkit;
-Double_t chi2=0;
-Int_t npoints=0;
-TVectorD fitParamA0;
-TVectorD fitParamA1;
-TVectorD fitParamC0;
-TVectorD fitParamC1;
-TMatrixD covMatrix;
-
-
-TString fstring="";
-//
-fstring+="mtheta++";
-fstring+="ca++";
-fstring+="sa++";
-fstring+="ca*mtheta++";
-fstring+="sa*mtheta++";
-//
-fstring+="side++";
-fstring+="side*mtheta++";
-fstring+="side*ca++";
-fstring+="side*sa++";
-fstring+="side*ca*mtheta++";
-fstring+="side*sa*mtheta++";
-
-
-TString *strTheta0 = toolkit.FitPlane(chain,"dtheta",fstring->Data(), "hpt&&!crossI&&!crossO", chi2,npoints,fitParamA0,covMatrix,0.8);
-chain->SetAlias("dtheta0",strTheta0.Data());
-strTheta0->Tokenize("+")->Print();
-
-
-//fstring+="mtheta++";
-//fstring+="mtheta^2++";
-//fstring+="ca*mtheta^2++";
-//fstring+="sa*mtheta^2++";
-
-
+void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
+ //
+ // Update track 1 with track 2
+ //
+ //
+ //
+ TMatrixD vecXk(5,1); // X vector
+ TMatrixD covXk(5,5); // X covariance
+ TMatrixD matHk(5,5); // vector to mesurement
+ TMatrixD measR(5,5); // measurement error
+ TMatrixD vecZk(5,1); // measurement
+ //
+ TMatrixD vecYk(5,1); // Innovation or measurement residual
+ TMatrixD matHkT(5,5);
+ TMatrixD matSk(5,5); // Innovation (or residual) covariance
+ TMatrixD matKk(5,5); // Optimal Kalman gain
+ TMatrixD mat1(5,5); // update covariance matrix
+ TMatrixD covXk2(5,5); //
+ TMatrixD covOut(5,5);
+ //
+ Double_t *param1=(Double_t*) track1.GetParameter();
+ Double_t *covar1=(Double_t*) track1.GetCovariance();
+ Double_t *param2=(Double_t*) track2.GetParameter();
+ Double_t *covar2=(Double_t*) track2.GetCovariance();
+ //
+ // copy data to the matrix
+ for (Int_t ipar=0; ipar<5; ipar++){
+ for (Int_t jpar=0; jpar<5; jpar++){
+ covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
+ measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
+ matHk(ipar,jpar)=0;
+ mat1(ipar,jpar)=0;
+ }
+ vecXk(ipar,0) = param1[ipar];
+ vecZk(ipar,0) = param2[ipar];
+ matHk(ipar,ipar)=1;
+ mat1(ipar,ipar)=0;
+ }
+ //
+ //
+ //
+ //
+ //
+ vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
+ matHkT=matHk.T(); matHk.T();
+ matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
+ matSk.Invert();
+ matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
+ vecXk += matKk*vecYk; // updated vector
+ covXk2 = (mat1-(matKk*matHk));
+ covOut = covXk2*covXk;
+ //
+ //
+ //
+ // copy from matrix to parameters
+ if (0) {
+ vecXk.Print();
+ vecZk.Print();
+ //
+ measR.Print();
+ covXk.Print();
+ covOut.Print();
+ //
+ track1.Print();
+ track2.Print();
+ }
+ for (Int_t ipar=0; ipar<5; ipar++){
+ param1[ipar]= vecXk(ipar,0) ;
+ for (Int_t jpar=0; jpar<5; jpar++){
+ covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
+ }
+ }
}
-*/
-
-
-
-
-/*
- void PosCorrection()
-
-
-
-
- TStatToolkit toolkit;
- Double_t chi2=0;
- Int_t npoints=0;
- TVectorD fitParam;
- TMatrixD covMatrix;
-
- //Theta
-chain->SetAlias("dthe","(Tr0.fP[3]+Tr1.fP[3])");
-chain->SetAlias("sign","(-1+(Tr0.fP[1]>0)*2)");
-chain->SetAlias("di","(1-abs(Tr0.fP[1])/250)");
-//
-//
-TString strFit="";
-//
-strFit+="sign++"; //1
-strFit+="Tr0.fP[3]++"; //2
-//
-strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]>0)++"; //3
-strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]<0)++"; //4
-//
-strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]>0)++"; //5
-strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]<0)++"; //6
-//
-strFit+="sin(Tr0.fAlpha)*Tr0.fP[3]++"; //7
-strFit+="cos(Tr0.fAlpha)*Tr0.fP[3]++"; //8
-
-
- //
- TString * thetaParam = toolkit.FitPlane(chain,"dthe", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix,0.8,0,10000)
-
- chain->SetAlias("corTheta",thetaParam->Data());
- chain->Draw("dthe:Tr0.fP[1]","","",50000);
-
-
-
-*/
-
-
-
-/*
-
-void AliTPCcalibCosmic::dEdxCorrection(){
- TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01"); // OK
- TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2"); // OK
- TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
- TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
- TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
- TCut cutA=cutT+cutD+cutPt+cutN+cutS;
-
-
- .x ~/UliStyle.C
- gSystem->Load("libANALYSIS");
- gSystem->Load("libTPCcalib");
- gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
- gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
- AliXRDPROOFtoolkit tool;
- TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
- chain->Lookup();
-
- chain->Draw(">>listEL",cutA,"entryList");
- TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
- chain->SetEntryList(elist);
-
- .x ~/rootlogon.C
- gSystem->Load("libSTAT.so");
- TStatToolkit toolkit;
- Double_t chi2=0;
- Int_t npoints=0;
- TVectorD fitParam;
- TMatrixD covMatrix;
-
- chain->Draw("Tr0.fP[4]+Tr1.fP[4]","OK"+cutA);
-
- TString strFit;
- strFit+="(Tr0.fP[1]/250)++";
- strFit+="(Tr0.fP[1]/250)^2++";
- strFit+="(Tr0.fP[3])++";
- strFit+="(Tr0.fP[3])^2++";
-
- TString * ptParam = TStatToolkit::FitPlane(chain,"Tr0.fP[4]+Tr1.fP[4]", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix)
-
-
-
-*/
-
-
-/*
-gSystem->Load("libANALYSIS");
-gSystem->Load("libSTAT");
-gSystem->Load("libTPCcalib");
-
-TStatToolkit toolkit;
-Double_t chi2;
-TVectorD fitParam;
-TMatrixD covMatrix;
-Int_t npoints;
-//
-TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03"); // OK
-TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5"); // OK
-TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.2&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
-TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>110");
-TCut cutA=cutT+cutD+cutPt+cutN;
-
-
-
-TTree * chain = Track0;
-
-
-chain->SetAlias("norm","signalTot0.fElements[3]/signalTot1.fElements[3]");
-//
-chain->SetAlias("dr1","(signalTot0.fElements[1]/signalTot0.fElements[3])");
-chain->SetAlias("dr2","(signalTot0.fElements[2]/signalTot0.fElements[3])");
-chain->SetAlias("dr4","(signalTot0.fElements[4]/signalTot0.fElements[3])");
-chain->SetAlias("dr5","(signalTot0.fElements[5]/signalTot0.fElements[3])");
-
-TString fstring="";
-fstring+="dr1++";
-fstring+="dr2++";
-fstring+="dr4++";
-fstring+="dr5++";
-//
-fstring+="dr1*dr2++";
-fstring+="dr1*dr4++";
-fstring+="dr1*dr5++";
-fstring+="dr2*dr4++";
-fstring+="dr2*dr5++";
-fstring+="dr4*dr5++";
-
-
-
-TString *strqdedx = toolkit.FitPlane(chain,"norm",fstring->Data(), cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
-
-chain->SetAlias("corQT",strqdedx->Data());
-
-*/
-/*
- chain->SetProof(kTRUE);
- chain->Draw("Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE):Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)",""+cutA,"",1000);
-
-
-chain->Draw("Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE)/Seed1.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE)>>his(100,0.5,1.5)","min(Orig0.fTPCncls,Orig1.fTPCncls)>130"+cutA,"",50000);
-
-*/
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff