X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TPC%2FAliTPCcalibTime.cxx;h=e25ca3916635a744f4a1fe9b410a52ad4e49903c;hp=e06510db53da3a7c6cc7ccdb970ced90345945fd;hb=d8171ef5bd201fca88ec174418e26c88926df1d6;hpb=b9908d0b8e9e96b9f6cf84fd98a9907fd2a4e5c6 diff --git a/TPC/AliTPCcalibTime.cxx b/TPC/AliTPCcalibTime.cxx index e06510db53d..e25ca391663 100644 --- a/TPC/AliTPCcalibTime.cxx +++ b/TPC/AliTPCcalibTime.cxx @@ -96,12 +96,15 @@ Comments to be written here: #include "AliTPCseed.h" #include "AliTrackPointArray.h" #include "AliTracker.h" +#include "AliKFVertex.h" +#include ClassImp(AliTPCcalibTime) AliTPCcalibTime::AliTPCcalibTime() - :AliTPCcalibBase(), + :AliTPCcalibBase(), + fMemoryMode(1), // 0 -do not fill THnSparse with residuals 1- fill only important QA THn 2 - Fill all THnsparse for calibration fLaser(0), // pointer to laser calibration fDz(0), // current delta z fCutMaxD(3), // maximal distance in rfi ditection @@ -109,10 +112,13 @@ AliTPCcalibTime::AliTPCcalibTime() fCutTheta(0.03), // maximal distan theta fCutMinDir(-0.99), // direction vector products fCutTracks(100), + fArrayLaserA(0), //laser fit parameters C + fArrayLaserC(0), //laser fit parameters A fArrayDz(0), //NEW! Tmap of V drifts for different triggers fAlignITSTPC(0), //alignemnt array ITS TPC match fAlignTRDTPC(0), //alignemnt array TRD TPC match fAlignTOFTPC(0), //alignemnt array TOF TPC match + fTimeKalmanBin(60*15), //time bin width for kalman - 15 minutes default fTimeBins(0), fTimeStart(0), fTimeEnd(0), @@ -144,12 +150,27 @@ AliTPCcalibTime::AliTPCcalibTime() fResHistoTPCTRD[i]=0; fResHistoTPCTOF[i]=0; fResHistoTPCvertex[i]=0; + fTPCVertex[i]=0; + } + for (Int_t i=0;i<12;i++) { + fTPCVertex[i]=0; + } + for (Int_t i=0;i<5;i++) { + fTPCVertexCorrelation[i]=0; + } + static Int_t counter=0; + if (1) { + TTimeStamp s; + Int_t time=s; + AliInfo(Form("Counter Constructor\t%d\t%d",counter,time)); + counter++; } } -AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeVdrift) +AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeVdrift, Int_t memoryMode) :AliTPCcalibBase(), + fMemoryMode(memoryMode), // 0 -do not fill THnSparse with residuals 1- fill only important QA THn 2 - Fill all THnsparse for calibration fLaser(0), // pointer to laser calibration fDz(0), // current delta z fCutMaxD(5*0.5356), // maximal distance in rfi ditection @@ -157,10 +178,13 @@ AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t fCutTheta(5*0.004644),// maximal distan theta fCutMinDir(-0.99), // direction vector products fCutTracks(100), + fArrayLaserA(new TObjArray(1000)), //laser fit parameters C + fArrayLaserC(new TObjArray(1000)), //laser fit parameters A fArrayDz(0), //Tmap of V drifts for different triggers fAlignITSTPC(0), //alignemnt array ITS TPC match fAlignTRDTPC(0), //alignemnt array TRD TPC match fAlignTOFTPC(0), //alignemnt array TOF TPC match + fTimeKalmanBin(60*15), //time bin width for kalman - 15 minutes default fTimeBins(0), fTimeStart(0), fTimeEnd(0), @@ -266,13 +290,27 @@ AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t fCosmiMatchingHisto[7]=new TH1F("Cosmics matching","p2-isPair",100,-10*0.01134 ,10*0.01134 ); fCosmiMatchingHisto[8]=new TH1F("Cosmics matching","p3-isPair",100,-10*0.004644,10*0.004644); fCosmiMatchingHisto[9]=new TH1F("Cosmics matching","p4-isPair",100,-10*0.03773 ,10*0.03773 ); + for (Int_t i=0;i<12;i++) { + fTPCVertex[i]=0; + } + for (Int_t i=0;i<5;i++) { + fTPCVertexCorrelation[i]=0; + } BookDistortionMaps(); + } AliTPCcalibTime::~AliTPCcalibTime(){ // // Virtual Destructor // + static Int_t counter=0; + if (1) { + TTimeStamp s; + Int_t time=s; + AliInfo(Form("Counter Destructor\t%s\t%d\t%d",GetName(),counter,time)); + counter++; + } for(Int_t i=0;i<3;i++){ if(fHistVdriftLaserA[i]){ delete fHistVdriftLaserA[i]; @@ -309,7 +347,13 @@ AliTPCcalibTime::~AliTPCcalibTime(){ fResHistoTPCvertex[i]=0; } - + if (fTPCVertex) { + for (Int_t i=0;i<12;i++) delete fTPCVertex[i]; + } + if (fTPCVertexCorrelation) { + for (Int_t i=0;i<5;i++) delete fTPCVertexCorrelation[i]; + } + fAlignITSTPC->SetOwner(kTRUE); fAlignTRDTPC->SetOwner(kTRUE); fAlignTOFTPC->SetOwner(kTRUE); @@ -444,6 +488,20 @@ void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){ npointsC= TMath::Nint(vdriftC[3]); chi2C= vdriftC[4]; + if (npointsA>kMinTracksSide || npointsC>kMinTracksSide){ + TVectorD *fitA = new TVectorD(6); + TVectorD *fitC = new TVectorD(6); + for (Int_t ipar=0; ipar<5; ipar++){ + (*fitA)[ipar]=vdriftA[ipar]; + (*fitC)[ipar]=vdriftC[ipar]; + } + (*fitA)[5]=fTime; + (*fitC)[5]=fTime; + fArrayLaserA->AddLast(fitA); + fArrayLaserC->AddLast(fitC); + } + // + TTimeStamp tstamp(fTime); Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0); Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1); @@ -473,33 +531,6 @@ void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){ fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA); fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC); } - -// THnSparse* curHist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift); -// TString shortName=curHist->ClassName(); -// shortName+="_MEAN_DRIFT_LASER_"; -// delete curHist; -// curHist=NULL; -// TString name=""; - -// name=shortName; -// name+=event->GetFiredTriggerClasses(); -// name.ToUpper(); -// curHist=(THnSparseF*)fArrayDz->FindObject(name); -// if(!curHist){ -// curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift); -// fArrayDz->AddLast(curHist); -// } -// curHist->Fill(vecDrift); - -// name=shortName; -// name+="ALL"; -// name.ToUpper(); -// curHist=(THnSparseF*)fArrayDz->FindObject(name); -// if(!curHist){ -// curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift); -// fArrayDz->AddLast(curHist); -// } -// curHist->Fill(vecDrift); } void AliTPCcalibTime::ProcessCosmic(const AliESDEvent *const event){ @@ -745,10 +776,232 @@ void AliTPCcalibTime::ProcessCosmic(const AliESDEvent *const event){ if (GetDebugLevel()>20) printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data()); } -void AliTPCcalibTime::ProcessBeam(const AliESDEvent *const /*event*/){ +void AliTPCcalibTime::ProcessBeam(const AliESDEvent *const event){ // - // Not special treatment yet - the same for cosmic and physic event // + // + const Int_t kMinClusters =80; + const Int_t kMinTracks =2; // minimal number of tracks to define the vertex + const Int_t kMinTracksVertex=30; // minimal number of tracks to define the cumulative vertex + const Double_t kMaxTgl =1.2; // maximal Tgl (z angle) + const Double_t kMinPt =0.2; // minimal pt + const Double_t kMaxD0 =5.; // cut on distance to the primary vertex first guess + const Double_t kMaxZ0 =20; + const Double_t kMaxD =2.5; // cut on distance to the primary vertex + const Double_t kMaxZ =4; // maximal z distance between tracks form the same side + const Double_t kMaxChi2 =15; // maximal chi2 of the TPCvertex + const Double_t kCumulCovarXY=0.003; //increase the error of cumul vertex 30 microns profile + const Double_t kCumulCovarZ=250.; //increase the error of cumul vertex + const Double_t kMaxDvertex = 1.0; // cut to accept the vertex; + // + Int_t flags=0; + const Int_t kBuffSize=100; + static Double_t deltaZ[kBuffSize]={0}; + static Int_t counterZ=0; + static AliKFVertex cumulVertexA, cumulVertexC, cumulVertexAC; // cumulative vertex + AliKFVertex vertexA, vertexC; + + Float_t dca0[2]={0,0}; + Double_t dcaVertex[2]={0,0}; + Int_t ntracks=event->GetNumberOfTracks(); + if (ntracks==0) return; + // + AliESDfriend *esdFriend=(AliESDfriend*)(((AliESDEvent*)event)->FindListObject("AliESDfriend")); + // + // Divide tracks to A and C side tracks - using the cluster indexes + TObjArray tracksA(ntracks); + TObjArray tracksC(ntracks); + // + AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD(); + AliESDVertex *vertex = (AliESDVertex *)event->GetPrimaryVertex(); + AliESDVertex *vertexTracks = (AliESDVertex *)event->GetPrimaryVertexTracks(); + Double_t vertexZA[10000], vertexZC[10000]; + // + Int_t ntracksA= 0; + Int_t ntracksC= 0; + // + for (Int_t itrack=0;itrackGetTrack(itrack); + AliESDfriendTrack *friendTrack = esdFriend->GetTrack(itrack); + if (!friendTrack) continue; + if (TMath::Abs(track->GetTgl())>kMaxTgl) continue; + if (TMath::Abs(track->Pt())GetInnerParam(); + TObject *calibObject=0; + AliTPCseed *seed = 0; + Int_t nA=0, nC=0; + for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) if ((seed=dynamic_cast(calibObject))) break; + if (seed) { + for (Int_t irow=159;irow>0;irow--) { + AliTPCclusterMI *cl=seed->GetClusterPointer(irow); + if (!cl) continue; + if ((cl->GetDetector()%36)<18) nA++; + if ((cl->GetDetector()%36)>=18) nC++; + } + if ((nA>kMinClusters || nC>kMinClusters) && (nA*nC==0) ){ + track->GetImpactParameters(dca0[0],dca0[1]); + if (TMath::Abs(dca0[0])>kMaxD0) continue; + if (TMath::Abs(dca0[1])>kMaxZ0) continue; + AliExternalTrackParam pTPCvertex(*trackIn); + if (!AliTracker::PropagateTrackToBxByBz(&pTPCvertex,4.+4.*TMath::Abs(dca0[0]),0.1,2,kTRUE)) continue; + pTPCvertex.PropagateToDCA(vertex,AliTracker::GetBz(), kMaxD, dcaVertex,0); + if (TMath::Abs(dcaVertex[0])>kMaxD) continue; + if (nA>kMinClusters &&nC==0) { tracksA.AddLast(pTPCvertex.Clone()); vertexZA[ntracksA++] = pTPCvertex.GetZ();} + if (nC>kMinClusters &&nA==0) {tracksC.AddLast(pTPCvertex.Clone()); vertexZC[ntracksC++] = pTPCvertex.GetZ();} + } + } + } + Double_t medianZA=TMath::Median(ntracksA, vertexZA); // tracks median + Double_t medianZC=TMath::Median(ntracksC, vertexZC); // tracks median + // + ntracksA= tracksA.GetEntriesFast(); + ntracksC= tracksC.GetEntriesFast(); + if (ntracksA>kMinTracks && ntracksC>kMinTracks){ + deltaZ[counterZ%kBuffSize]=medianZA-medianZC; + counterZ+=1; + Double_t medianDelta=(counterZ>=kBuffSize)? TMath::Median(kBuffSize, deltaZ): TMath::Median(counterZ, deltaZ); + if (TMath::Abs(medianDelta-(medianZA-medianZC))>kMaxZ) flags+=16; + // increse the error of cumulative vertex at the beginning of event + cumulVertexA.Covariance(0,0)+=kCumulCovarXY*kCumulCovarXY; + cumulVertexA.Covariance(1,1)+=kCumulCovarXY*kCumulCovarXY; + cumulVertexA.Covariance(2,2)+=kCumulCovarZ*kCumulCovarZ; + cumulVertexC.Covariance(0,0)+=kCumulCovarXY*kCumulCovarXY; + cumulVertexC.Covariance(1,1)+=kCumulCovarXY*kCumulCovarXY; + cumulVertexC.Covariance(2,2)+=kCumulCovarZ*kCumulCovarZ; + cumulVertexAC.Covariance(0,0)+=kCumulCovarXY*kCumulCovarXY; + cumulVertexAC.Covariance(1,1)+=kCumulCovarXY*kCumulCovarXY; + cumulVertexAC.Covariance(2,2)+=kCumulCovarZ*kCumulCovarZ; + // + for (Int_t iA=0; iAGetZ()-medianZA)>kMaxZ) continue; + AliKFParticle part(*aliTrack,211); + vertexA+=part; + } + for (Int_t iC=0; iCGetZ()-medianZC)>kMaxZ) continue; + AliKFParticle part(*aliTrack,211); + vertexC+=part; + } + // + if (vertexA.GetNDF()kMaxZ) flags+=1; //apply cuts + if (TMath::Abs(vertexC.Z()-medianZC)>kMaxZ) flags+=2; + if (TMath::Abs(vertexA.GetChi2()/vertexA.GetNDF()+vertexC.GetChi2()/vertexC.GetNDF())> kMaxChi2) flags+=4; + // + if (flags==0){ + for (Int_t iA=0; iAGetZ()-medianZA)>kMaxZ) continue; + AliKFParticle part(*aliTrack,211); + cumulVertexA+=part; + cumulVertexAC+=part; + } + for (Int_t iC=0; iCGetZ()-medianZC)>kMaxZ) continue; + AliKFParticle part(*aliTrack,211); + cumulVertexC+=part; + cumulVertexAC+=part; + } + // + if (TMath::Abs(cumulVertexA.X()-vertexA.X())>kMaxDvertex) flags+=64; + if (TMath::Abs(cumulVertexA.Y()-vertexA.Y())>kMaxDvertex) flags+=64; + if (TMath::Abs(cumulVertexA.Z()-vertexA.Z())>kMaxDvertex) flags+=64; + // + if (TMath::Abs(cumulVertexC.X()-vertexC.X())>kMaxDvertex) flags+=64; + if (TMath::Abs(cumulVertexC.Y()-vertexC.Y())>kMaxDvertex) flags+=64; + if (TMath::Abs(cumulVertexC.Z()-vertexC.Z())>kMaxDvertex) flags+=64; + + + if ( flags==0 && cumulVertexC.GetNDF()>kMinTracksVertex&&cumulVertexA.GetNDF()>kMinTracksVertex){ + Double_t cont[2]={0,fTime}; + // + cont[0]= cumulVertexA.X(); + fTPCVertex[0]->Fill(cont); + cont[0]= cumulVertexC.X(); + fTPCVertex[1]->Fill(cont); + cont[0]= 0.5*(cumulVertexA.X()-cumulVertexC.X()); + fTPCVertex[2]->Fill(cont); + cont[0]= 0.5*(cumulVertexA.X()+cumulVertexC.X())-vertexSPD->GetX(); + fTPCVertex[3]->Fill(cont); + // + cont[0]= cumulVertexA.Y(); + fTPCVertex[4]->Fill(cont); + cont[0]= cumulVertexC.Y(); + fTPCVertex[5]->Fill(cont); + cont[0]= 0.5*(cumulVertexA.Y()-cumulVertexC.Y()); + fTPCVertex[6]->Fill(cont); + cont[0]= 0.5*(cumulVertexA.Y()+cumulVertexC.Y())-vertexSPD->GetY(); + fTPCVertex[7]->Fill(cont); + // + // + cont[0]= 0.5*(cumulVertexA.Z()+cumulVertexC.Z()); + fTPCVertex[8]->Fill(cont); + cont[0]= 0.5*(cumulVertexA.Z()-cumulVertexC.Z()); + fTPCVertex[9]->Fill(cont); + cont[0]= 0.5*(cumulVertexA.Z()-cumulVertexC.Z()); + fTPCVertex[10]->Fill(cont); + cont[0]= 0.5*(cumulVertexA.Z()+cumulVertexC.Z())-vertexSPD->GetZ(); + fTPCVertex[11]->Fill(cont); + // + Double_t correl[2]={0,0}; + // + correl[0]=cumulVertexC.Z(); + correl[1]=cumulVertexA.Z(); + fTPCVertexCorrelation[0]->Fill(correl); // fill A side :TPC + correl[0]=cumulVertexA.Z(); + correl[1]=cumulVertexC.Z(); + fTPCVertexCorrelation[1]->Fill(correl); // fill C side :TPC + // + correl[0]=vertexSPD->GetZ(); + correl[1]=cumulVertexA.Z()-correl[0]; + fTPCVertexCorrelation[2]->Fill(correl); // fill A side :ITS + correl[1]=cumulVertexC.Z()-correl[0]; + fTPCVertexCorrelation[3]->Fill(correl); // fill C side :ITS + correl[1]=0.5*(cumulVertexA.Z()+cumulVertexC.Z())-correl[0]; + fTPCVertexCorrelation[4]->Fill(correl); // fill C side :ITS + } + } + TTreeSRedirector *cstream = GetDebugStreamer(); + if (cstream){ + /* + TCut cutChi2= "sqrt(vA.fChi2/vA.fNDF+vC.fChi2/vC.fNDF)<10"; // chi2 Cut e.g 10 + TCut cutXY= "sqrt((vA.fP[0]-vC.fP[0])^2+(vA.fP[0]-vC.fP[1])^2)<5"; // vertex Cut + TCut cutZ= "abs(vA.fP[2]-mZA)<3&&abs(vC.fP[2]-mZC)<5"; // vertex Cut + tree->Draw("sqrt(vA.fChi2/vA.fNDF)","sqrt(vA.fChi2/vA.fNDF)<100","") + + */ + //vertexA.Print(); + //vertexC.Print(); + (*cstream)<<"vertexTPC"<< + "flags="<MakeIterator(); AliTPCcalibTime* cal = 0; - + // while ((cal = (AliTPCcalibTime*)iter->Next())) { if (!cal->InheritsFrom(AliTPCcalibTime::Class())) { Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName()); @@ -863,35 +1116,61 @@ Long64_t AliTPCcalibTime::Merge(TCollection *const li) { } } // - for (Int_t imeas=0; imeas<5; imeas++){ - if (cal->GetResHistoTPCCE(imeas) && cal->GetResHistoTPCCE(imeas)){ - fResHistoTPCCE[imeas]->Add(cal->fResHistoTPCCE[imeas]); - }else{ - fResHistoTPCCE[imeas]=(THnSparse*)cal->fResHistoTPCCE[imeas]->Clone(); + if (fTPCVertexCorrelation && cal->fTPCVertexCorrelation){ + for (Int_t imeas=0; imeas<5; imeas++){ + if (fTPCVertexCorrelation[imeas] && cal->fTPCVertexCorrelation[imeas]) fTPCVertexCorrelation[imeas]->Add(cal->fTPCVertexCorrelation[imeas]); + } + } + + if (fTPCVertex && cal->fTPCVertex) + for (Int_t imeas=0; imeas<12; imeas++){ + if (fTPCVertex[imeas] && cal->fTPCVertex[imeas]) fTPCVertex[imeas]->Add(cal->fTPCVertex[imeas]); } - if (cal->GetResHistoTPCITS(imeas) && cal->GetResHistoTPCITS(imeas)){ - fResHistoTPCITS[imeas]->Add(cal->fResHistoTPCITS[imeas]); - fResHistoTPCvertex[imeas]->Add(cal->fResHistoTPCvertex[imeas]); + + if (fMemoryMode>0) for (Int_t imeas=0; imeas<5; imeas++){ + if (fMemoryMode>1){ + if ( cal->GetResHistoTPCCE(imeas) && cal->GetResHistoTPCCE(imeas)){ + fResHistoTPCCE[imeas]->Add(cal->fResHistoTPCCE[imeas]); + }else{ + fResHistoTPCCE[imeas]=(THnSparse*)cal->fResHistoTPCCE[imeas]->Clone(); + } + } + // + if ((fMemoryMode>0) &&cal->GetResHistoTPCITS(imeas) && cal->GetResHistoTPCITS(imeas)){ + if (fMemoryMode>1 || (imeas%2)==1) fResHistoTPCITS[imeas]->Add(cal->fResHistoTPCITS[imeas]); + if (fMemoryMode>1) fResHistoTPCvertex[imeas]->Add(cal->fResHistoTPCvertex[imeas]); } - if (cal->fResHistoTPCTRD[imeas]){ + // + if ((fMemoryMode>1) && cal->fResHistoTPCTRD[imeas]){ if (fResHistoTPCTRD[imeas]) fResHistoTPCTRD[imeas]->Add(cal->fResHistoTPCTRD[imeas]); else fResHistoTPCTRD[imeas]=(THnSparse*)cal->fResHistoTPCTRD[imeas]->Clone(); } - if (cal->fResHistoTPCTOF[imeas]){ + // + if ((fMemoryMode>1) && cal->fResHistoTPCTOF[imeas]){ if (fResHistoTPCTOF[imeas]) fResHistoTPCTOF[imeas]->Add(cal->fResHistoTPCTOF[imeas]); else fResHistoTPCTOF[imeas]=(THnSparse*)cal->fResHistoTPCTOF[imeas]->Clone(); } + // + if (cal->fArrayLaserA){ + fArrayLaserA->Expand(fArrayLaserA->GetEntriesFast()+cal->fArrayLaserA->GetEntriesFast()); + fArrayLaserC->Expand(fArrayLaserC->GetEntriesFast()+cal->fArrayLaserC->GetEntriesFast()); + for (Int_t ical=0; icalfArrayLaserA->GetEntriesFast(); ical++){ + if (cal->fArrayLaserA->UncheckedAt(ical)) fArrayLaserA->AddLast(cal->fArrayLaserA->UncheckedAt(ical)->Clone()); + if (cal->fArrayLaserC->UncheckedAt(ical)) fArrayLaserC->AddLast(cal->fArrayLaserC->UncheckedAt(ical)->Clone()); + } + } + } TObjArray* addArray=cal->GetHistoDrift(); if(!addArray) return 0; TIterator* iterator = addArray->MakeIterator(); iterator->Reset(); THnSparse* addHist=NULL; - while((addHist=(THnSparseF*)iterator->Next())){ + if ((fMemoryMode>1)) while((addHist=(THnSparseF*)iterator->Next())){ if(!addHist) continue; addHist->Print(); THnSparse* localHist=(THnSparseF*)fArrayDz->FindObject(addHist->GetName()); @@ -1073,8 +1352,8 @@ void AliTPCcalibTime::ProcessSame(AliESDtrack *const track, AliESDfriendTrack * Double_t xyz[3]={0,0.,0.0}; Double_t bz =0; Int_t nclIn=0,nclOut=0; - trackIn.ResetCovariance(10.); - trackOut.ResetCovariance(10.); + trackIn.ResetCovariance(1000.); + trackOut.ResetCovariance(1000.); // //2.a Refit inner // @@ -1209,41 +1488,40 @@ void AliTPCcalibTime::ProcessAlignITS(AliESDtrack *const track, AliESDfriendTra const Int_t kMinITS = 3; // minimal number of ITS cluster const Double_t kMinZ = 10; // maximal dz distance const Double_t kMaxDy = 2.; // maximal dy distance - const Double_t kMaxAngle= 0.015; // maximal angular distance + const Double_t kMaxAngle= 0.07; // maximal angular distance const Double_t kSigmaCut= 5; // maximal sigma distance to median const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction + const Double_t kT0Err = 3.; // initial uncertainty of the T0 time const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman const Double_t kMinPt = 0.3; // minimal pt + const Double_t kMax1Pt=0.5; //maximal 1/pt distance const Int_t kN=50; // deepnes of history static Int_t kglast=0; static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]}; - /* - 0. Standrd cuts: - TCut cut="abs(pTPC.fP[2]-pITS.fP[2])<0.01&&abs(pTPC.fP[3]-pITS.fP[3])<0.01&&abs(pTPC.fP[2]-pITS.fP[2])<1"; - */ // // 0. Apply standard cuts - // + // Int_t dummycl[1000]; if (track->GetTPCNcls()GetITSclusters(dummycl)IsOn(AliESDtrack::kTPCrefit)) return; - if (!friendTrack->GetITSOut()) return; if (!track->GetInnerParam()) return; if (!track->GetOuterParam()) return; if (track->GetInnerParam()->Pt()GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return; - if (TMath::Abs(track->GetInnerParam()->GetZ())GetInnerParam()->GetZ())GetInnerParam()->GetX()>90) return; // AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetInnerParam())); - AliExternalTrackParam pITS(*(friendTrack->GetITSOut())); // ITS standalone if possible - AliExternalTrackParam pITS2(*(friendTrack->GetITSOut())); //TPC-ITS track - pITS2.Rotate(pTPC.GetAlpha()); - // pITS2.PropagateTo(pTPC.GetX(),fMagF); - AliTracker::PropagateTrackToBxByBz(&pITS2,pTPC.GetX(),0.1,0.1,kFALSE); + // + AliExternalTrackParam pITS; // ITS standalone if possible + AliExternalTrackParam pITS2; //TPC-ITS track + if (friendTrack->GetITSOut()){ + pITS2=(*(friendTrack->GetITSOut())); //TPC-ITS track - snapshot ITS out + pITS2.Rotate(pTPC.GetAlpha()); + AliTracker::PropagateTrackToBxByBz(&pITS2,pTPC.GetX(),0.1,0.1,kFALSE); + } AliESDfriendTrack *itsfriendTrack=0; // @@ -1252,21 +1530,27 @@ void AliTPCcalibTime::ProcessAlignITS(AliESDtrack *const track, AliESDfriendTra Bool_t hasAlone=kFALSE; Int_t ntracks=event->GetNumberOfTracks(); for (Int_t i=0; iGetTrack(i); - if (trackS->GetTPCNcls()>0) continue; //continue if has TPC info + AliESDtrack * trackITS = event->GetTrack(i); + if (!trackITS) continue; + if (trackITS->GetITSclusters(dummycl)GetTrack(i); if (!itsfriendTrack) continue; if (!itsfriendTrack->GetITSOut()) continue; - if (TMath::Abs(pITS2.GetTgl()-itsfriendTrack->GetITSOut()->GetTgl())> kMaxAngle) continue; + + if (TMath::Abs(pTPC.GetTgl()-itsfriendTrack->GetITSOut()->GetTgl())> kMaxAngle) continue; + if (TMath::Abs(pTPC.GetSigned1Pt()-itsfriendTrack->GetITSOut()->GetSigned1Pt())> kMax1Pt) continue; pITS=(*(itsfriendTrack->GetITSOut())); // pITS.Rotate(pTPC.GetAlpha()); - //pITS.PropagateTo(pTPC.GetX(),fMagF); AliTracker::PropagateTrackToBxByBz(&pITS,pTPC.GetX(),0.1,0.1,kFALSE); - if (TMath::Abs(pITS2.GetY()-pITS.GetY())> kMaxDy) continue; + if (TMath::Abs(pTPC.GetY()-pITS.GetY())> kMaxDy) continue; + if (TMath::Abs(pTPC.GetSnp()-pITS.GetSnp())> kMaxAngle) continue; hasAlone=kTRUE; } - if (!hasAlone) pITS=pITS2; + if (!hasAlone) { + if (track->GetITSclusters(dummycl)kMaxDy) return; if (TMath::Abs(pITS.GetSnp()-pTPC.GetSnp())>kMaxAngle) return; @@ -1301,7 +1585,7 @@ void AliTPCcalibTime::ProcessAlignITS(AliESDtrack *const track, AliESDfriendTra // // 3. Update alignment // - Int_t htime = fTime/3600; //time in hours + Int_t htime = (fTime-fTimeKalmanBin/2)/fTimeKalmanBin; //time bins number if (fAlignITSTPC->GetEntriesFast()Expand(htime*2+20); } @@ -1311,6 +1595,7 @@ void AliTPCcalibTime::ProcessAlignITS(AliESDtrack *const track, AliESDfriendTra align=new AliRelAlignerKalman(); align->SetRunNumber(fRun); (*align->GetStateCov())(6,6)=kVdErr*kVdErr; + (*align->GetStateCov())(7,7)=kT0Err*kT0Err; (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr; align->SetOutRejSigma(kOutCut+kOutCut*kN); align->SetRejectOutliers(kFALSE); @@ -1320,7 +1605,12 @@ void AliTPCcalibTime::ProcessAlignITS(AliESDtrack *const track, AliESDfriendTra fAlignITSTPC->AddAt(align,htime); } align->AddTrackParams(&pITS,&pTPC); - align->SetTimeStamp(fTime); + Double_t averageTime = fTime; + if (align->GetTimeStamp()>0&&align->GetNUpdates()>0){ + averageTime=((Double_t(align->GetTimeStamp())*Double_t(align->GetNUpdates())+Double_t(fTime)))/(Double_t(align->GetNUpdates())+1.); + } + align->SetTimeStamp(Int_t(averageTime)); + align->SetRunNumber(fRun ); Float_t dca[2],cov[3]; track->GetImpactParameters(dca,cov); @@ -1386,6 +1676,7 @@ void AliTPCcalibTime::ProcessAlignTRD(AliESDtrack *const track, AliESDfriendTra const Double_t kMaxAngle= 0.1; // maximal angular distance const Double_t kSigmaCut= 10; // maximal sigma distance to median const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction + const Double_t kT0Err = 3.; // initial uncertainty of the T0 time const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman const Double_t kRefX = 275; // reference X @@ -1451,7 +1742,8 @@ void AliTPCcalibTime::ProcessAlignTRD(AliESDtrack *const track, AliESDfriendTra // // 3. Update alignment // - Int_t htime = fTime/3600; //time in hours + //Int_t htime = fTime/3600; //time in hours + Int_t htime = (Int_t)(fTime-fTimeKalmanBin/2)/fTimeKalmanBin; //time in half hour if (fAlignTRDTPC->GetEntriesFast()Expand(htime*2+20); } @@ -1461,6 +1753,7 @@ void AliTPCcalibTime::ProcessAlignTRD(AliESDtrack *const track, AliESDfriendTra align=new AliRelAlignerKalman(); align->SetRunNumber(fRun); (*align->GetStateCov())(6,6)=kVdErr*kVdErr; + (*align->GetStateCov())(7,7)=kT0Err*kT0Err; (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr; align->SetOutRejSigma(kOutCut+kOutCut*kN); align->SetRejectOutliers(kFALSE); @@ -1469,7 +1762,16 @@ void AliTPCcalibTime::ProcessAlignTRD(AliESDtrack *const track, AliESDfriendTra fAlignTRDTPC->AddAt(align,htime); } align->AddTrackParams(&pTRD,&pTPC); - align->SetTimeStamp(fTime); + //align->SetTimeStamp(fTime); + Double_t averageTime = fTime; + if (align->GetTimeStamp()>0 && align->GetNUpdates()>0) { + averageTime = (((Double_t)fTime) + ((Double_t)align->GetTimeStamp())*align->GetNUpdates()) / (align->GetNUpdates() + 1.); + //printf("align->GetTimeStamp() %d, align->GetNUpdates() %d \n", align->GetTimeStamp(), align->GetNUpdates()); + } + align->SetTimeStamp((Int_t)averageTime); + + //printf("fTime %d, averageTime %d \n", fTime, (Int_t)averageTime); + align->SetRunNumber(fRun ); Float_t dca[2],cov[3]; track->GetImpactParameters(dca,cov); @@ -1527,9 +1829,10 @@ void AliTPCcalibTime::ProcessAlignTOF(AliESDtrack *const track, AliESDfriendTra const Int_t kMinTPC = 80; // minimal number of TPC cluster // const Double_t kMinZ = 10; // maximal dz distance const Double_t kMaxDy = 5.; // maximal dy distance - const Double_t kMaxAngle= 0.015; // maximal angular distance + const Double_t kMaxAngle= 0.05; // maximal angular distance const Double_t kSigmaCut= 5; // maximal sigma distance to median const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction + const Double_t kT0Err = 3.; // initial uncertainty of the T0 time const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman @@ -1616,7 +1919,8 @@ void AliTPCcalibTime::ProcessAlignTOF(AliESDtrack *const track, AliESDfriendTra // // 3. Update alignment // - Int_t htime = fTime/3600; //time in hours + //Int_t htime = fTime/3600; //time in hours + Int_t htime = (Int_t)(fTime-fTimeKalmanBin)/fTimeKalmanBin; //time bin if (fAlignTOFTPC->GetEntriesFast()Expand(htime*2+20); } @@ -1626,6 +1930,7 @@ void AliTPCcalibTime::ProcessAlignTOF(AliESDtrack *const track, AliESDfriendTra align=new AliRelAlignerKalman(); align->SetRunNumber(fRun); (*align->GetStateCov())(6,6)=kVdErr*kVdErr; + (*align->GetStateCov())(7,7)=kT0Err*kT0Err; (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr; align->SetOutRejSigma(kOutCut+kOutCut*kN); align->SetRejectOutliers(kFALSE); @@ -1639,7 +1944,16 @@ void AliTPCcalibTime::ProcessAlignTOF(AliESDtrack *const track, AliESDfriendTra if (TMath::Abs(dca[0])SetTimeStamp(fTime); + //align->SetTimeStamp(fTime); + Double_t averageTime = fTime; + if (align->GetTimeStamp()>0 && align->GetNUpdates()>0) { + averageTime = (((Double_t)fTime) + ((Double_t)align->GetTimeStamp())*align->GetNUpdates()) / (align->GetNUpdates() + 1.); + //printf("align->GetTimeStamp() %d, align->GetNUpdates() %d \n", align->GetTimeStamp(), align->GetNUpdates()); + } + align->SetTimeStamp((Int_t)averageTime); + + //printf("fTime %d, averageTime %d \n", fTime, (Int_t)averageTime); + align->SetRunNumber(fRun ); // Int_t nupdates=align->GetNUpdates(); @@ -1693,12 +2007,12 @@ void AliTPCcalibTime::BookDistortionMaps(){ axisName[0] ="#Delta"; axisTitle[0] ="#Delta"; // - binsTrack[1] =20; - xminTrack[1] =-1.5; xmaxTrack[1]=1.5; + binsTrack[1] =44; + xminTrack[1] =-1.1; xmaxTrack[1]=1.1; axisName[1] ="tanTheta"; axisTitle[1] ="tan(#Theta)"; // - binsTrack[2] =90; + binsTrack[2] =180; xminTrack[2] =-TMath::Pi(); xmaxTrack[2]=TMath::Pi(); axisName[2] ="phi"; axisTitle[2] ="#phi"; @@ -1770,6 +2084,35 @@ void AliTPCcalibTime::BookDistortionMaps(){ } } } + // + // Book vertex: time histograms + // + Int_t binsVertex[2]={500, fTimeBins}; + Double_t aminVertex[2]={-5,fTimeStart}; + Double_t amaxVertex[2]={5, fTimeEnd}; + const char* hnames[12]={"TPCXAside", "TPCXCside","TPCXACdiff","TPCXAPCdiff", + "TPCYAside", "TPCYCside","TPCYACdiff","TPCYAPCdiff", + "TPCZAPCside", "TPCZAMCside","TPCZACdiff","TPCZAPCdiff"}; + const char* anames[12]={"x (cm) - A side ", "x (cm) - C side","#Delta_{x} (cm) - TPC-A-C","#Delta_{x} (cm) - TPC-Common", + "y (cm) - A side ", "y (cm) - C side","#Delta_{x} (cm) - TPC-A-C","#Delta_{y} (cm) - TPC-Common", + "z (cm)", "#Delta_{Z} (cm) A-C side","#Delta_{x} (cm) - TPC-A-C","#Delta_{Z} (cm) TPC-common"}; + for (Int_t ihis=0; ihis<12; ihis++) { + if (ihis>=8) aminVertex[0]=-20.; + if (ihis>=8) amaxVertex[0]=20.; + fTPCVertex[ihis]=new THnSparseF(hnames[ihis],hnames[ihis],2,binsVertex,aminVertex,amaxVertex); + fTPCVertex[ihis]->GetAxis(1)->SetTitle("Time"); + fTPCVertex[ihis]->GetAxis(0)->SetTitle(anames[ihis]); + } + + Int_t binsVertexC[2]={40, 300}; + Double_t aminVertexC[2]={-20,-30}; + Double_t amaxVertexC[2]={20,30}; + const char* hnamesC[5]={"TPCA_TPC","TPCC_TPC","TPCA_ITS","TPCC_ITS","TPC_ITS"}; + for (Int_t ihis=0; ihis<5; ihis++) { + fTPCVertexCorrelation[ihis]=new THnSparseF(hnamesC[ihis],hnamesC[ihis],2,binsVertexC,aminVertexC,amaxVertexC); + fTPCVertexCorrelation[ihis]->GetAxis(1)->SetTitle("z (cm)"); + fTPCVertexCorrelation[ihis]->GetAxis(0)->SetTitle("z (cm)"); + } } @@ -1778,6 +2121,7 @@ void AliTPCcalibTime::FillResHistoTPCCE(const AliExternalTrackParam * pTP // fill residual histograms pTPCOut-pTPCin - trac crossing CE // Histogram // + if (fMemoryMode<2) return; Double_t histoX[5]; Double_t xyz[3]; pTPCIn->GetXYZ(xyz); @@ -1823,6 +2167,7 @@ void AliTPCcalibTime::FillResHistoTPC(const AliESDtrack * pTrack){ // fill residual histograms pTPC - vertex // Histogram is filled only for primary tracks // + if (fMemoryMode<2) return; Double_t histoX[4]; const AliExternalTrackParam * pTPCIn = pTrack->GetInnerParam(); AliExternalTrackParam pTPCvertex(*(pTrack->GetInnerParam())); @@ -1854,6 +2199,7 @@ void AliTPCcalibTime::FillResHistoTPCTRD(const AliExternalTrackParam * pT // // fill resuidual histogram TPCout-TRDin // + if (fMemoryMode<2) return; Double_t histoX[4]; Double_t xyz[3]; pTPCOut->GetXYZ(xyz); @@ -1878,6 +2224,7 @@ void AliTPCcalibTime::FillResHistoTPCTOF(const AliExternalTrackParam * pT // // fill resuidual histogram TPCout-TOFin // track propagated to the TOF position + if (fMemoryMode<2) return; Double_t histoX[4]; Double_t xyz[3];