* provided "as is" without express or implied warranty. *
**************************************************************************/
+/*
+ Comments to be written here:
+ 1. What do we calibrate.
+ 2. How to interpret results
+ 3. Simple example
+ 4. Analysis using debug streamers.
+
+
+
+ 3.Simple example
+ // To make cosmic scan the user interaction neccessary
+ //
+ .x ~/UliStyle.C
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libTPCcalib");
+ TFile fcalib("CalibObjects.root");
+ TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
+ AliTPCcalibCosmic * cosmic = ( AliTPCcalibCosmic *)array->FindObject("cosmicTPC");
+
+
+
+*/
+
+
+
#include "Riostream.h"
#include "TChain.h"
#include "TTree.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TFile.h"
+#include "TF1.h"
+#include "AliTPCclusterMI.h"
#include "AliTPCseed.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliESDfriend.h"
#include "AliESDInputHandler.h"
+#include "AliAnalysisManager.h"
#include "AliTracker.h"
#include "AliMagFMaps.h"
+#include "AliTPCCalROC.h"
#include "AliLog.h"
AliTPCcalibCosmic::AliTPCcalibCosmic()
:AliTPCcalibBase(),
+ fGainMap(0),
fHistNTracks(0),
fClusters(0),
fModules(0),
fHistPt(0),
- fPtResolution(0),
fDeDx(0),
+ fDeDxMIP(0),
+ fMIPvalue(1),
fCutMaxD(5), // maximal distance in rfi ditection
fCutTheta(0.03), // maximal distan theta
fCutMinDir(-0.99) // direction vector products
{
- AliInfo("Defualt Constructor");
+ AliInfo("Default Constructor");
}
AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
:AliTPCcalibBase(),
+ fGainMap(0),
fHistNTracks(0),
fClusters(0),
fModules(0),
fHistPt(0),
- fPtResolution(0),
fDeDx(0),
+ fDeDxMIP(0),
+ fMIPvalue(1),
fCutMaxD(5), // maximal distance in rfi ditection
fCutTheta(0.03), // maximal distan theta
fCutMinDir(-0.99) // direction vector products
SetTitle(title);
AliMagFMaps * field = new AliMagFMaps("dummy1", "dummy2",0,5,0);
AliTracker::SetFieldMap(field, kTRUE);
- fHistNTracks = new TH1F("ntracks","Number of Tracks per Event",501,-0.5,500.5);
- fClusters = new TH1F("signal","Number of Clusters per track",160,0,160);
- fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-1200,1200,600,-1000,1000);
- fHistPt = new TH1F("Pt","Pt distribution",2000,0,50);
- fPtResolution = new TH1F("PtResolution","Pt resolution",100,-50,50);
- fDeDx = new TH2F("DeDx","dEdx",500,0.01,20.,500,0.,500);
+
+ 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);
+ fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
+ fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
+ 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);
+
AliInfo("Non Default Constructor");
+ //
}
AliTPCcalibCosmic::~AliTPCcalibCosmic(){
Printf("ERROR: ESDfriend not available");
return;
}
- FindPairs(event);
+
+ FindPairs(event); // nearly everything takes place in find pairs...
- if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
+ if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
Int_t ntracks=event->GetNumberOfTracks();
fHistNTracks->Fill(ntracks);
- TObjArray tpcSeeds(ntracks);
if (ntracks==0) return;
- //
- //track loop
- //
- for (Int_t i=0;i<ntracks;++i) {
- AliESDtrack *track = event->GetTrack(i);
- fClusters->Fill(track->GetTPCNcls());
- AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
-
- AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
- TObject *calibObject;
- AliTPCseed *seed = 0;
- for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
- if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
- }
- if (seed) tpcSeeds.AddAt(seed,i);
- if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));
- }
- if (ntracks<2) return;
+}
- // dE/dx,pt and ACORDE study --> studies which need the pair selection
- for (Int_t i=0;i<ntracks;++i) {
- AliESDtrack *track1 = event->GetTrack(i);
-
- Double_t d1[3];
- track1->GetDirection(d1);
-
- for (Int_t j=i+1;j<ntracks;++j) {
- AliESDtrack *track2 = event->GetTrack(j);
- Double_t d2[3];
- track2->GetDirection(d2);
-
- if (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2] < -0.999) {
-
- /*___________________________________ Pt resolution ________________________________________*/
- if (track1->Pt() != 0 && track1->GetTPCNcls() > 80 && track2->GetTPCNcls() > 80) {
- Double_t res = (track1->Pt() - track2->Pt());
- res = res/(2*(track1->Pt() + track2->Pt()));
- fPtResolution->Fill(100*res);
- }
-
- /*_______________________________ Propagation to ACORDE ___________________________________*/
- 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
-
- if (d1[1] > 0 && d2[1] < 0 && track1->GetTPCNcls() > 50) {
- Double_t r[3];
- track1->GetXYZ(r);
- Double_t x,z;
- z = r[2] + (d1[2]/d1[1])*(AcordePlane - r[1]);
- x = r[0] + (d1[0]/d1[1])*(AcordePlane - r[1]);
-
- if (x > roof) {
- x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
- z = z - TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
- }
- if (x < -roof) {
- x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
- z = z - TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
- }
-
- fModules->Fill(z, x);
- }
-
- if (d2[1] > 0 && d1[1] < 0 && track2->GetTPCNcls() > 50) {
- Double_t r[3];
- track2->GetXYZ(r);
- Double_t x,z;
- z = r[2] + (d2[2]/d2[1])*(AcordePlane - r[1]);
- x = r[0] + (d2[0]/d2[1])*(AcordePlane - r[1]);
-
- if (x > roof) {
- x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
- z = z - TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
- }
- if (x < -roof) {
- x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
- z = z - TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
- }
-
- fModules->Fill(z, x);
- }
-
- // AliExternalTrackParam * trackOut = new AliExternalTrackParam(*track2->GetOuterParam());
-// AliTracker::PropagateTrackTo(trackOut,850.,105.658,30);
-// delete trackOut;
-
-
+void AliTPCcalibCosmic::Analyze() {
+
+ fMIPvalue = CalculateMIPvalue(fDeDxMIP);
+
+ return;
+
+}
+
- break;
- }
- }
- }
-
-
-
-
-}
void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
//
// Track0 is choosen in upper TPC part
// Track1 is choosen in lower TPC part
//
- if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
+ if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
Int_t ntracks=event->GetNumberOfTracks();
TObjArray tpcSeeds(ntracks);
//
//track loop
//
- for (Int_t i=0;i<ntracks;++i) {
- AliESDtrack *track = event->GetTrack(i);
- fClusters->Fill(track->GetTPCNcls());
- AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
+ for (Int_t i=0;i<ntracks;++i) {
+ AliESDtrack *track = event->GetTrack(i);
+ fClusters->Fill(track->GetTPCNcls());
+
+ const AliExternalTrackParam * trackIn = track->GetInnerParam();
+ const AliExternalTrackParam * trackOut = track->GetOuterParam();
+ if (!trackIn) continue;
+ if (!trackOut) continue;
AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
TObject *calibObject;
if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
}
if (seed) tpcSeeds.AddAt(seed,i);
- if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));
+
+ 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));
+ //
+ if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
+ //
+ 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 (ntracks<2) return;
//
// Find pairs
if (!track0) continue;
if (!track0->GetOuterParam()) continue;
if (track0->GetOuterParam()->GetAlpha()<0) continue;
- Double_t d1[3];
- track0->GetDirection(d1);
+ Double_t dir0[3];
+ track0->GetDirection(dir0);
for (Int_t j=0;j<ntracks;++j) {
if (i==j) continue;
AliESDtrack *track1 = event->GetTrack(j);
if (!track1->GetOuterParam()) continue;
if (track1->GetOuterParam()->GetAlpha()>0) continue;
//
- Double_t d2[3];
- track1->GetDirection(d2);
- printf("My stream level=%d\n",fStreamLevel);
+ Double_t dir1[3];
+ track1->GetDirection(dir1);
+
AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
if (! seed0) continue;
if (! seed1) continue;
- Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0);
- Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0);
- Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
+ 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 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 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 dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
Float_t d0 = track0->GetLinearD(0,0);
Float_t d1 = track1->GetLinearD(0,0);
//
param1.GetXYZ(xyz1);
Bool_t isPair = IsPair(¶m0,¶m1);
//
+ if (isPair) FillAcordeHist(track0);
+ //
if (fStreamLevel>0){
TTreeSRedirector * cstream = GetDebugStreamer();
printf("My stream=%p\n",(void*)cstream);
+ AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
+ AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
+ AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
+ AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
+ Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
+ 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]);
+ Int_t runNr = event->GetRunNumber();
+ UInt_t time = event->GetTimeStamp();
if (cstream) {
(*cstream) << "Track0" <<
+ "runNr="<<runNr<< // run number
+ "time="<<time<< // time stamp of event
"dir="<<dir<< // direction
- "OK="<<isPair<< // will be accepted
+ "OK="<<isPair<< // will be accepted
"b0="<<b0<< // propagate status
"b1="<<b1<< // propagate status
+ "crossI="<<isCrossI<< // cross inner
+ "crossO="<<isCrossO<< // cross outer
+ //
"Orig0.=" << track0 << // original track 0
"Orig1.=" << track1 << // original track 1
"Tr0.="<<¶m0<< // track propagated to the DCA 0,0
"Tr1.="<<¶m1<< // track propagated to the DCA 0,0
+ "Ip0.="<<ip0<< // inner param - upper
+ "Ip1.="<<ip1<< // inner param - lower
+ "Op0.="<<op0<< // outer param - upper
+ "Op1.="<<op1<< // outer param - lower
+ //
"v00="<<dvertex0[0]<< // distance using kalman
"v01="<<dvertex0[1]<< //
"v10="<<dvertex1[0]<< //
"d0="<<d0<< // linear distance to 0,0
"d1="<<d1<< // linear distance to 0,0
//
- "x00="<<xyz0[0]<<
+ //
+ //
+ "x00="<<xyz0[0]<< // global position close to vertex
"x01="<<xyz0[1]<<
"x02="<<xyz0[2]<<
//
- "x10="<<xyz1[0]<<
+ "x10="<<xyz1[0]<< // global position close to vertex
"x11="<<xyz1[1]<<
"x12="<<xyz1[2]<<
//
- "Seed0.=" << track0 << // original seed 0
- "Seed1.=" << track1 << // original seed 1
- "dedx0="<<dedx0<< // dedx0
- "dedx1="<<dedx1<< // dedx1
+ "alpha0="<<alpha0<<
+ "alpha1="<<alpha1<<
+ "dir00="<<dir0[0]<< // direction upper
+ "dir01="<<dir0[1]<<
+ "dir02="<<dir0[2]<<
+ //
+ "dir10="<<dir1[0]<< // direction lower
+ "dir11="<<dir1[1]<<
+ "dir12="<<dir1[2]<<
+ //
+ //
+ "Seed0.=" << seed0 << // original seed 0
+ "Seed1.=" << seed1 << // original seed 1
+ //
+ "dedx0="<<dedx0<< // dedx0 - all
+ "dedx1="<<dedx1<< // dedx1 - all
+ //
+ "dedx0I="<<dedx0I<< // dedx0 - inner ROC
+ "dedx1I="<<dedx1I<< // dedx1 - inner ROC
+ //
+ "dedx0O="<<dedx0O<< // dedx0 - outer ROC
+ "dedx1O="<<dedx1O<< // dedx1 - outer ROC
"\n";
}
}
}
}
-/*
-void AliTPCcalibCosmic::dEdxCorrection(){
- TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
- TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
- TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
+void AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
+
+ // 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 roof = 210.5; // distance from x =0 to end of magnet roof
+
+ Double_t r[3];
+ upperTrack->GetXYZ(r);
+ 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]);
+
+ 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]);
+ }
+ 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]);
+ }
+
+ fModules->Fill(z, x);
+
}
-*/
+Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
+
+ TIterator* iter = li->MakeIterator();
+ AliTPCcalibCosmic* cal = 0;
-Long64_t AliTPCcalibCosmic::Merge(TCollection */*li*/) {
+ 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());
+ return -1;
+ }
+
+ fHistNTracks->Add(cal->GetHistNTracks());
+ fClusters->Add(cal-> GetHistClusters());
+ fModules->Add(cal->GetHistAcorde());
+ fHistPt->Add(cal->GetHistPt());
+ fDeDx->Add(cal->GetHistDeDx());
+ fDeDxMIP->Add(cal->GetHistMIP());
+
+ }
+
+ return 0;
}
+
Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
//
//
//
return kTRUE;
}
+
+
+
+Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
+
+ TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
+ 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));
+
+ delete funcDoubleGaus;
+
+ return MIPvalue;
+
+}
+
+
+void AliTPCcalibCosmic::CalculateBetheParams(TH2F *hist, Double_t * initialParam) {
+
+ return;
+
+}
+
void AliTPCcalibCosmic::BinLogX(TH1 *h) {
}
+AliExternalTrackParam *AliTPCcalibCosmic::Invert(AliExternalTrackParam *input)
+{
+ //
+ // Invert paramerameter - not covariance yet
+ //
+ AliExternalTrackParam *output = new AliExternalTrackParam(*input);
+ Double_t * param = (Double_t*)output->GetParameter();
+ param[0]*=-1;
+ param[3]*=-1;
+ param[4]*=-1;
+ //
+ return output;
+}
+
+/*
+
+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])<4"); // 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("dtheta","(Tr0.fP[3]+Tr1.fP[3])");
+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.))");
+
+}
+
+
+*/
+
+
+/*
+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+="ca++";
+fstring+="sa++";
+fstring+="ca*mtheta++";
+fstring+="sa*mtheta++";
+fstring+="side++";
+fstring+="side*ca++";
+fstring+="side*sa++";
+fstring+="side*ca*mtheta++";
+fstring+="side*sa*mtheta++";
+
+
+TString *strTheta0 = toolkit.FitPlane(chain,"dtheta",fstring->Data(), "hpt"+cutS, 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 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="abs(norm-1)<0.3"+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());
+
+*/
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff