X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TPC%2FAliTPCcalibTime.cxx;h=a9aa6ef2d9610835b771540b3a4b5633ce2f1995;hp=0fc8d1cb3b269384c09be55bf4190715d10e2782;hb=0a8abcad68b30cf11bb0c83a3a42cfb293f5480d;hpb=817766d592506123e9d5f3ab7b72964a7d6d841c diff --git a/TPC/AliTPCcalibTime.cxx b/TPC/AliTPCcalibTime.cxx index 0fc8d1cb3b2..a9aa6ef2d96 100644 --- a/TPC/AliTPCcalibTime.cxx +++ b/TPC/AliTPCcalibTime.cxx @@ -22,107 +22,62 @@ Comments to be written here: AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04); -2. How to interpret results - -3. Simple example - - a) determine the required time range: - - AliXRDPROOFtoolkit tool; - TChain * chain = tool.MakeChain("pass2.txt","esdTree",0,6000); - chain->Draw("GetTimeStamp()") - - b) analyse calibration object on Proof in calibration train - - AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime", StartTimeStamp, EndTimeStamp, IntegrationTimeVdrift); - - c) plot results - .x ~/NimStyle.C - gSystem->Load("libANALYSIS"); - gSystem->Load("libTPCcalib"); - - TFile f("CalibObjectsTrain1.root"); - AliTPCcalibTime *calib = (AliTPCcalibTime *)f->Get("calibTime"); - calib->GetHistoDrift("all")->Projection(2,0)->Draw() - calib->GetFitDrift("all")->Draw("lp") - -4. Analysis using debug streamers. - - gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros"); - gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+") - AliXRDPROOFtoolkit tool; - AliXRDPROOFtoolkit::FilterList("timeitstpc.txt","* itstpc",1) - AliXRDPROOFtoolkit::FilterList("time.txt","* trackInfo",1) - AliXRDPROOFtoolkit::FilterList("timelaser.txt","* laserInfo",1) - - TChain * chainTPCITS = tool.MakeChainRandom("timeitstpc.txt.Good","itstpc",0,10000); - TChain * chainTPCTOF = tool.MakeChainRandom("time.txt","pointMatch",0,10000); - TChain * chainTime = tool.MakeChainRandom("time.txt.Good","trackInfo",0,10000); - TChain * chainLaser = tool.MakeChainRandom("timelaser.txt.Good","laserInfo",0,10000); - chainTime->Lookup(); - chainLaser->Lookup(); */ #include "Riostream.h" -#include "TChain.h" -#include "TTree.h" +#include "TDatabasePDG.h" +#include "TGraphErrors.h" #include "TH1F.h" -#include "TH2F.h" -#include "TH3F.h" #include "THnSparse.h" #include "TList.h" #include "TMath.h" -#include "TCanvas.h" -#include "TFile.h" -#include "TF1.h" +#include "TTimeStamp.h" +#include "TTree.h" #include "TVectorD.h" -#include "TProfile.h" -#include "TGraphErrors.h" -#include "TCanvas.h" -#include "AliTPCclusterMI.h" -#include "AliTPCseed.h" -#include "AliESDVertex.h" +//#include "TChain.h" +//#include "TFile.h" + +#include "AliDCSSensor.h" +#include "AliDCSSensorArray.h" #include "AliESDEvent.h" -#include "AliESDfriend.h" #include "AliESDInputHandler.h" -#include "AliAnalysisManager.h" - -#include "AliTracker.h" -#include "AliMagF.h" -#include "AliTPCCalROC.h" - +#include "AliESDVertex.h" +#include "AliESDfriend.h" #include "AliLog.h" - -#include "AliTPCcalibTime.h" #include "AliRelAlignerKalman.h" - -#include "TTreeStream.h" +#include "AliTPCCalROC.h" +#include "AliTPCParam.h" #include "AliTPCTracklet.h" -#include "TTimeStamp.h" #include "AliTPCcalibDB.h" #include "AliTPCcalibLaser.h" -#include "AliDCSSensorArray.h" -#include "AliDCSSensor.h" - -#include "TDatabasePDG.h" +#include "AliTPCcalibTime.h" +#include "AliTPCclusterMI.h" +#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 fCutMaxDz(25), // maximal distance in rfi ditection fCutTheta(0.03), // maximal distan theta fCutMinDir(-0.99), // direction vector products - fCutTracks(10), + fCutTracks(2500), + 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), @@ -135,10 +90,10 @@ AliTPCcalibTime::AliTPCcalibTime() fRunBins(0), fRunStart(0), fRunEnd(0) -// fBinsVdrift(fTimeBins,fPtBins,fVdriftBins), -// fXminVdrift(fTimeStart,fPtStart,fVdriftStart), -// fXmaxVdrift(fTimeEnd,fPtEnd,fVdriftEnd) { + // + // default constructor + // AliInfo("Default Constructor"); for (Int_t i=0;i<3;i++) { fHistVdriftLaserA[i]=0; @@ -147,21 +102,48 @@ AliTPCcalibTime::AliTPCcalibTime() for (Int_t i=0;i<10;i++) { fCosmiMatchingHisto[i]=0; } + // + for (Int_t i=0;i<5;i++) { + fResHistoTPCCE[i]=0; + fResHistoTPCITS[i]=0; + 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 fCutMaxDz(40), // maximal distance in rfi ditection fCutTheta(5*0.004644),// maximal distan theta fCutMinDir(-0.99), // direction vector products - fCutTracks(10), + fCutTracks(2500), + 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), @@ -175,6 +157,10 @@ AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t fRunStart(0), fRunEnd(0) { + // + // Non deafaul constructor - to be used in the Calibration setups + // + SetName(name); SetTitle(title); for (Int_t i=0;i<3;i++) { @@ -182,6 +168,15 @@ AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t fHistVdriftLaserC[i]=0; } + for (Int_t i=0;i<5;i++) { + fResHistoTPCCE[i]=0; + fResHistoTPCITS[i]=0; + fResHistoTPCTRD[i]=0; + fResHistoTPCTOF[i]=0; + fResHistoTPCvertex[i]=0; + } + + AliInfo("Non Default Constructor"); fTimeBins =(EndTime-StartTime)/deltaIntegrationTimeVdrift; fTimeStart =StartTime; //(((TObjString*)(mapGRP->GetValue("fAliceStartTime")))->GetString()).Atoi(); @@ -192,9 +187,9 @@ AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t fVdriftBins = 500; fVdriftStart= -0.1; fVdriftEnd = 0.1; - fRunBins = 100001; + fRunBins = 1000001; fRunStart = -1.5; - fRunEnd = 99999.5; + fRunEnd = 999999.5; Int_t binsVdriftLaser[4] = {fTimeBins , fPtBins , fVdriftBins*20, fRunBins }; Double_t xminVdriftLaser[4] = {fTimeStart, fPtStart, fVdriftStart , fRunStart}; @@ -239,13 +234,10 @@ AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t fAlignITSTPC = new TObjArray; //alignemnt array ITS TPC match fAlignTRDTPC = new TObjArray; //alignemnt array ITS TPC match fAlignTOFTPC = new TObjArray; //alignemnt array ITS TPC match - - // fArrayDz->AddLast(fHistVdriftLaserA[0]); -// fArrayDz->AddLast(fHistVdriftLaserA[1]); -// fArrayDz->AddLast(fHistVdriftLaserA[2]); -// fArrayDz->AddLast(fHistVdriftLaserC[0]); -// fArrayDz->AddLast(fHistVdriftLaserC[1]); -// fArrayDz->AddLast(fHistVdriftLaserC[2]); + fAlignITSTPC->SetOwner(kTRUE); + fAlignTRDTPC->SetOwner(kTRUE); + fAlignTOFTPC->SetOwner(kTRUE); + fCosmiMatchingHisto[0]=new TH1F("Cosmics matching","p0-all" ,100,-10*0.5356 ,10*0.5356 ); fCosmiMatchingHisto[1]=new TH1F("Cosmics matching","p1-all" ,100,-10*4.541 ,10*4.541 ); @@ -257,18 +249,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 ); -// Char_t nameHisto[3]={'p','0','\n'}; -// for (Int_t i=0;i<10;i++){ -// fCosmiMatchingHisto[i]=new TH1F("Cosmics matching",nameHisto,8192,0,0); -// nameHisto[1]++; -// if(i==4) nameHisto[1]='0'; -// } + 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(){ // - // Destructor + // 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]; @@ -291,30 +292,87 @@ AliTPCcalibTime::~AliTPCcalibTime(){ fCosmiMatchingHisto[i]=NULL; } } + + for (Int_t i=0;i<5;i++) { + delete fResHistoTPCCE[i]; + delete fResHistoTPCITS[i]; + delete fResHistoTPCTRD[i]; + delete fResHistoTPCTOF[i]; + delete fResHistoTPCvertex[i]; + fResHistoTPCCE[i]=0; + fResHistoTPCITS[i]=0; + fResHistoTPCTRD[i]=0; + fResHistoTPCTOF[i]=0; + fResHistoTPCvertex[i]=0; + } + + if (fTPCVertex[0]) { + for (Int_t i=0;i<12;i++) delete fTPCVertex[i]; + } + if (fTPCVertexCorrelation[0]) { + for (Int_t i=0;i<5;i++) delete fTPCVertexCorrelation[i]; + } + + fAlignITSTPC->SetOwner(kTRUE); + fAlignTRDTPC->SetOwner(kTRUE); + fAlignTOFTPC->SetOwner(kTRUE); + fAlignITSTPC->Delete(); fAlignTRDTPC->Delete(); fAlignTOFTPC->Delete(); + delete fAlignITSTPC; + delete fAlignTRDTPC; + delete fAlignTOFTPC; } -Bool_t AliTPCcalibTime::IsLaser(AliESDEvent */*event*/){ - return kTRUE; //More accurate creteria to be added -} -Bool_t AliTPCcalibTime::IsCosmics(AliESDEvent */*event*/){ - return kTRUE; //More accurate creteria to be added -} -Bool_t AliTPCcalibTime::IsBeam(AliESDEvent */*event*/){ - return kTRUE; //More accurate creteria to be added -} +// Bool_t AliTPCcalibTime::IsLaser(const AliESDEvent *const /*event*/) const{ +// // +// // Indicator is laser event not yet implemented - to be done using trigger info or event specie +// // +// return kTRUE; //More accurate creteria to be added +// } +// Bool_t AliTPCcalibTime::IsCosmics(const AliESDEvent *const /*event*/){ +// // +// // Indicator is cosmic event not yet implemented - to be done using trigger info or event specie +// // + +// return kTRUE; //More accurate creteria to be added +// } +// Bool_t AliTPCcalibTime::IsBeam(const AliESDEvent *const /*event*/) const{ +// // +// // Indicator is physic event not yet implemented - to be done using trigger info or event specie +// // + +// return kTRUE; //More accurate creteria to be added +// } void AliTPCcalibTime::ResetCurrent(){ + // + //ResetCurrent + // fDz=0; //Reset current dz } + + + void AliTPCcalibTime::Process(AliESDEvent *event){ + // + // main function to make calibration + // if(!event) return; - if (event->GetNumberOfTracks()<2) return; + if (event->GetNumberOfTracks()<2) return; + AliESDfriend *ESDfriend=static_cast(event->FindListObject("AliESDfriend")); + if (!ESDfriend) { + return; + } + if (ESDfriend->TestSkipBit()) return; + ResetCurrent(); - if(IsLaser (event)) ProcessLaser (event); - if(IsCosmics(event)) ProcessCosmic(event); - if(IsBeam (event)) ProcessBeam (event); + //if(IsLaser (event)) + ProcessLaser (event); + //if(IsCosmics(event)) + ProcessCosmic(event); + //if(IsBeam (event)) + ProcessBeam (event); } void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){ @@ -371,34 +429,12 @@ void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){ TTreeSRedirector *cstream = GetDebugStreamer(); if (cstream){ TTimeStamp tstamp(fTime); - Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0); - Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1); - Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0); - Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1); - Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0); - Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1); - TVectorD vecGoofie(20); - AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun); - if (goofieArray){ - for (Int_t isensor=0; isensorNumSensors();isensor++){ - AliDCSSensor *gsensor = goofieArray->GetSensor(isensor); - if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp); - } - } (*cstream)<<"laserInfo"<< "run="<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); + //if (isReject[0]==0) fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA); + //if (isReject[1]==0) fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC); + fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA); + fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC); + } } -void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){ +void AliTPCcalibTime::ProcessCosmic(const AliESDEvent *const event){ + // + // process Cosmic event - track matching A side C side + // if (!event) { Printf("ERROR: ESD not available"); return; @@ -495,12 +523,14 @@ void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){ // Track0 is choosen in upper TPC part // Track1 is choosen in lower TPC part // + const Int_t kMinClustersCross =30; + const Int_t kMinClusters =80; Int_t ntracks=event->GetNumberOfTracks(); if (ntracks==0) return; if (ntracks > fCutTracks) return; - 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=(AliESDfriend*)(((AliESDEvent*)event)->FindListObject("AliESDfriend")); TObjArray tpcSeeds(ntracks); Double_t vtxx[3]={0,0,0}; @@ -509,7 +539,11 @@ void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){ // // track loop // + TArrayI clusterSideA(ntracks); + TArrayI clusterSideC(ntracks); for (Int_t i=0;iGetTrack(i); const AliExternalTrackParam * trackIn = track->GetInnerParam(); @@ -517,13 +551,27 @@ void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){ if (!trackIn) continue; if (!trackOut) continue; - AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i); - if (friendTrack) ProcessAlignITS(track,friendTrack); + AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i); + if (!friendTrack) continue; + if (friendTrack) ProcessSame(track,friendTrack,event); + if (friendTrack) ProcessAlignITS(track,friendTrack,event,esdFriend); + if (friendTrack) ProcessAlignTRD(track,friendTrack); if (friendTrack) ProcessAlignTOF(track,friendTrack); 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) { + tpcSeeds.AddAt(seed,i); + Int_t nA=0, nC=0; + 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++; + } + clusterSideA[i]=nA; + clusterSideC[i]=nC; + } } if (ntracks<2) return; // @@ -544,7 +592,15 @@ void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){ //track 1 lower part if (!track1) continue; if (!track1->GetOuterParam()) continue; - if (track1->GetOuterParam()->GetAlpha()>0) continue; + if (track0->GetTPCNcls()+ track1->GetTPCNcls()< kMinClusters) continue; + Int_t nAC = TMath::Max( TMath::Min(clusterSideA[i], clusterSideC[j]), + TMath::Min(clusterSideC[i], clusterSideA[j])); + if (nACGetOuterParam()->GetAlpha()>0) continue; // Double_t d2[3]; track1->GetDirection(d2); @@ -559,8 +615,8 @@ void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){ // // conservative cuts - convergence to be guarantied // applying before track propagation - if (TMath::Abs(dist0+dist1)>fCutMaxD) continue; // distance to the 0,0 - if (dir>fCutMinDir) continue; // direction vector product + if (TMath::Abs(TMath::Abs(dist0)-TMath::Abs(dist1))>fCutMaxD) continue; // distance to the 0,0 + if (TMath::Abs(dir)GetParticle("e-")->Mass(),3,kTRUE); + AliTracker::PropagateTrackTo(¶m1,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE); // // Propagate rest to the 0,0 DCA - z should be ignored // @@ -602,10 +658,11 @@ void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){ delete hist; hist=NULL; - if(isSame || (isCross && isPair)){ - if (track0->GetTPCNcls() > 80) { + if((isSame) || (isCross && isPair)){ + if (track0->GetTPCNcls()+ track1->GetTPCNcls()> 80) { fDz = param0.GetZ() - param1.GetZ(); - if(track0->GetOuterParam()->GetZ()<0) fDz=-fDz; + Double_t sign=(nA0>nA1)? 1:-1; + fDz*=sign; TTimeStamp tstamp(fTime); Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0); Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1); @@ -652,6 +709,11 @@ void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){ "tr1.="<GetFiredTriggerClasses().Data()); + if (GetDebugLevel()>20) printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data()); } -void AliTPCcalibTime::ProcessBeam(AliESDEvent */*event*/){ +void AliTPCcalibTime::ProcessBeam(const AliESDEvent *const event){ + // + // Process beam data - calculates vartex + // from A side and C side + // Histogram the differences + // + 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; + if (ntracks > fCutTracks) 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(); iterator->Reset(); TString newName=name; @@ -720,7 +997,7 @@ THnSparse* AliTPCcalibTime::GetHistoDrift(const char* name){ THnSparse* newHist=new THnSparseF(newName,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift); THnSparse* addHist=NULL; while((addHist=(THnSparseF*)iterator->Next())){ - if(!addHist) continue; + // if(!addHist) continue; TString histName=addHist->GetName(); if(!histName.Contains(newName)) continue; addHist->Print(); @@ -729,11 +1006,18 @@ THnSparse* AliTPCcalibTime::GetHistoDrift(const char* name){ return newHist; } -TObjArray* AliTPCcalibTime::GetHistoDrift(){ +TObjArray* AliTPCcalibTime::GetHistoDrift() const +{ + // + // return array of histograms + // return fArrayDz; } TGraphErrors* AliTPCcalibTime::GetGraphDrift(const char* name){ + // + // Make a drift velocity (delta Z) graph + // THnSparse* histoDrift=GetHistoDrift(name); TGraphErrors* graphDrift=NULL; if(histoDrift){ @@ -750,6 +1034,9 @@ TGraphErrors* AliTPCcalibTime::GetGraphDrift(const char* name){ } TObjArray* AliTPCcalibTime::GetGraphDrift(){ + // + // make a array of drift graphs + // TObjArray* arrayGraphDrift=new TObjArray(); TIterator* iterator=fArrayDz->MakeIterator(); iterator->Reset(); @@ -759,6 +1046,9 @@ TObjArray* AliTPCcalibTime::GetGraphDrift(){ } AliSplineFit* AliTPCcalibTime::GetFitDrift(const char* name){ + // + // Make a fit AliSplinefit of drift velocity + // TGraph* graphDrift=GetGraphDrift(name); AliSplineFit* fitDrift=NULL; if(graphDrift && graphDrift->GetN()){ @@ -780,19 +1070,14 @@ AliSplineFit* AliTPCcalibTime::GetFitDrift(const char* name){ return fitDrift; } -//TObjArray* AliTPCcalibTime::GetFitDrift(){ -// TObjArray* arrayFitDrift=new TObjArray(); -// TIterator* iterator = fArrayDz->MakeIterator(); -// iterator->Reset(); -// THnSparse* addHist=NULL; -// while((addHist=(THnSparseF*)iterator->Next())) arrayFitDrift->AddLast(GetFitDrift(addHist->GetName())); -// return arrayFitDrift; -//} -Long64_t AliTPCcalibTime::Merge(TCollection *li) { +Long64_t AliTPCcalibTime::Merge(TCollection *const li) { + // + // Object specific merging procedure + // TIterator* iter = li->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()); @@ -804,13 +1089,63 @@ Long64_t AliTPCcalibTime::Merge(TCollection *li) { fHistVdriftLaserC[imeas]->Add(cal->GetHistVdriftLaserC(imeas)); } } + // + if (fTPCVertexCorrelation[0] && cal->fTPCVertexCorrelation[0]){ + for (Int_t imeas=0; imeas<5; imeas++){ + if (fTPCVertexCorrelation[imeas] && cal->fTPCVertexCorrelation[imeas]) fTPCVertexCorrelation[imeas]->Add(cal->fTPCVertexCorrelation[imeas]); + } + } + + if (fTPCVertex[0] && cal->fTPCVertex[0]) + for (Int_t imeas=0; imeas<12; imeas++){ + if (fTPCVertex[imeas] && cal->fTPCVertex[imeas]) fTPCVertex[imeas]->Add(cal->fTPCVertex[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 ((fMemoryMode>1) && cal->fResHistoTPCTRD[imeas]){ + if (fResHistoTPCTRD[imeas]) + fResHistoTPCTRD[imeas]->Add(cal->fResHistoTPCTRD[imeas]); + else + fResHistoTPCTRD[imeas]=(THnSparse*)cal->fResHistoTPCTRD[imeas]->Clone(); + } + // + 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(!addHist) continue; + if ((fMemoryMode>1)) while((addHist=(THnSparseF*)iterator->Next())){ + // if(!addHist) continue; addHist->Print(); THnSparse* localHist=(THnSparseF*)fArrayDz->FindObject(addHist->GetName()); if(!localHist){ @@ -819,28 +1154,39 @@ Long64_t AliTPCcalibTime::Merge(TCollection *li) { } localHist->Add(addHist); } -// TMap * addMap=cal->GetHistoDrift(); -// if(!addMap) return 0; -// TIterator* iterator = addMap->MakeIterator(); -// iterator->Reset(); -// TPair* addPair=0; -// while((addPair=(TPair *)(addMap->FindObject(iterator->Next())))){ -// THnSparse* addHist=dynamic_cast(addPair->Value()); -// if (!addHist) continue; -// addHist->Print(); -// THnSparse* localHist=dynamic_cast(fMapDz->GetValue(addHist->GetName())); -// if(!localHist){ -// localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift); -// fMapDz->Add(new TObjString(addHist->GetName()),localHist); -// } -// localHist->Add(addHist); -// } + for(Int_t i=0;i<10;i++) if (cal->GetCosmiMatchingHisto(i)) fCosmiMatchingHisto[i]->Add(cal->GetCosmiMatchingHisto(i)); + // + // Merge alignment + // + for (Int_t itype=0; itype<3; itype++){ + // + // + TObjArray *arr0= 0; + TObjArray *arr1= 0; + if (itype==0) {arr0=fAlignITSTPC; arr1=cal->fAlignITSTPC;} + if (itype==1) {arr0=fAlignTRDTPC; arr1=cal->fAlignTRDTPC;} + if (itype==2) {arr0=fAlignTOFTPC; arr1=cal->fAlignTOFTPC;} + if (!arr1) continue; + if (!arr0) arr0=new TObjArray(arr1->GetEntriesFast()); + if (arr1->GetEntriesFast()>arr0->GetEntriesFast()){ + arr0->Expand(arr1->GetEntriesFast()); + } + for (Int_t i=0;iGetEntriesFast(); i++){ + AliRelAlignerKalman *kalman1 = (AliRelAlignerKalman *)arr1->UncheckedAt(i); + AliRelAlignerKalman *kalman0 = (AliRelAlignerKalman *)arr0->UncheckedAt(i); + if (!kalman1) continue; + if (!kalman0) {arr0->AddAt(new AliRelAlignerKalman(*kalman1),i); continue;} + kalman0->SetRejectOutliers(kFALSE); + kalman0->Merge(kalman1); + } + } + } return 0; } -Bool_t AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){ +Bool_t AliTPCcalibTime::IsPair(const AliExternalTrackParam *tr0, const AliExternalTrackParam *tr1){ /* // 0. Same direction - OPOSITE - cutDir +cutT TCut cutDir("cutDir","dir<-0.99") @@ -877,109 +1223,381 @@ Bool_t AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackPara return kTRUE; } -Bool_t AliTPCcalibTime::IsCross(AliESDtrack *tr0, AliESDtrack *tr1){ - return tr0->GetOuterParam()->GetZ()*tr1->GetOuterParam()->GetZ()<0 && tr0->GetInnerParam()->GetZ()*tr1->GetInnerParam()->GetZ()<0 && tr0->GetOuterParam()->GetZ()*tr0->GetInnerParam()->GetZ()>0 && tr1->GetOuterParam()->GetZ()*tr1->GetInnerParam()->GetZ()>0; +Bool_t AliTPCcalibTime::IsCross(const AliESDtrack *const tr0, const AliESDtrack *const tr1){ + // + // check if the cosmic pair of tracks crossed A/C side + // + Bool_t result= tr0->GetOuterParam()->GetZ()*tr1->GetOuterParam()->GetZ()<0; + if (result==kFALSE) return result; + result=kTRUE; + return result; } -Bool_t AliTPCcalibTime::IsSame(AliESDtrack */*tr0*/, AliESDtrack */*tr1*/){ - // To be implemented - return kFALSE; +Bool_t AliTPCcalibTime::IsSame(const AliESDtrack *const tr0, const AliESDtrack *const tr1){ + // + // track crossing the CE + // 0. minimal number of clusters + // 1. Same sector +-1 + // 2. Inner and outer track param on opposite side + // 3. Outer and inner track parameter close each to other + // 3. + Bool_t result=kTRUE; + // + // inner and outer on opposite sides in z + // + const Int_t knclCut0 = 30; + const Double_t kalphaCut = 0.4; + // + // 0. minimal number of clusters + // + if (tr0->GetTPCNcls()GetTPCNcls()GetOuterParam()->GetAlpha()-tr1->GetOuterParam()->GetAlpha())>kalphaCut) return kFALSE; + // + // 2. Z crossing + // + if (tr0->GetOuterParam()->GetZ()*tr0->GetInnerParam()->GetZ()>0) result&=kFALSE; + if (tr1->GetOuterParam()->GetZ()*tr1->GetInnerParam()->GetZ()>0) result&=kFALSE; + if (result==kFALSE){ + return result; + } + // + // + const Double_t *p0I = tr0->GetInnerParam()->GetParameter(); + const Double_t *p1I = tr1->GetInnerParam()->GetParameter(); + const Double_t *p0O = tr0->GetOuterParam()->GetParameter(); + const Double_t *p1O = tr1->GetOuterParam()->GetParameter(); + // + if (TMath::Abs(p0I[0]-p1I[0])>fCutMaxD) result&=kFALSE; + if (TMath::Abs(p0I[1]-p1I[1])>fCutMaxDz) result&=kFALSE; + if (TMath::Abs(p0I[2]-p1I[2])>fCutTheta) result&=kFALSE; + if (TMath::Abs(p0I[3]-p1I[3])>fCutTheta) result&=kFALSE; + if (TMath::Abs(p0O[0]-p1O[0])>fCutMaxD) result&=kFALSE; + if (TMath::Abs(p0O[1]-p1O[1])>fCutMaxDz) result&=kFALSE; + if (TMath::Abs(p0O[2]-p1O[2])>fCutTheta) result&=kFALSE; + if (TMath::Abs(p0O[3]-p1O[3])>fCutTheta) result&=kFALSE; + if (result==kTRUE){ + result=kTRUE; // just to put break point here + } + return result; } -/* -chainDrift->Draw("p0.fP[0]+p1.fP[0]","isPair"); - mean ~-0.02 ~-0.03913 - RMS ~ 0.5 ~ 0.5356 --> 3 (fCutMaxD) - -chainDrift->Draw("p0.fP[1]-p1.fP[1]","isPair"); - mean ~ 0.1855 - RMS ~ 4.541 -->25 (fCutMaxDz) - -chainDrift->Draw("p1.fAlpha-p0.fAlpha+pi","isPair"); -//chainDrift->Draw("p1.fAlpha+p0.fAlpha","isPair"); -//chainDrift->Draw("p1.fP[2]-p0.fP[2]+pi","isPair"); -//chainDrift->Draw("p1.fP[2]+p0.fP[2]","isPair"); - mean ~ 0 ~ 0.001898 - RMS ~ 0.009 ~ 0.01134 --> 0.06 - -chainDrift->Draw("p0.fP[3]+p1.fP[3]","isPair"); - mean ~ 0.0013 ~ 0.001539 - RMS ~ 0.003 ~ 0.004644 --> 0.03 (fCutTheta) - -chainDrift->Draw("p1.fP[4]+p0.fP[4]>>his(100,-0.2,0.2)","isPair") - mean ~ 0.012 ~-0.0009729 - RMS ~ 0.036 ~ 0.03773 --> 0.2 -*/ +void AliTPCcalibTime::ProcessSame(const AliESDtrack *const track, AliESDfriendTrack *const friendTrack, const AliESDEvent *const event){ + // + // Process TPC tracks crossing CE + // + // 0. Select only track crossing the CE + // 1. Cut on the track length + // 2. Refit the the track on A and C side separatelly + // 3. Fill time histograms + const Int_t kMinNcl=100; + const Int_t kMinNclS=25; // minimul number of clusters on the sides + const Double_t pimass=TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); + const Double_t kMaxDy=1; // maximal distance in y + const Double_t kMaxDsnp=0.05; // maximal distance in snp + const Double_t kMaxDtheta=0.05; // maximal distance in theta + + if (!friendTrack->GetTPCOut()) return; + // + // 0. Select only track crossing the CE + // + if (track->GetInnerParam()->GetZ()*friendTrack->GetTPCOut()->GetZ()>0) return; + // + // 1. cut on track length + // + if (track->GetTPCNcls()GetCalibObject(l));++l) { + if ((seed=dynamic_cast(calibObject))) break; + } + if (!seed) return; + // + AliExternalTrackParam trackIn(*track->GetInnerParam()); + AliExternalTrackParam trackOut(*track->GetOuterParam()); + Double_t cov[3]={0.01,0.,0.01}; //use the same errors + Double_t xyz[3]={0,0.,0.0}; + Double_t bz =0; + Int_t nclIn=0,nclOut=0; + trackIn.ResetCovariance(1000.); + trackOut.ResetCovariance(1000.); + // + //2.a Refit inner + // + Int_t sideIn=0; + for (Int_t irow=0;irow<159;irow++) { + AliTPCclusterMI *cl=seed->GetClusterPointer(irow); + if (!cl) continue; + if (cl->GetX()<80) continue; + if (sideIn==0){ + if (cl->GetDetector()%36<18) sideIn=1; + if (cl->GetDetector()%36>=18) sideIn=-1; + } + if (sideIn== -1 && (cl->GetDetector()%36)<18) break; + if (sideIn== 1 &&(cl->GetDetector()%36)>=18) break; + Int_t sector = cl->GetDetector(); + Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha(); + if (TMath::Abs(dalpha)>0.01){ + if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break; + } + Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()}; + trackIn.GetXYZ(xyz); + bz = AliTracker::GetBz(xyz); + AliTracker::PropagateTrackToBxByBz(&trackIn,r[0],1.,pimass,kFALSE); + if (!trackIn.PropagateTo(r[0],bz)) break; + nclIn++; + trackIn.Update(&r[1],cov); + } + // + //2.b Refit outer + // + Int_t sideOut=0; + for (Int_t irow=159;irow>0;irow--) { + AliTPCclusterMI *cl=seed->GetClusterPointer(irow); + if (!cl) continue; + if (cl->GetX()<80) continue; + if (sideOut==0){ + if (cl->GetDetector()%36<18) sideOut=1; + if (cl->GetDetector()%36>=18) sideOut=-1; + if (sideIn==sideOut) break; + } + if (sideOut== -1 && (cl->GetDetector()%36)<18) break; + if (sideOut== 1 &&(cl->GetDetector()%36)>=18) break; + // + Int_t sector = cl->GetDetector(); + Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha(); + if (TMath::Abs(dalpha)>0.01){ + if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break; + } + Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()}; + trackOut.GetXYZ(xyz); + bz = AliTracker::GetBz(xyz); + AliTracker::PropagateTrackToBxByBz(&trackOut,r[0],1.,pimass,kFALSE); + if (!trackOut.PropagateTo(r[0],bz)) break; + nclOut++; + trackOut.Update(&r[1],cov); + } + trackOut.Rotate(trackIn.GetAlpha()); + Double_t meanX = (trackIn.GetX()+trackOut.GetX())*0.5; + trackIn.PropagateTo(meanX,bz); + trackOut.PropagateTo(meanX,bz); + if (TMath::Abs(trackIn.GetY()-trackOut.GetY())>kMaxDy) return; + if (TMath::Abs(trackIn.GetSnp()-trackOut.GetSnp())>kMaxDsnp) return; + if (TMath::Abs(trackIn.GetTgl()-trackOut.GetTgl())>kMaxDtheta) return; + if (TMath::Min(nclIn,nclOut)>kMinNclS){ + FillResHistoTPCCE(&trackIn,&trackOut); + } + TTreeSRedirector *cstream = GetDebugStreamer(); + if (cstream){ + TVectorD gxyz(3); + trackIn.GetXYZ(gxyz.GetMatrixArray()); + TTimeStamp tstamp(fTime); + (*cstream)<<"tpctpc"<< + "run="<Draw("(tIn.fP[1]-tOut.fP[1])*sign(-tIn.fP[3]):tIn.fP[3]","min(nclIn,nclOut)>30","") + if (TMath::Min(nclIn,nclOut)>kMinNclS){ + fDz = trackOut.GetZ()-trackIn.GetZ(); + if (trackOut.GetTgl()<0) fDz*=-1.; + TTimeStamp tstamp(fTime); + Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0); + Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1); + Double_t vecDrift[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()}; + // + // fill histograms per trigger class and itegrated + // + THnSparse* curHist=NULL; + for (Int_t itype=0; itype<2; itype++){ + TString name="MEAN_VDRIFT_CROSS_"; + if (itype==0){ + name+=event->GetFiredTriggerClasses(); + name.ToUpper(); + }else{ + name+="ALL"; + } + 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::ProcessAlignITS(AliESDtrack* track, AliESDfriendTrack *friendTrack){ +void AliTPCcalibTime::ProcessAlignITS(AliESDtrack *const track, const AliESDfriendTrack *const friendTrack, const AliESDEvent *const event, AliESDfriend *const esdFriend){ // - // Process track - // Update TPC-ITS alignment + // Process track - Update TPC-ITS alignment + // Updates: + // 0. Apply standartd cuts + // 1. Recalucluate the current statistic median/RMS + // 2. Apply median+-rms cut + // 3. Update kalman filter // - const Int_t kMinTPC = 80; - const Int_t kMinITS = 3; - const Double_t kMinZ = 10; - const Double_t kMaxDy = 2; - const Double_t kMaxAngle= 0.02; + const Int_t kMinTPC = 80; // minimal number of TPC cluster + 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.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. Apply standard cuts + // Int_t dummycl[1000]; - if (track->GetITSclusters(dummycl)GetTPCNcls()GetITSOut()) return; + if (!track->IsOn(AliESDtrack::kTPCrefit)) 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())); + // + 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; + // + // try to find standalone ITS track corresponing to the TPC if possible + // + Bool_t hasAlone=kFALSE; + Int_t ntracks=event->GetNumberOfTracks(); + for (Int_t i=0; iGetTrack(i); + if (!trackITS) continue; + if (trackITS->GetITSclusters(dummycl)GetTrack(i); + if (!itsfriendTrack) continue; + if (!itsfriendTrack->GetITSOut()) 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()); + AliTracker::PropagateTrackToBxByBz(&pITS,pTPC.GetX(),0.1,0.1,kFALSE); + if (TMath::Abs(pTPC.GetY()-pITS.GetY())> kMaxDy) continue; + if (TMath::Abs(pTPC.GetSnp()-pITS.GetSnp())> kMaxAngle) continue; + hasAlone=kTRUE; + } + if (!hasAlone) { + if (track->GetITSclusters(dummycl)kMaxDy) return; + if (TMath::Abs(pITS.GetSnp()-pTPC.GetSnp())>kMaxAngle) return; + if (TMath::Abs(pITS.GetTgl()-pTPC.GetTgl())>kMaxAngle) return; // + // 1. Update median and RMS info // - Int_t htime = fTime/3600; //time in hours - if (fAlignITSTPC->GetEntries()0)? 1.:-1.; + vecDelta[4]=0; + for (Int_t i=0;i<4;i++){ + vecDelta[i]=(pITS.GetParameter()[i]-pTPC.GetParameter()[i])*sign; + kgdP[i][kglast%kN]=vecDelta[i]; + } + kglast=(kglast+1); + Int_t entries=(kglast0.){ + vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i]; + vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals + } + } + // + // 2. Apply median+-rms cut + // + if (kglast<3) return; //median and RMS to be defined + if ( vecDeltaN[4]/4.>kSigmaCut) return; + // + // 3. Update alignment + // + Int_t htime = (fTime-fTimeKalmanBin/2)/fTimeKalmanBin; //time bins number + if (fAlignITSTPC->GetEntriesFast()Expand(htime*2+20); } AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignITSTPC->At(htime); if (!align){ + // make Alignment object if doesn't exist align=new AliRelAlignerKalman(); - align->SetOutRejSigma(2.); - //align->SetRejectOutliers(kFALSE); + 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); + + align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.); align->SetMagField(fMagF); fAlignITSTPC->AddAt(align,htime); } - pITS.Rotate(pTPC.GetAlpha()); - pITS.PropagateTo(pTPC.GetX(),fMagF); - if (TMath::Abs(pITS.GetY()-pTPC.GetY())>kMaxDy) return; - if (TMath::Abs(pITS.GetSnp()-pTPC.GetSnp())>kMaxAngle) return; - if (TMath::Abs(pITS.GetTgl()-pTPC.GetTgl())>kMaxAngle) return; 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); + if (TMath::Abs(dca[0])GetNUpdates(); - // align->SetRunNumber(fRun ); - static Int_t entry=-1; - entry++; - align->SetOutRejSigma(1.+1./Double_t(nupdates)); + align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates)); + align->SetRejectOutliers(kFALSE); TTreeSRedirector *cstream = GetDebugStreamer(); if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){ - TTimeStamp tstamp(fTime); - Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0); - Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1); - Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0); - Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1); - Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0); - Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1); - TVectorD vecGoofie(20); - AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun); - if (goofieArray){ - for (Int_t isensor=0; isensorNumSensors();isensor++){ - AliDCSSensor *gsensor = goofieArray->GetSensor(isensor); - if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp); - } - } TVectorD gpTPC(3), gdTPC(3); TVectorD gpITS(3), gdITS(3); pTPC.GetXYZ(gpTPC.GetMatrixArray()); @@ -992,123 +1610,610 @@ void AliTPCcalibTime::ProcessAlignITS(AliESDtrack* track, AliESDfriendTrack *fr "time="<GetTPCNcls()GetOuterParam()==0) return; - if (track->GetInnerParam()==0) return; - if (track->GetTOFsignal()<=0) return; + const Int_t kMinTPC = 80; // minimal number of TPC cluster + const Int_t kMinTRD = 50; // minimal number of TRD cluster + const Double_t kMinZ = 20; // maximal dz distance + const Double_t kMaxDy = 5.; // maximal dy distance + 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 + 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. Apply standard cuts + // + Int_t dummycl[1000]; + if (track->GetTRDclusters(dummycl)GetTPCNcls()GetTRDIn()) return; + if (!track->IsOn(AliESDtrack::kTRDrefit)) return; + if (!track->IsOn(AliESDtrack::kTRDout)) return; + if (!track->GetInnerParam()) return; + if (!friendTrack->GetTPCOut()) return; + // exclude crossing track + if (friendTrack->GetTPCOut()->GetZ()*track->GetInnerParam()->GetZ()<0) return; + if (TMath::Abs(track->GetInnerParam()->GetZ())GetTPCOut())); + AliTracker::PropagateTrackToBxByBz(&pTPC,kRefX,0.1,0.1,kFALSE); + AliExternalTrackParam pTRD(*(friendTrack->GetTRDIn())); + pTRD.Rotate(pTPC.GetAlpha()); + // pTRD.PropagateTo(pTPC.GetX(),fMagF); + AliTracker::PropagateTrackToBxByBz(&pTRD,pTPC.GetX(),0.1,0.1,kFALSE); + + ((Double_t*)pTRD.GetCovariance())[2]+=3.*3.; // increas sys errors + ((Double_t*)pTRD.GetCovariance())[9]+=0.1*0.1; // increse sys errors + + if (TMath::Abs(pTRD.GetY()-pTPC.GetY()) >kMaxDy) return; + if (TMath::Abs(pTRD.GetSnp()-pTPC.GetSnp())>kMaxAngle) return; + // if (TMath::Abs(pTRD.GetTgl()-pTPC.GetTgl())>kMaxAngle) return; + // + // 1. Update median and RMS info + // + TVectorD vecDelta(5),vecMedian(5), vecRMS(5); + TVectorD vecDeltaN(5); + Double_t sign=(pTRD.GetParameter()[1]>0)? 1.:-1.; + vecDelta[4]=0; + for (Int_t i=0;i<4;i++){ + vecDelta[i]=(pTRD.GetParameter()[i]-pTPC.GetParameter()[i])*sign; + kgdP[i][kglast%kN]=vecDelta[i]; + } + kglast=(kglast+1); + Int_t entries=(kglast0.){ + vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i]; + vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals + } + } + // + // 2. Apply median+-rms cut + // + if (kglast<3) return; //median and RMS to be defined + if ( vecDeltaN[4]/4.>kSigmaCut) return; + // + // 3. Update alignment + // + //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); + } + AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTRDTPC->At(htime); + if (!align){ + // make Alignment object if doesn't exist + 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); + align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.); + align->SetMagField(fMagF); + fAlignTRDTPC->AddAt(align,htime); + } + align->AddTrackParams(&pTRD,&pTPC); + //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); + if (TMath::Abs(dca[0])GetOuterParam())); - AliExternalTrackParam *param=0; + Int_t nupdates=align->GetNUpdates(); + align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates)); + align->SetRejectOutliers(kFALSE); + TTreeSRedirector *cstream = GetDebugStreamer(); + if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){ + TVectorD gpTPC(3), gdTPC(3); + TVectorD gpTRD(3), gdTRD(3); + pTPC.GetXYZ(gpTPC.GetMatrixArray()); + pTPC.GetDirection(gdTPC.GetMatrixArray()); + pTRD.GetXYZ(gpTRD.GetMatrixArray()); + pTRD.GetDirection(gdTRD.GetMatrixArray()); + (*cstream)<<"trdtpc"<< + "run="<GetTOFsignal()<=0) return; + if (!friendTrack->GetTPCOut()) return; + if (!track->GetInnerParam()) return; + if (!friendTrack->GetTPCOut()) return; const AliTrackPointArray *points=friendTrack->GetTrackPointArray(); if (!points) return; + AliExternalTrackParam pTPC(*(friendTrack->GetTPCOut())); + AliExternalTrackParam pTOF(pTPC); + Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass(); Int_t npoints = points->GetNPoints(); AliTrackPoint point; - //Double_t alpha= - Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass(); - TTreeSRedirector * cstream = GetDebugStreamer(); - // - // + Int_t naccept=0; // for (Int_t ipoint=0;ipointGetPoint(point,ipoint); Float_t xyz[3]; point.GetXYZ(xyz); Double_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); - if (r<300) continue; + if (r<350) continue; if (r>400) continue; - param=paramOut; - if (!param) continue; - AliTracker::PropagateTrackToBxByBz(param,r,mass,2.,kTRUE); - AliTracker::PropagateTrackToBxByBz(param,r,mass,0.1,kTRUE); - AliTrackPoint lpoint = point.Rotate(param->GetAlpha()); - param->PropagateTo(lpoint.GetX(),fMagF); - // - // - // this is ugly - we need AliCluster constructor - // - AliExternalTrackParam &pTPC=*param; - AliExternalTrackParam pTOF(*param); + AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,2.,kTRUE); + AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,0.1,kTRUE); + AliTrackPoint lpoint = point.Rotate(pTPC.GetAlpha()); + pTPC.PropagateTo(lpoint.GetX(),fMagF); + pTOF=pTPC; ((Double_t*)pTOF.GetParameter())[0] =lpoint.GetY(); ((Double_t*)pTOF.GetParameter())[1] =lpoint.GetZ(); - pTOF.ResetCovariance(20); - ((Double_t*)pTOF.GetCovariance())[0]+=3.*3.; - ((Double_t*)pTOF.GetCovariance())[2]+=3.*3.; - if (TMath::Abs(pTOF.GetY()-pTPC.GetY())>kMaxDy) continue; - if (TMath::Abs(pTOF.GetZ()-pTPC.GetZ())>kMaxDz) continue; - // - Int_t htime = fTime/3600; //time in hours - - if (fAlignTOFTPC->GetEntries()Expand(htime*2+20); - } - AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTOFTPC->At(htime); - if (!align){ - align=new AliRelAlignerKalman(); - align->SetOutRejSigma(2.); - //align->SetRejectOutliers(kFALSE); - align->SetMagField(fMagF); - fAlignTOFTPC->AddAt(align,htime); + ((Double_t*)pTOF.GetCovariance())[0]+=3.*3./12.; + ((Double_t*)pTOF.GetCovariance())[2]+=3.*3./12.; + ((Double_t*)pTOF.GetCovariance())[5]+=0.1*0.1; + ((Double_t*)pTOF.GetCovariance())[9]+=0.1*0.1; + naccept++; + } + if (naccept==0) return; // no tof match clusters + // + // 0. Apply standard cuts + // + if (track->GetTPCNcls()GetTPCOut()->GetZ()*track->GetInnerParam()->GetZ()<0) return; + // + if (TMath::Abs(pTOF.GetY()-pTPC.GetY()) >kMaxDy) return; + if (TMath::Abs(pTOF.GetSnp()-pTPC.GetSnp())>kMaxAngle) return; + if (TMath::Abs(pTOF.GetTgl()-pTPC.GetTgl())>kMaxAngle) return; + // + // 1. Update median and RMS info + // + TVectorD vecDelta(5),vecMedian(5), vecRMS(5); + TVectorD vecDeltaN(5); + Double_t sign=(pTOF.GetParameter()[1]>0)? 1.:-1.; + vecDelta[4]=0; + for (Int_t i=0;i<4;i++){ + vecDelta[i]=(pTOF.GetParameter()[i]-pTPC.GetParameter()[i])*sign; + kgdP[i][kglast%kN]=vecDelta[i]; + } + kglast=(kglast+1); + Int_t entries=(kglast0.){ + vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/(vecRMS[i]+1.); + vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals + if (TMath::Abs(vecDeltaN[i])>kSigmaCut) isOK=kFALSE; } - pTOF.Rotate(pTPC.GetAlpha()); - pTOF.PropagateTo(pTPC.GetX(),fMagF); - align->AddTrackParams(&pTOF,&pTPC); - align->SetTimeStamp(fTime); - Int_t nupdates=align->GetNUpdates(); - static Int_t entry=-1; - entry++; - align->SetOutRejSigma(1.+1./Double_t(nupdates)); - - // - // - if (cstream) { - (*cstream) << "pointMatch" << - "run="<GetEntriesFast()Expand(htime*2+20); + } + AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTOFTPC->At(htime); + if (!align){ + // make Alignment object if doesn't exist + 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); + align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.); + align->SetMagField(fMagF); + fAlignTOFTPC->AddAt(align,htime); + } + align->AddTrackParams(&pTOF,&pTPC); + Float_t dca[2],cov[3]; + track->GetImpactParameters(dca,cov); + if (TMath::Abs(dca[0])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(); + align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates)); + align->SetRejectOutliers(kFALSE); + TTreeSRedirector *cstream = GetDebugStreamer(); + if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){ + TVectorD gpTPC(3), gdTPC(3); + TVectorD gpTOF(3), gdTOF(3); + pTPC.GetXYZ(gpTPC.GetMatrixArray()); + pTPC.GetDirection(gdTPC.GetMatrixArray()); + pTOF.GetXYZ(gpTOF.GetMatrixArray()); + pTOF.GetDirection(gdTOF.GetMatrixArray()); + (*cstream)<<"toftpc"<< + "run="<GetAxis(ivar2)->SetName(axisName[ivar2].Data()); + fResHistoTPCCE[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data()); + if (ivar2<4){ + fResHistoTPCITS[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data()); + fResHistoTPCITS[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data()); + fResHistoTPCTRD[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data()); + fResHistoTPCTRD[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data()); + fResHistoTPCvertex[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data()); + fResHistoTPCvertex[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data()); + } } + } + // + // 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)"); + } +} +void AliTPCcalibTime::FillResHistoTPCCE(const AliExternalTrackParam * pTPCIn, const AliExternalTrackParam * pTPCOut ){ + // + // fill residual histograms pTPCOut-pTPCin - trac crossing CE + // Histogram + // + if (fMemoryMode<2) return; + Double_t histoX[5]; + Double_t xyz[3]; + pTPCIn->GetXYZ(xyz); + Double_t phi= TMath::ATan2(xyz[1],xyz[0]); + histoX[1]= pTPCIn->GetTgl(); + histoX[2]= phi; + histoX[3]= pTPCIn->GetSnp(); + histoX[4]= pTPCIn->GetX(); + AliExternalTrackParam lout(*pTPCOut); + lout.Rotate(pTPCIn->GetAlpha()); + lout.PropagateTo(pTPCIn->GetX(),fMagF); + // + for (Int_t ihisto=0; ihisto<5; ihisto++){ + histoX[0]=lout.GetParameter()[ihisto]-pTPCIn->GetParameter()[ihisto]; + fResHistoTPCCE[ihisto]->Fill(histoX); + } +} +void AliTPCcalibTime::FillResHistoTPCITS(const AliExternalTrackParam * pTPCIn, const AliExternalTrackParam * pITSOut ){ + // + // fill residual histograms pTPCIn-pITSOut + // Histogram is filled only for primary tracks + // + Double_t histoX[4]; + Double_t xyz[3]; + pTPCIn->GetXYZ(xyz); + Double_t phi= TMath::ATan2(xyz[1],xyz[0]); + histoX[1]= pTPCIn->GetTgl(); + histoX[2]= phi; + histoX[3]= pTPCIn->GetSnp(); + AliExternalTrackParam lits(*pITSOut); + lits.Rotate(pTPCIn->GetAlpha()); + lits.PropagateTo(pTPCIn->GetX(),fMagF); + // + for (Int_t ihisto=0; ihisto<5; ihisto++){ + histoX[0]=pTPCIn->GetParameter()[ihisto]-lits.GetParameter()[ihisto]; + fResHistoTPCITS[ihisto]->Fill(histoX); + } +} + +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())); + // + AliExternalTrackParam lits(*pTrack); + if (TMath::Abs(pTrack->GetY())>3) return; // beam pipe + pTPCvertex.Rotate(lits.GetAlpha()); + //pTPCvertex.PropagateTo(pTPCvertex->GetX(),fMagF); + AliTracker::PropagateTrackToBxByBz(&pTPCvertex,lits.GetX(),0.1,2,kFALSE); + AliTracker::PropagateTrackToBxByBz(&pTPCvertex,lits.GetX(),0.1,0.1,kFALSE); + Double_t xyz[3]; + pTPCIn->GetXYZ(xyz); + Double_t phi= TMath::ATan2(xyz[1],xyz[0]); + histoX[1]= pTPCIn->GetTgl(); + histoX[2]= phi; + histoX[3]= pTPCIn->GetSnp(); + // + Float_t dca[2], cov[3]; + pTrack->GetImpactParametersTPC(dca,cov); + for (Int_t ihisto=0; ihisto<5; ihisto++){ + histoX[0]=pTPCvertex.GetParameter()[ihisto]-lits.GetParameter()[ihisto]; + // if (ihisto<2) histoX[0]=dca[ihisto]; + fResHistoTPCvertex[ihisto]->Fill(histoX); } - delete paramOut; +} + + +void AliTPCcalibTime::FillResHistoTPCTRD(const AliExternalTrackParam * pTPCOut, const AliExternalTrackParam * pTRDIn ){ + // + // fill resuidual histogram TPCout-TRDin + // + if (fMemoryMode<2) return; + Double_t histoX[4]; + Double_t xyz[3]; + pTPCOut->GetXYZ(xyz); + Double_t phi= TMath::ATan2(xyz[1],xyz[0]); + histoX[1]= pTPCOut->GetTgl(); + histoX[2]= phi; + histoX[3]= pTPCOut->GetSnp(); + // + AliExternalTrackParam ltrd(*pTRDIn); + ltrd.Rotate(pTPCOut->GetAlpha()); + // ltrd.PropagateTo(pTPCOut->GetX(),fMagF); + AliTracker::PropagateTrackToBxByBz(<rd,pTPCOut->GetX(),0.1,0.1,kFALSE); + + for (Int_t ihisto=0; ihisto<5; ihisto++){ + histoX[0]=pTPCOut->GetParameter()[ihisto]-ltrd.GetParameter()[ihisto]; + fResHistoTPCTRD[ihisto]->Fill(histoX); + } + +} + +void AliTPCcalibTime::FillResHistoTPCTOF(const AliExternalTrackParam * pTPCOut, const AliExternalTrackParam * pTOFIn ){ + // + // fill resuidual histogram TPCout-TOFin + // track propagated to the TOF position + if (fMemoryMode<2) return; + Double_t histoX[4]; + Double_t xyz[3]; + + AliExternalTrackParam ltpc(*pTPCOut); + ltpc.Rotate(pTOFIn->GetAlpha()); + AliTracker::PropagateTrackToBxByBz(<pc,pTOFIn->GetX(),0.1,0.1,kFALSE); + // + ltpc.GetXYZ(xyz); + Double_t phi= TMath::ATan2(xyz[1],xyz[0]); + histoX[1]= ltpc.GetTgl(); + histoX[2]= phi; + histoX[3]= ltpc.GetSnp(); + // + for (Int_t ihisto=0; ihisto<2; ihisto++){ + histoX[0]=ltpc.GetParameter()[ihisto]-pTOFIn->GetParameter()[ihisto]; + fResHistoTPCTOF[ihisto]->Fill(histoX); + } + }