X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCcalibCosmic.cxx;h=66e960d128b25ef08fb4712aabcb1d12ba47fd2a;hb=b3849ff32a9559a0639dd295e964d37c857f3692;hp=7af7ec4f2c738026352119f3cc3e0c1188f9146d;hpb=29620e4bd38a9cf0db9f86a791bc6043cf629ddc;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCcalibCosmic.cxx b/TPC/AliTPCcalibCosmic.cxx index 7af7ec4f2c7..66e960d128b 100644 --- a/TPC/AliTPCcalibCosmic.cxx +++ b/TPC/AliTPCcalibCosmic.cxx @@ -1,3 +1,5 @@ + + /************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * @@ -13,6 +15,23 @@ * 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 + // + + */ + + + #include "Riostream.h" #include "TChain.h" #include "TTree.h" @@ -22,221 +41,439 @@ #include "TMath.h" #include "TCanvas.h" #include "TFile.h" +#include "TF1.h" +#include "THnSparse.h" +#include "TDatabasePDG.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 "AliMagF.h" +#include "AliTPCCalROC.h" +#include "AliTPCParam.h" #include "AliLog.h" #include "AliTPCcalibCosmic.h" - #include "TTreeStream.h" #include "AliTPCTracklet.h" - +//#include "AliESDcosmic.h" +#include "AliRecoParam.h" +#include "AliMultiplicity.h" +#include "AliTPCTransform.h" +#include "AliTPCcalibDB.h" +#include "AliTPCseed.h" +#include "AliGRPObject.h" +#include "AliTPCCorrection.h" ClassImp(AliTPCcalibCosmic) 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 + fCutMaxDz(40), // maximal distance in z ditection fCutTheta(0.03), // maximal distan theta - fCutMinDir(-0.99) // direction vector products + fCutMinDir(-0.99), // direction vector products + fCosmicTree(0) // tree with cosmic data { - AliInfo("Defualt Constructor"); - TFile f("/u/miranov/calibKr.root"); - AliTPCCalPad *gainMap = (AliTPCCalPad *)f.Get("spectrMean"); - if (gainMap) { - gainMap->Multiply(1/gainMap->GetMedian()); - fGainMap =gainMap; + // + // CONSTRUCTOR - SEE COMMENTS ABOVE + // + 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), fHistPt(0), - fPtResolution(0), fDeDx(0), - fCutMaxD(5), // maximal distance in rfi ditection + 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 + fCutMinDir(-0.99), // direction vector products + fCosmicTree(0) // tree with cosmic data { + // + // cONSTRUCTOR - SEE COMENTS ABOVE + // 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",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); + Init(); AliInfo("Non Default Constructor"); // - TFile f("/u/miranov/calibKr.root"); - AliTPCCalPad *gainMap = (AliTPCCalPad *)f.Get("spectrMean"); - gainMap->Multiply(1/gainMap->GetMedian()); - fGainMap =gainMap; } AliTPCcalibCosmic::~AliTPCcalibCosmic(){ // + // destructor + // + 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 + // 9 - is corss indicator + Int_t ndim=10; + Double_t xminTrack[10], xmaxTrack[10]; + Int_t binsTrack[10]; + TString axisName[10]; + // + 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] =50; + xminTrack[7] =1; xmaxTrack[7]=1000; // + axisName[7] ="pt (GeV)"; + // + binsTrack[8] =18; + xminTrack[8] =0; xmaxTrack[8]=TMath::Pi(); // + axisName[8] ="alpha"; + // + binsTrack[9] =3; + xminTrack[9] =-0.1; xmaxTrack[9]=2.1; // + axisName[9] ="cross"; + // + // delta y + xminTrack[0] =-1; xmaxTrack[0]=1; // + fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", ndim, binsTrack,xminTrack, xmaxTrack); + xminTrack[0] =-5; xmaxTrack[0]=5; // + fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", ndim, binsTrack,xminTrack, xmaxTrack); + // + // delta z + xminTrack[0] =-1; xmaxTrack[0]=1; // + fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", ndim, binsTrack,xminTrack, xmaxTrack); + xminTrack[0] =-5; xmaxTrack[0]=5; // + fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", ndim, binsTrack,xminTrack, xmaxTrack); + // + // delta P2 + xminTrack[0] =-10; xmaxTrack[0]=10; // + fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", ndim, binsTrack,xminTrack, xmaxTrack); + xminTrack[0] =-5; xmaxTrack[0]=5; // + fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", ndim, binsTrack,xminTrack, xmaxTrack); + // + // delta P3 + xminTrack[0] =-10; xmaxTrack[0]=10; // + fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", ndim, binsTrack,xminTrack, xmaxTrack); + xminTrack[0] =-5; xmaxTrack[0]=5; // + fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", ndim, 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)", ndim, binsTrack,xminTrack, xmaxTrack); + xminTrack[0] =-5; xmaxTrack[0]=5; // + fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", ndim, 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}", ndim, binsTrack,xminTrack, xmaxTrack); + xminTrack[0] =-5; xmaxTrack[0]=5; // + fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", ndim, 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), + ndim, binsTrack,xminTrack, xmaxTrack); + fHistodEdxTot[idedx] = new THnSparseS(Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx), + Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx), + ndim, binsTrack,xminTrack, xmaxTrack); + } + + + + for (Int_t ivar=0;ivar<6;ivar++){ + for (Int_t ivar2=0;ivar2GetAxis(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){ + // + // merge the content of the cosmic componentnts + // + 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]); + } + } + if (cosmic->fCosmicTree){ + if (!fCosmicTree) { + fCosmicTree = new TTree("pairs","pairs"); + fCosmicTree->SetDirectory(0); + } + AliTPCcalibCosmic::AddTree(fCosmicTree,cosmic->fCosmicTree); + } +} + void AliTPCcalibCosmic::Process(AliESDEvent *event) { // - // + // Process of the ESD event - fill calibration components // if (!event) { Printf("ERROR: ESD not available"); return; } - AliESDfriend *ESDfriend=static_cast(event->FindListObject("AliESDfriend")); - if (!ESDfriend) { - Printf("ERROR: ESDfriend not available"); + AliESDfriend *esdFriend=static_cast(event->FindListObject("AliESDfriend")); + if (!esdFriend) { + Printf("ERROR: esdFriend not available"); return; } - FindPairs(event); + + // + //Int_t isOK=kTRUE; + // COSMIC not signed properly + // UInt_t specie = event->GetEventSpecie(); // select only cosmic events + //if (specie==AliRecoParam::kCosmic || specie==AliRecoParam::kCalib) { + // isOK = kTRUE; + //} + //if (!isOK) return; + // Work around + FindCosmicPairs(event); + const AliMultiplicity *multiplicity = event->GetMultiplicity(); + Int_t ntracklets = multiplicity->GetNumberOfTracklets(); + if (ntracklets>6) return; // filter out "normal" event with high multiplicity + const TString &trigger = event->GetFiredTriggerClasses(); + if (trigger.Contains("C0OB0")==0) return; + - if (GetDebugLevel()>1) printf("Hallo world: Im here\n"); + FindPairs(event); // nearly everything takes place in find pairs... + + 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; + +} + + +void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, const AliExternalTrackParam *par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam */*inner1*/, AliTPCseed *seed0, AliTPCseed *seed1, const AliExternalTrackParam *param0Combined , Int_t cross){ // - //track loop + // 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 // - for (Int_t i=0;iGetTrack(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(calibObject))) break; - } - if (seed) tpcSeeds.AddAt(seed,i); - if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0,0,159,fGainMap)); + Int_t kMinCldEdx =20; + Int_t ncl0 = seed0->GetNumberOfClusters(); + Int_t ncl1 = seed1->GetNumberOfClusters(); + const Double_t kpullCut = 10; + Double_t x[10]; + 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; + x[9] = cross; + // 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 (ntracks<2) return; - + // - // dE/dx,pt and ACORDE study --> studies which need the pair selection - for (Int_t i=0;iGetTrack(i); - - Double_t d1[3]; - track1->GetDirection(d1); - - for (Int_t j=i+1;jGetTrack(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; - + if (isOK) for (Int_t ivar=0;ivar<6;ivar++){ + x[0]= delta[ivar]; // Modifiation 10.10 use not normalized deltas + if (ivar==2 || ivar ==3) x[0]*=1000; // angles in radian + 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 (minClCookdEdxAnalytical(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); + } - - break; - } - } - } - - - -} +} -void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) { +void AliTPCcalibCosmic::FindPairs(const AliESDEvent *event){ // // Find cosmic pairs // // Track0 is choosen in upper TPC part // Track1 is choosen in lower TPC part // - if (GetDebugLevel()>1) printf("Hallo world: Im here\n"); - AliESDfriend *ESDfriend=static_cast(event->FindListObject("AliESDfriend")); + if (GetDebugLevel()>20) printf("Hallo world: Im here\n"); + AliESDfriend *esdFriend=static_cast(event->FindListObject("AliESDfriend")); Int_t ntracks=event->GetNumberOfTracks(); TObjArray tpcSeeds(ntracks); if (ntracks==0) return; @@ -246,21 +483,42 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) { // //track loop // - for (Int_t i=0;iGetTrack(i); - fClusters->Fill(track->GetTPCNcls()); - if (!track->GetInnerParam()) continue; - AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam()); - - AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i); + for (Int_t i=0;iGetTrack(i); + fClusters->Fill(track->GetTPCNcls()); + + const AliExternalTrackParam * trackIn = track->GetInnerParam(); + const AliExternalTrackParam * trackOut = track->GetOuterParam(); + if (!trackIn) continue; + if (!trackOut) continue; + if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // filter laser + + + AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i); + if (!friendTrack) continue; TObject *calibObject; - AliTPCseed *seed = 0; + AliTPCseed *seed = 0; for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) { if ((seed=dynamic_cast(calibObject))) break; } if (seed) tpcSeeds.AddAt(seed,i); - if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0,0,159,fGainMap)); + + 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)); + // + 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)>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 @@ -271,8 +529,8 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) { 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;jGetTrack(j); @@ -281,23 +539,23 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) { 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,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 = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]); + 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); // @@ -321,8 +579,11 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) { // // 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 // @@ -331,26 +592,71 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) { // 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]; param0.GetXYZ(xyz0); param1.GetXYZ(xyz1); Bool_t isPair = IsPair(¶m0,¶m1); + // + if (isPair) FillAcordeHist(track0); + if (isPair &¶m0.Pt()>1) { + const TString &trigger = event->GetFiredTriggerClasses(); + UInt_t specie = event->GetEventSpecie(); + printf("COSMIC ?\t%s\t%d\t%f\t%f\n", trigger.Data(),specie, param0.GetZ(), param1.GetZ()); + } + // + // combined track params + // + AliExternalTrackParam *par0U=MakeCombinedTrack(¶m0,¶m1); + AliExternalTrackParam *par1U=MakeCombinedTrack(¶m1,¶m0); + // if (fStreamLevel>0){ TTreeSRedirector * cstream = GetDebugStreamer(); - printf("My stream=%p\n",(void*)cstream); + //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 cross =0; // 0 no cross, 2 cross on both sides + if (isCrossI) cross+=1; + if (isCrossO) cross+=1; + FillHistoPerformance(¶m0, ¶m1, ip0, ip1, seed0, seed1,par0U, cross); if (cstream) { (*cstream) << "Track0" << + "run="<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 *const li) { + // + // component merging + // 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; } - fHistNTracks->Add(cal->fHistNTracks); - fClusters->Add(cal->fClusters); - fModules->Add(cal->fModules); - fHistPt->Add(cal->fHistPt); - fPtResolution->Add(cal->fPtResolution); - fDeDx->Add(cal->fDeDx); - + fHistNTracks->Add(cal->GetHistNTracks()); + fClusters->Add(cal-> GetHistClusters()); + fModules->Add(cal->GetHistAcorde()); + fHistPt->Add(cal->GetHistPt()); + fDeDx->Add(cal->GetHistDeDx()); + fDeDxMIP->Add(cal->GetHistMIP()); + Add(cal); } - return 0; } -Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){ + +Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1) const{ // // /* @@ -431,7 +786,9 @@ Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackPa 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); @@ -439,10 +796,62 @@ Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackPa // return kTRUE; } + + + +Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) { + // + // Calculate the MIP value - gaussian fit used + // + 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 aMIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4)); + + delete funcDoubleGaus; + + return aMIPvalue; + +} -void AliTPCcalibCosmic::BinLogX(TH1 *h) { + + +void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) { + // + // Not implemented yet + // + return; + +} + + +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 @@ -451,78 +860,727 @@ void AliTPCcalibCosmic::BinLogX(TH1 *h) { 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::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){ + // + // Make a atrack using the kalman update of track0 and track1 + // + AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1); + par1R->Rotate(track0->GetAlpha()); + par1R->PropagateTo(track0->GetX(),AliTracker::GetBz()); + // + // + 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; +} +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 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(); + } -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])<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; + 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); + } + } +} - .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(); - .x ~/rootlogon.C - gSystem->Load("libSTAT.so"); +void AliTPCcalibCosmic::FindCosmicPairs(const AliESDEvent * event) { + // + // find cosmic pairs trigger by random trigger + // + // + AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD(); + AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC(); + AliESDfriend *esdFriend=static_cast(event->FindListObject("AliESDfriend")); + const Double_t kMinPt=1; + const Double_t kMinPtMax=0.8; + const Double_t kMinNcl=50; + const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1}; + Int_t ntracks=event->GetNumberOfTracks(); + // Float_t dcaTPC[2]={0,0}; + // Float_t covTPC[3]={0,0,0}; - Double_t chi2=0; - Int_t npoints=0; - TVectorD fitParam; - TMatrixD covMatrix; + UInt_t specie = event->GetEventSpecie(); // skip laser events + if (specie==AliRecoParam::kCalib) return; - 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(),cutA, chi2,npoints,fitParam,covMatrix) -strFit+="(Tr0.fP[1]/250)++"; -strFit+="(Tr0.fP[1]/250)^2++"; -strFit+="(Tr0.fP[3])++"; -strFit+="(Tr0.fP[3])^2++"; -strFit+="(Tr0.fP[1]/250)^2*Tr0.fP[3]++"; -strFit+="(Tr0.fP[1]/250)^2*Tr0.fP[3]^2++"; -// + for (Int_t itrack0=0;itrack0GetTrack(itrack0); + if (!track0) continue; + if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue; + + if (TMath::Abs(AliTracker::GetBz())>1&&track0->Pt()GetTPCncls()GetY())GetKinkIndex(0)>0) continue; + const Double_t * par0=track0->GetParameter(); //track param at rhe DCA + //rm primaries + // + //track0->GetImpactParametersTPC(dcaTPC,covTPC); + //if (TMath::Abs(dcaTPC[0])GetInnerParam(); + for (Int_t itrack1=itrack0+1;itrack1GetTrack(itrack1); + if (!track1) continue; + if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue; + if (track1->GetKinkIndex(0)>0) continue; + if (TMath::Abs(AliTracker::GetBz())>1&&track1->Pt()GetTPCncls()1&&TMath::Max(track1->Pt(), track0->Pt())GetY())GetImpactParametersTPC(dcaTPC,covTPC); + // if (TMath::Abs(dcaTPC[0])GetParameter(); //track param at rhe DCA + // + Bool_t isPair=kTRUE; + for (Int_t ipar=0; ipar<5; ipar++){ + if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off + if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE; + } + if (!isPair) continue; + if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE; + //delta with correct sign + /* + TCut cut0="abs(t1.fP[0]+t0.fP[0])<2" + TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02" + TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2" + */ + if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign + if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign + if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign + if (!isPair) continue; + // const AliExternalTrackParam * trackIn1 = track1->GetInnerParam(); + // + // + TTreeSRedirector * pcstream = GetDebugStreamer(); + Int_t ntracksSPD = vertexSPD->GetNContributors(); + Int_t ntracksTPC = vertexTPC->GetNContributors(); + // + AliESDfriendTrack *friendTrack0 = esdFriend->GetTrack(itrack0); + if (!friendTrack0) continue; + AliESDfriendTrack *friendTrack1 = esdFriend->GetTrack(itrack1); + if (!friendTrack1) continue; + TObject *calibObject; + AliTPCseed *seed0 = 0; + AliTPCseed *seed1 = 0; + // + for (Int_t l=0;(calibObject=friendTrack0->GetCalibObject(l));++l) { + if ((seed0=dynamic_cast(calibObject))) break; + } + for (Int_t l=0;(calibObject=friendTrack1->GetCalibObject(l));++l) { + if ((seed1=dynamic_cast(calibObject))) break; + } + // + if (pcstream){ + (*pcstream)<<"pairs"<< + "run="<SetDirectory(0); + } + if (fCosmicTree->GetEntries()==0){ + // + fCosmicTree->SetDirectory(0); + fCosmicTree->Branch("t0.",&track0); + fCosmicTree->Branch("t1.",&track1); + fCosmicTree->Branch("ft0.",&friendTrack0); + fCosmicTree->Branch("ft1.",&friendTrack1); + }else{ + fCosmicTree->SetBranchAddress("t0.",&track0); + fCosmicTree->SetBranchAddress("t1.",&track1); + fCosmicTree->SetBranchAddress("ft0.",&friendTrack0); + fCosmicTree->SetBranchAddress("ft1.",&friendTrack1); + } + fCosmicTree->Fill(); + } + } +} + + +void AliTPCcalibCosmic::Terminate(){ + // + // copy the cosmic tree to memory resident tree + // + static Int_t counter=0; + printf("AliTPCcalibCosmic::Terminate\t%d\n",counter); + counter++; + AliTPCcalibBase::Terminate(); +} + + +void AliTPCcalibCosmic::AddTree(TTree* treeOutput, TTree * treeInput){ + // + // Add the content of tree: + // Notice automatic copy of tree in ROOT does not work for such complicated tree + // + AliESDtrack *track0=new AliESDtrack; + AliESDtrack *track1=new AliESDtrack; + AliESDfriendTrack *ftrack0=new AliESDfriendTrack; + AliESDfriendTrack *ftrack1=new AliESDfriendTrack; + treeInput->SetBranchAddress("t0.",&track0); + treeInput->SetBranchAddress("t1.",&track1); + treeInput->SetBranchAddress("ft0.",&ftrack0); + treeInput->SetBranchAddress("ft1.",&ftrack1); + if (treeOutput->GetEntries()==0){ + // + treeOutput->SetDirectory(0); + treeOutput->Branch("t0.",&track0); + treeOutput->Branch("t1.",&track1); + treeOutput->Branch("ft0.",&ftrack0); + treeOutput->Branch("ft1.",&ftrack1); + }else{ + treeOutput->SetBranchAddress("t0.",&track0); + treeOutput->SetBranchAddress("t1.",&track1); + treeOutput->SetBranchAddress("ft0.",&ftrack0); + treeOutput->SetBranchAddress("ft1.",&ftrack1); + } + Int_t entries= treeInput->GetEntries(); + for (Int_t i=0; iGetEntry(i); + treeOutput->Fill(); + } +} + + + +void AliTPCcalibCosmic::MakeFitTree(TTree * treeInput, TTreeSRedirector *pcstream, const TObjArray * corrArray, Int_t step, Int_t run){ + // + // Make fit tree + // refit the tracks with original points + corrected points for each correction + // Input: + // treeInput - tree with cosmic tracks + // pcstream - debug output + + // Algorithm: + // Loop over pair of cosmic tracks: + // 1. Find trigger offset between cosmic event and "physic" trigger + // 2. Refit tracks with current transformation + // 3. Refit tracks using additional "primitive" distortion on top + // Best correction estimated as linear combination of corrections + // minimizing the observed distortions + // Observed distortions - matching in the y,z, snp, theta and 1/pt + // + const Double_t kResetCov=20.; + const Double_t kMaxDelta[5]={1,40,0.03,0.01,0.2}; + const Double_t kSigma=2.; + const Double_t kMaxTime=1050; + const Double_t kMaxSnp=0.8; + Int_t ncorr=corrArray->GetEntries(); + AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; + AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters(); + AliGRPObject* grp = AliTPCcalibDB::Instance()->GetGRP(run); + Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd()); + transform->SetCurrentRun(run); + transform->SetCurrentTimeStamp(TMath::Nint(time)); + Double_t covar[15]; + for (Int_t i=0;i<15;i++) covar[i]=0; + covar[0]=kSigma*kSigma; + covar[2]=kSigma*kSigma; + covar[5]=kSigma*kSigma/Float_t(150*150); + covar[9]=kSigma*kSigma/Float_t(150*150); + covar[14]=0.2*0.2; + Double_t *distortions = new Double_t[ncorr+1]; + + AliESDtrack *track0=new AliESDtrack; + AliESDtrack *track1=new AliESDtrack; + AliESDfriendTrack *ftrack0=new AliESDfriendTrack; + AliESDfriendTrack *ftrack1=new AliESDfriendTrack; + treeInput->SetBranchAddress("t0.",&track0); + treeInput->SetBranchAddress("t1.",&track1); + treeInput->SetBranchAddress("ft0.",&ftrack0); + treeInput->SetBranchAddress("ft1.",&ftrack1); + Int_t entries= treeInput->GetEntries(); + for (Int_t i=0; iGetEntry(i); + if (i%20==0) printf("%d\n",i); + if (!ftrack0->GetTPCOut()) continue; + if (!ftrack1->GetTPCOut()) continue; + AliTPCseed *seed0=0; + AliTPCseed *seed1=0; + TObject *calibObject; + for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) { + if ((seed0=dynamic_cast(calibObject))) break; + } + for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) { + if ((seed1=dynamic_cast(calibObject))) break; + } + if (!seed0) continue; + if (!seed1) continue; + if (TMath::Abs(seed0->GetSnp())>kMaxSnp) continue; + if (TMath::Abs(seed1->GetSnp())>kMaxSnp) continue; + // + // + Int_t nclA0=0, nclC0=0; // number of clusters + Int_t nclA1=0, nclC1=0; // number of clusters + Int_t ncl0=0,ncl1=0; + Double_t rmin0=300, rmax0=-300; // variables to estimate the time0 - trigger offset + Double_t rmin1=300, rmax1=-300; + Double_t tmin0=2000, tmax0=-2000; + Double_t tmin1=2000, tmax1=-2000; + // + // + // calculate trigger offset usig "missing clusters" + for (Int_t irow=0; irow<159; irow++){ + AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow); + if (cluster0 &&cluster0->GetX()>10){ + if (cluster0->GetX()GetX(); tmin0=cluster0->GetTimeBin();} + if (cluster0->GetX()>rmax0) { rmax0=cluster0->GetX(); tmax0=cluster0->GetTimeBin();} + ncl0++; + if (cluster0->GetDetector()%36<18) nclA0++; + if (cluster0->GetDetector()%36>=18) nclC0++; + } + AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow); + if (cluster1&&cluster1->GetX()>10){ + if (cluster1->GetX()GetX(); tmin1=cluster1->GetTimeBin();} + if (cluster1->GetX()>rmax1) { rmax1=cluster1->GetX(); tmax1=cluster1->GetTimeBin();} + ncl1++; + if (cluster1->GetDetector()%36<18) nclA1++; + if (cluster1->GetDetector()%36>=18) nclC1++; + } + } + Int_t cosmicType=0; // types of cosmic topology + if ((nclA0>nclC0) && (nclA1>nclC1)) cosmicType=0; // AA side + if ((nclA0nclC0) && (nclA1nclC1)) cosmicType=3; // CA side + //if ((nclA0>nclC0) && (nclA1nclC0) && (nclA11)&&TMath::Abs(track1->GetZ()-track0->GetZ())>4){ + cosmicType+=4; + deltaTime=0.5*(track1->GetZ()-track0->GetZ())/param->GetZWidth(); + if (nclA0>nclC0) deltaTime*=-1; // if A side track + } + // + TVectorD vectorDT(8); + Int_t crossCounter=0; + Double_t deltaTimeCross = AliTPCcalibCosmic::GetDeltaTime(rmin0, rmax0, rmin1, rmax1, tmin0, tmax0, tmin1, tmax1, TMath::Abs(track0->GetY()),vectorDT); + Bool_t isOKTrigger=kTRUE; + for (Int_t ic=0; ic<6;ic++) { + if (TMath::Abs(vectorDT[ic])>0) { + if (vectorDT[ic]+vectorDT[6]<0) isOKTrigger=kFALSE; + if (vectorDT[ic]+vectorDT[7]>kMaxTime) isOKTrigger=kFALSE; + if (isOKTrigger){ + crossCounter++; + } + } + } + Double_t deltaTimeCluster=deltaTime; + if ((cosmicType==0 || cosmicType==1) && crossCounter>0){ + deltaTimeCluster=deltaTimeCross; + cosmicType+=8; + } + if (nclA0*nclC0>0 || nclA1*nclC1>0) cosmicType+=16; // mixed A side C side - bad for visualization + // + // Apply current transformation + // + // + for (Int_t irow=0; irow<159; irow++){ + AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow); + if (cluster0 &&cluster0->GetX()>10){ + Double_t x0[3]={cluster0->GetRow(),cluster0->GetPad(),cluster0->GetTimeBin()+deltaTimeCluster}; + Int_t index0[1]={cluster0->GetDetector()}; + transform->Transform(x0,index0,0,1); + cluster0->SetX(x0[0]); + cluster0->SetY(x0[1]); + cluster0->SetZ(x0[2]); + // + } + AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow); + if (cluster1&&cluster1->GetX()>10){ + Double_t x1[3]={cluster1->GetRow(),cluster1->GetPad(),cluster1->GetTimeBin()+deltaTimeCluster}; + Int_t index1[1]={cluster1->GetDetector()}; + transform->Transform(x1,index1,0,1); + cluster1->SetX(x1[0]); + cluster1->SetY(x1[1]); + cluster1->SetZ(x1[2]); + } + } + // + // + Double_t alpha=track0->GetAlpha(); // rotation frame + Double_t cos = TMath::Cos(alpha); + Double_t sin = TMath::Sin(alpha); + Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass(); + AliExternalTrackParam btrack0=*(ftrack0->GetTPCOut()); + AliExternalTrackParam btrack1=*(ftrack1->GetTPCOut()); + btrack0.Rotate(alpha); + btrack1.Rotate(alpha); + // change the sign for track 1 + Double_t* par1=(Double_t*)btrack0.GetParameter(); + par1[3]*=-1; + par1[4]*=-1; + btrack0.AddCovariance(covar); + btrack1.AddCovariance(covar); + btrack0.ResetCovariance(kResetCov); + btrack1.ResetCovariance(kResetCov); + Bool_t isOK=kTRUE; + Bool_t isOKT=kTRUE; + TObjArray tracks0(ncorr+1); + TObjArray tracks1(ncorr+1); + // + Double_t dEdx0Tot=seed0->CookdEdxAnalytical(0.02,0.6,kTRUE); + Double_t dEdx1Tot=seed1->CookdEdxAnalytical(0.02,0.6,kTRUE); + Double_t dEdx0Max=seed0->CookdEdxAnalytical(0.02,0.6,kFALSE); + Double_t dEdx1Max=seed1->CookdEdxAnalytical(0.02,0.6,kFALSE); + //if (TMath::Abs((dEdx0Max+1)/(dEdx0Tot+1)-1.)>0.1) isOK=kFALSE; + //if (TMath::Abs((dEdx1Max+1)/(dEdx1Tot+1)-1.)>0.1) isOK=kFALSE; + ncl0=0; ncl1=0; + for (Int_t icorr=-1; icorr=0) corr = (AliTPCCorrection*)corrArray->At(icorr); + // + for (Int_t irow=159; irow>0; irow--){ + AliTPCclusterMI *cluster=seed0->GetClusterPointer(irow); + if (!cluster) continue; + if (!isOKT) break; + Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()}; + transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD); // transform to global + Float_t r[3]={rD[0],rD[1],rD[2]}; + if (corr){ + corr->DistortPoint(r, cluster->GetDetector()); + } + // + Double_t cov[3]={0.01,0.,0.01}; + Double_t lx =cos*r[0]+sin*r[1]; + Double_t ly =-sin*r[0]+cos*r[1]; + rD[1]=ly; rD[0]=lx; rD[2]=r[2]; //transform to track local + if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, lx,mass,1.,kFALSE)) isOKT=kFALSE;; + if (!rtrack0.Update(&rD[1],cov)) isOKT =kFALSE; + if (icorr<0) ncl0++; + } + // + for (Int_t irow=159; irow>0; irow--){ + AliTPCclusterMI *cluster=seed1->GetClusterPointer(irow); + if (!cluster) continue; + if (!isOKT) break; + Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()}; + transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD); + Float_t r[3]={rD[0],rD[1],rD[2]}; + if (corr){ + corr->DistortPoint(r, cluster->GetDetector()); + } + // + Double_t cov[3]={0.01,0.,0.01}; + Double_t lx =cos*r[0]+sin*r[1]; + Double_t ly =-sin*r[0]+cos*r[1]; + rD[1]=ly; rD[0]=lx; rD[2]=r[2]; + if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, lx,mass,1.,kFALSE)) isOKT=kFALSE; + if (!rtrack1.Update(&rD[1],cov)) isOKT=kFALSE; + if (icorr<0) ncl1++; + } + if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,10.,kFALSE)) isOKT=kFALSE; + if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,10.,kFALSE)) isOKT=kFALSE; + if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,1.,kFALSE)) isOKT=kFALSE; + if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,1.,kFALSE)) isOKT=kFALSE; + const Double_t *param0=rtrack0.GetParameter(); + const Double_t *param1=rtrack1.GetParameter(); + for (Int_t ipar=0; ipar<4; ipar++){ + if (TMath::Abs(param1[ipar]-param0[ipar])>kMaxDelta[ipar]) isOK=kFALSE; + } + tracks0.AddAt(rtrack0.Clone(), icorr+1); + tracks1.AddAt(rtrack1.Clone(), icorr+1); + AliExternalTrackParam out0=*(ftrack0->GetTPCOut()); + AliExternalTrackParam out1=*(ftrack1->GetTPCOut()); + Int_t nentries=TMath::Min(ncl0,ncl1); -strFit+="sign(Tr0.fP[1])++" -strFit+="sign(Tr0.fP[1])*(1-abs(Tr0.fP[1]/250))" - -TString * thetaParam = TStatToolkit::FitPlane(chain,"Tr0.fP[3]+Tr1.fP[3]", strFit.Data(),cutA, chi2,npoints,fitParam,covMatrix) + if (icorr<0) { + (*pcstream)<<"cosmic"<< + "isOK="<=0) corr = (AliTPCCorrection*)corrArray->At(icorr); + // + AliExternalTrackParam *param0=(AliExternalTrackParam *) tracks0.At(icorr+1); + AliExternalTrackParam *param1=(AliExternalTrackParam *) tracks1.At(icorr+1); + distortions[icorr+1]=param1->GetParameter()[ipar]-param0->GetParameter()[ipar]; + if (icorr>=0){ + distortions[icorr+1]-=distortions[0]; + } + // + if (icorr<0){ + Double_t bz=AliTrackerBase::GetBz(); + Double_t gxyz[3]; + param0->GetXYZ(gxyz); + Int_t dtype=20; + Double_t theta=param0->GetParameter()[3]; + Double_t phi = alpha; + Double_t snp = track0->GetInnerParam()->GetSnp(); + Double_t mean= distortions[0]; + Int_t index = param0->GetIndex(ipar,ipar); + Double_t rms=TMath::Sqrt(param1->GetCovariance()[index]+param1->GetCovariance()[index]); + if (crossCounter<1) rms*=1; + Double_t sector=9*phi/TMath::Pi(); + Double_t dsec = sector-TMath::Nint(sector+0.5); + Double_t gx=gxyz[0],gy=gxyz[1],gz=gxyz[2]; + Double_t refX=TMath::Sqrt(gx*gx+gy*gy); + Double_t dRrec=0; + // Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5; + Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5; + (*pcstream)<<"fit"<< // dump valus for fit + "run="<tmin => deltaT=Tcmax-tmax + // 1.b) rmin>Rcmin && tmintmax => deltaT=Tcmax-tmin + // // case the z matching gives proper time + // 1.c) rmaxTcmin && tmax deltaT=-tmax + // + // check algorithm: + // TCut cutStop = "(min(rmax0,rmax1)<235||abs(rmin0-rmin1)>10)"; // tracks not registered for full time + // Combinations: + // 0-1 - forbidden + // 0-2 - forbidden + // 0-3 - occur - wrong correlation + // 1-2 - occur - wrong correlation + // 1-3 - forbidden + // 2-3 - occur - small number of outlyers -20% + // Frequency: + // 0 - 106 + // 1 - 265 + // 2 - 206 + // 3 - 367 + // + const Double_t kMaxRCut=235; // max radius + const Double_t kMinRCut=TMath::Max(dcaR,90.); // min radius + const Double_t kMaxDCut=30; // max distance for minimal radius + const Double_t kMinTime=110; + const Double_t kMaxTime=950; + Double_t deltaT=0; + Int_t counter=0; + vectorDT[6]=TMath::Min(TMath::Min(tmin0,tmin1),TMath::Min(tmax0,tmax1)); + vectorDT[7]=TMath::Max(TMath::Max(tmin0,tmin1),TMath::Max(tmax0,tmax1)); + if (TMath::Min(rmax0,rmax1)0 + if (rmax0tmin0) vectorDT[0]=kMaxTime-tmax0; // disapear at CE + if (rmax1tmin1) vectorDT[1]=kMaxTime-tmax1; // disapear at CE + // min cross - deltaT<0 - OK they are correlated + if (rmax0kMinTime && tmax0kMinTime && tmax1 kMinRCut+kMaxDCut && tmin0 tmax0) vectorDT[4]=kMaxTime-tmin0; + if (rmin1> kMinRCut+kMaxDCut && tmin1 tmax1) vectorDT[5]=kMaxTime-tmin1; + Bool_t isOK=kTRUE; + for (Int_t i=0; i<6;i++) { + if (TMath::Abs(vectorDT[i])>0) { + counter++; + if (vectorDT[i]+vectorDT[6]<0) isOK=kFALSE; + if (vectorDT[i]+vectorDT[7]>kMaxTime) isOK=kFALSE; + if (isOK) deltaT=vectorDT[i]; + } + } + return deltaT; +}