+void AliTPCcalibCosmic::MakeFitTree(TTree * treeInput, TTreeSRedirector *pcstream, const TObjArray * corrArray, Int_t step, Int_t run){
+ //
+ // Make fit tree
+ // refit the tracks with original points + corrected points for each correction
+ // Input:
+ // treeInput - tree with cosmic tracks
+ // pcstream - debug output
+
+ // Algorithm:
+ // Loop over pair of cosmic tracks:
+ // 1. Find trigger offset between cosmic event and "physic" trigger
+ // 2. Refit tracks with current transformation
+ // 3. Refit tracks using additional "primitive" distortion on top
+ // Best correction estimated as linear combination of corrections
+ // minimizing the observed distortions
+ // Observed distortions - matching in the y,z, snp, theta and 1/pt
+ //
+ const Double_t kResetCov=20.;
+ const Double_t kMaxDelta[5]={1,40,0.03,0.01,0.2};
+ const Double_t kSigma=2.;
+ const Double_t kMaxTime=1050;
+ const Double_t kMaxSnp=0.8;
+ Int_t ncorr=corrArray->GetEntries();
+ AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+ AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
+ AliGRPObject* grp = AliTPCcalibDB::Instance()->GetGRP(run);
+ Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
+ transform->SetCurrentRun(run);
+ transform->SetCurrentTimeStamp(TMath::Nint(time));
+ Double_t covar[15];
+ for (Int_t i=0;i<15;i++) covar[i]=0;
+ covar[0]=kSigma*kSigma;
+ covar[2]=kSigma*kSigma;
+ covar[5]=kSigma*kSigma/Float_t(150*150);
+ covar[9]=kSigma*kSigma/Float_t(150*150);
+ covar[14]=0.2*0.2;
+ Double_t *distortions = new Double_t[ncorr+1];
+
+ AliESDtrack *track0=new AliESDtrack;
+ AliESDtrack *track1=new AliESDtrack;
+ AliESDfriendTrack *ftrack0=new AliESDfriendTrack;
+ AliESDfriendTrack *ftrack1=new AliESDfriendTrack;
+ treeInput->SetBranchAddress("t0.",&track0);
+ treeInput->SetBranchAddress("t1.",&track1);
+ treeInput->SetBranchAddress("ft0.",&ftrack0);
+ treeInput->SetBranchAddress("ft1.",&ftrack1);
+ Int_t entries= treeInput->GetEntries();
+ for (Int_t i=0; i<entries; i+=step){
+ treeInput->GetEntry(i);
+ if (i%20==0) printf("%d\n",i);
+ if (!ftrack0->GetTPCOut()) continue;
+ if (!ftrack1->GetTPCOut()) continue;
+ AliTPCseed *seed0=0;
+ AliTPCseed *seed1=0;
+ TObject *calibObject;
+ for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
+ if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
+ }
+ for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
+ if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
+ }
+ if (!seed0) continue;
+ if (!seed1) continue;
+ if (TMath::Abs(seed0->GetSnp())>kMaxSnp) continue;
+ if (TMath::Abs(seed1->GetSnp())>kMaxSnp) continue;
+ //
+ //
+ Int_t nclA0=0, nclC0=0; // number of clusters
+ Int_t nclA1=0, nclC1=0; // number of clusters
+ Int_t ncl0=0,ncl1=0;
+ Double_t rmin0=300, rmax0=-300; // variables to estimate the time0 - trigger offset
+ Double_t rmin1=300, rmax1=-300;
+ Double_t tmin0=2000, tmax0=-2000;
+ Double_t tmin1=2000, tmax1=-2000;
+ //
+ //
+ // calculate trigger offset usig "missing clusters"
+ for (Int_t irow=0; irow<159; irow++){
+ AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
+ if (cluster0 &&cluster0->GetX()>10){
+ if (cluster0->GetX()<rmin0) { rmin0=cluster0->GetX(); tmin0=cluster0->GetTimeBin();}
+ if (cluster0->GetX()>rmax0) { rmax0=cluster0->GetX(); tmax0=cluster0->GetTimeBin();}
+ ncl0++;
+ if (cluster0->GetDetector()%36<18) nclA0++;
+ if (cluster0->GetDetector()%36>=18) nclC0++;
+ }
+ AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
+ if (cluster1&&cluster1->GetX()>10){
+ if (cluster1->GetX()<rmin1) { rmin1=cluster1->GetX(); tmin1=cluster1->GetTimeBin();}
+ if (cluster1->GetX()>rmax1) { rmax1=cluster1->GetX(); tmax1=cluster1->GetTimeBin();}
+ ncl1++;
+ if (cluster1->GetDetector()%36<18) nclA1++;
+ if (cluster1->GetDetector()%36>=18) nclC1++;
+ }
+ }
+ Int_t cosmicType=0; // types of cosmic topology
+ if ((nclA0>nclC0) && (nclA1>nclC1)) cosmicType=0; // AA side
+ if ((nclA0<nclC0) && (nclA1<nclC1)) cosmicType=1; // CC side
+ if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=2; // AC side
+ if ((nclA0<nclC0) && (nclA1>nclC1)) cosmicType=3; // CA side
+ //if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=6; // AC side out of time
+ //if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=7; // CA side out of time
+ //
+ Double_t deltaTime = 0; // correction for trigger not in time - equalizing the time arival
+ if ((cosmicType>1)&&TMath::Abs(track1->GetZ()-track0->GetZ())>4){
+ cosmicType+=4;
+ deltaTime=0.5*(track1->GetZ()-track0->GetZ())/param->GetZWidth();
+ if (nclA0>nclC0) deltaTime*=-1; // if A side track
+ }
+ //
+ TVectorD vectorDT(8);
+ Int_t crossCounter=0;
+ Double_t deltaTimeCross = AliTPCcalibCosmic::GetDeltaTime(rmin0, rmax0, rmin1, rmax1, tmin0, tmax0, tmin1, tmax1, TMath::Abs(track0->GetY()),vectorDT);
+ Bool_t isOKTrigger=kTRUE;
+ for (Int_t ic=0; ic<6;ic++) {
+ if (TMath::Abs(vectorDT[ic])>0) {
+ if (vectorDT[ic]+vectorDT[6]<0) isOKTrigger=kFALSE;
+ if (vectorDT[ic]+vectorDT[7]>kMaxTime) isOKTrigger=kFALSE;
+ if (isOKTrigger){
+ crossCounter++;
+ }
+ }
+ }
+ Double_t deltaTimeCluster=deltaTime;
+ if ((cosmicType==0 || cosmicType==1) && crossCounter>0){
+ deltaTimeCluster=deltaTimeCross;
+ cosmicType+=8;
+ }
+ if (nclA0*nclC0>0 || nclA1*nclC1>0) cosmicType+=16; // mixed A side C side - bad for visualization
+ //
+ // Apply current transformation
+ //
+ //
+ for (Int_t irow=0; irow<159; irow++){
+ AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
+ if (cluster0 &&cluster0->GetX()>10){
+ Double_t x0[3]={cluster0->GetRow(),cluster0->GetPad(),cluster0->GetTimeBin()+deltaTimeCluster};
+ Int_t index0[1]={cluster0->GetDetector()};
+ transform->Transform(x0,index0,0,1);
+ cluster0->SetX(x0[0]);
+ cluster0->SetY(x0[1]);
+ cluster0->SetZ(x0[2]);
+ //
+ }
+ AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
+ if (cluster1&&cluster1->GetX()>10){
+ Double_t x1[3]={cluster1->GetRow(),cluster1->GetPad(),cluster1->GetTimeBin()+deltaTimeCluster};
+ Int_t index1[1]={cluster1->GetDetector()};
+ transform->Transform(x1,index1,0,1);
+ cluster1->SetX(x1[0]);
+ cluster1->SetY(x1[1]);
+ cluster1->SetZ(x1[2]);
+ }
+ }
+ //
+ //
+ Double_t alpha=track0->GetAlpha(); // rotation frame
+ Double_t cos = TMath::Cos(alpha);
+ Double_t sin = TMath::Sin(alpha);
+ Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
+ AliExternalTrackParam btrack0=*(ftrack0->GetTPCOut());
+ AliExternalTrackParam btrack1=*(ftrack1->GetTPCOut());
+ btrack0.Rotate(alpha);
+ btrack1.Rotate(alpha);
+ // change the sign for track 1
+ Double_t* par1=(Double_t*)btrack0.GetParameter();
+ par1[3]*=-1;
+ par1[4]*=-1;
+ btrack0.AddCovariance(covar);
+ btrack1.AddCovariance(covar);
+ btrack0.ResetCovariance(kResetCov);
+ btrack1.ResetCovariance(kResetCov);
+ Bool_t isOK=kTRUE;
+ Bool_t isOKT=kTRUE;
+ TObjArray tracks0(ncorr+1);
+ TObjArray tracks1(ncorr+1);
+ //
+ Double_t dEdx0Tot=seed0->CookdEdxAnalytical(0.02,0.6,kTRUE);
+ Double_t dEdx1Tot=seed1->CookdEdxAnalytical(0.02,0.6,kTRUE);
+ Double_t dEdx0Max=seed0->CookdEdxAnalytical(0.02,0.6,kFALSE);
+ Double_t dEdx1Max=seed1->CookdEdxAnalytical(0.02,0.6,kFALSE);
+ //if (TMath::Abs((dEdx0Max+1)/(dEdx0Tot+1)-1.)>0.1) isOK=kFALSE;
+ //if (TMath::Abs((dEdx1Max+1)/(dEdx1Tot+1)-1.)>0.1) isOK=kFALSE;
+ ncl0=0; ncl1=0;
+ for (Int_t icorr=-1; icorr<ncorr; icorr++){
+ AliExternalTrackParam rtrack0=btrack0;
+ AliExternalTrackParam rtrack1=btrack1;
+ AliTPCCorrection *corr = 0;
+ if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
+ //
+ for (Int_t irow=159; irow>0; irow--){
+ AliTPCclusterMI *cluster=seed0->GetClusterPointer(irow);
+ if (!cluster) continue;
+ if (!isOKT) break;
+ Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
+ transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD); // transform to global
+ Float_t r[3]={rD[0],rD[1],rD[2]};
+ if (corr){
+ corr->DistortPoint(r, cluster->GetDetector());
+ }
+ //
+ Double_t cov[3]={0.01,0.,0.01};
+ Double_t lx =cos*r[0]+sin*r[1];
+ Double_t ly =-sin*r[0]+cos*r[1];
+ rD[1]=ly; rD[0]=lx; rD[2]=r[2]; //transform to track local
+ if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, lx,mass,1.,kFALSE)) isOKT=kFALSE;;
+ if (!rtrack0.Update(&rD[1],cov)) isOKT =kFALSE;
+ if (icorr<0) ncl0++;
+ }
+ //
+ for (Int_t irow=159; irow>0; irow--){
+ AliTPCclusterMI *cluster=seed1->GetClusterPointer(irow);
+ if (!cluster) continue;
+ if (!isOKT) break;
+ Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
+ transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD);
+ Float_t r[3]={rD[0],rD[1],rD[2]};
+ if (corr){
+ corr->DistortPoint(r, cluster->GetDetector());
+ }
+ //
+ Double_t cov[3]={0.01,0.,0.01};
+ Double_t lx =cos*r[0]+sin*r[1];
+ Double_t ly =-sin*r[0]+cos*r[1];
+ rD[1]=ly; rD[0]=lx; rD[2]=r[2];
+ if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, lx,mass,1.,kFALSE)) isOKT=kFALSE;
+ if (!rtrack1.Update(&rD[1],cov)) isOKT=kFALSE;
+ if (icorr<0) ncl1++;
+ }
+ if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,10.,kFALSE)) isOKT=kFALSE;
+ if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,10.,kFALSE)) isOKT=kFALSE;
+ if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,1.,kFALSE)) isOKT=kFALSE;
+ if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,1.,kFALSE)) isOKT=kFALSE;
+ const Double_t *param0=rtrack0.GetParameter();
+ const Double_t *param1=rtrack1.GetParameter();
+ for (Int_t ipar=0; ipar<4; ipar++){
+ if (TMath::Abs(param1[ipar]-param0[ipar])>kMaxDelta[ipar]) isOK=kFALSE;
+ }
+ tracks0.AddAt(rtrack0.Clone(), icorr+1);
+ tracks1.AddAt(rtrack1.Clone(), icorr+1);
+ AliExternalTrackParam out0=*(ftrack0->GetTPCOut());
+ AliExternalTrackParam out1=*(ftrack1->GetTPCOut());
+ Int_t nentries=TMath::Min(ncl0,ncl1);
+
+ if (icorr<0) {
+ (*pcstream)<<"cosmic"<<
+ "isOK="<<isOK<< // correct all propagation update and also residuals
+ "isOKT="<<isOKT<< // correct all propagation update
+ "isOKTrigger="<<isOKTrigger<< // correct? estimate of trigger offset
+ "id="<<cosmicType<<
+ //
+ //
+ "cross="<<crossCounter<<
+ "vDT.="<<&vectorDT<<
+ //
+ "dTime="<<deltaTime<< // delta time using the A-c side cross
+ "dTimeCross="<<deltaTimeCross<< // delta time using missing clusters
+ //
+ "dEdx0Max="<<dEdx0Max<<
+ "dEdx0Tot="<<dEdx0Tot<<
+ "dEdx1Max="<<dEdx1Max<<
+ "dEdx1Tot="<<dEdx1Tot<<
+ //
+ "track0.="<<track0<< // original track 0
+ "track1.="<<track1<< // original track 1
+ "out0.="<<&out0<< // outer track 0
+ "out1.="<<&out1<< // outer track 1
+ "rtrack0.="<<&rtrack0<< // refitted track with current transform
+ "rtrack1.="<<&rtrack1<< //
+ "ncl0="<<ncl0<<
+ "ncl1="<<ncl1<<
+ "entries="<<nentries<< // number of clusters
+ "\n";
+ }
+ }
+ //
+
+ if (isOK){
+ Int_t nentries=TMath::Min(ncl0,ncl1);
+ for (Int_t ipar=0; ipar<5; ipar++){
+ for (Int_t icorr=-1; icorr<ncorr; icorr++){
+ AliTPCCorrection *corr = 0;
+ if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
+ //
+ AliExternalTrackParam *param0=(AliExternalTrackParam *) tracks0.At(icorr+1);
+ AliExternalTrackParam *param1=(AliExternalTrackParam *) tracks1.At(icorr+1);
+ distortions[icorr+1]=param1->GetParameter()[ipar]-param0->GetParameter()[ipar];
+ if (icorr>=0){
+ distortions[icorr+1]-=distortions[0];
+ }
+ //
+ if (icorr<0){
+ Double_t bz=AliTrackerBase::GetBz();
+ Double_t gxyz[3];
+ param0->GetXYZ(gxyz);
+ Int_t dtype=20;
+ Double_t theta=param0->GetParameter()[3];
+ Double_t phi = alpha;
+ Double_t snp = track0->GetInnerParam()->GetSnp();
+ Double_t mean= distortions[0];
+ Int_t index = param0->GetIndex(ipar,ipar);
+ Double_t rms=TMath::Sqrt(param1->GetCovariance()[index]+param1->GetCovariance()[index]);
+ if (crossCounter<1) rms*=1;
+ Double_t sector=9*phi/TMath::Pi();
+ Double_t dsec = sector-TMath::Nint(sector+0.5);
+ Double_t gx=gxyz[0],gy=gxyz[1],gz=gxyz[2];
+ Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
+ Double_t dRrec=0;
+ // Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5;
+ Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5;
+
+ (*pcstream)<<"fit"<< // dump valus for fit
+ "run="<<run<< //run number
+ "bz="<<bz<< // magnetic filed used
+ "dtype="<<dtype<< // detector match type 20
+ "ptype="<<ipar<< // parameter type
+ "theta="<<theta<< // theta
+ "phi="<<phi<< // phi
+ "snp="<<snp<< // snp
+ "mean="<<mean<< // mean dist value
+ "rms="<<rms<< // rms
+ "sector="<<sector<<
+ "dsec="<<dsec<<
+ //
+ "refX="<<refX<< // reference radius
+ "gx="<<gx<< // global position
+ "gy="<<gy<< // global position
+ "gz="<<gz<< // global position
+ "dRrec="<<dRrec<< // delta Radius in reconstruction
+ "pt="<<pt<< //1/pt
+ "id="<<cosmicType<< //type of the cosmic used
+ "entries="<<nentries;// number of entries in bin
+ (*pcstream)<<"cosmicDebug"<<
+ "p0.="<<param0<< // dump distorted track 0
+ "p1.="<<param1; // dump distorted track 1
+ }
+ if (icorr>=0){
+ (*pcstream)<<"fit"<<
+ Form("%s=",corr->GetName())<<distortions[icorr+1]; // dump correction value
+ (*pcstream)<<"cosmicDebug"<<
+ Form("%s=",corr->GetName())<<distortions[icorr+1]<< // dump correction value
+ Form("%sp0.=",corr->GetName())<<param0<< // dump distorted track 0
+ Form("%sp1.=",corr->GetName())<<param1; // dump distorted track 1
+ }
+ } //loop corrections
+ (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
+ (*pcstream)<<"cosmicDebug"<<"isOK="<<isOK<<"\n";
+ } //loop over parameters
+ } // dump results
+ }//loop tracks
+ delete [] distortions;
+}