3 /**************************************************************************
4 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
6 * Author: The ALICE Off-line Project. *
7 * Contributors are mentioned in the code where appropriate. *
9 * Permission to use, copy, modify and distribute this software and its *
10 * documentation strictly for non-commercial purposes is hereby granted *
11 * without fee, provided that the above copyright notice appears in all *
12 * copies and that both the copyright notice and this permission notice *
13 * appear in the supporting documentation. The authors make no claims *
14 * about the suitability of this software for any purpose. It is *
15 * provided "as is" without express or implied warranty. *
16 **************************************************************************/
19 Comments to be written here:
20 1. What do we calibrate.
21 2. How to interpret results
23 4. Analysis using debug streamers.
28 // To make cosmic scan the user interaction neccessary
35 #include "Riostream.h"
45 #include "THnSparse.h"
46 #include "TDatabasePDG.h"
48 #include "AliTPCclusterMI.h"
49 #include "AliTPCseed.h"
50 #include "AliESDVertex.h"
51 #include "AliESDEvent.h"
52 #include "AliESDfriend.h"
53 #include "AliESDInputHandler.h"
54 #include "AliAnalysisManager.h"
56 #include "AliTracker.h"
58 #include "AliTPCCalROC.h"
59 #include "AliTPCParam.h"
62 #include "AliTPCcalibCosmic.h"
63 #include "TTreeStream.h"
64 #include "AliTPCTracklet.h"
65 //#include "AliESDcosmic.h"
66 #include "AliRecoParam.h"
67 #include "AliMultiplicity.h"
68 #include "AliTPCTransform.h"
69 #include "AliTPCcalibDB.h"
70 #include "AliTPCseed.h"
71 #include "AliGRPObject.h"
72 #include "AliTPCCorrection.h"
73 ClassImp(AliTPCcalibCosmic)
76 AliTPCcalibCosmic::AliTPCcalibCosmic()
85 fCutMaxD(5), // maximal distance in rfi ditection
86 fCutMaxDz(40), // maximal distance in z ditection
87 fCutTheta(0.03), // maximal distan theta
88 fCutMinDir(-0.99), // direction vector products
89 fCosmicTree(0) // tree with cosmic data
92 // CONSTRUCTOR - SEE COMMENTS ABOVE
94 AliInfo("Default Constructor");
95 for (Int_t ihis=0; ihis<6;ihis++){
99 for (Int_t ihis=0; ihis<4;ihis++){
100 fHistodEdxMax[ihis] =0;
101 fHistodEdxTot[ihis] =0;
106 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
115 fCutMaxD(5), // maximal distance in rfi ditection
116 fCutMaxDz(40), // maximal distance in z ditection
117 fCutTheta(0.03), // maximal distan theta
118 fCutMinDir(-0.99), // direction vector products
119 fCosmicTree(0) // tree with cosmic data
122 // cONSTRUCTOR - SEE COMENTS ABOVE
127 fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
128 fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
129 fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
130 fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
131 fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
133 fDeDxMIP = new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
135 AliInfo("Non Default Constructor");
139 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
143 for (Int_t ihis=0; ihis<6;ihis++){
144 delete fHistoDelta[ihis];
145 delete fHistoPull[ihis];
147 for (Int_t ihis=0; ihis<4;ihis++){
148 delete fHistodEdxTot[ihis];
149 delete fHistodEdxMax[ihis];
152 delete fHistNTracks; // histogram showing number of ESD tracks per event
153 delete fClusters; // histogram showing the number of clusters per track
154 delete fModules; // 2d histogram of tracks which are propagated to the ACORDE scintillator array
155 delete fHistPt; // Pt histogram of reconstructed tracks
156 delete fDeDx; // dEdx spectrum showing the different particle types
157 delete fDeDxMIP; // TPC signal close to the MIP region of muons 0.4 < p < 0.45 GeV
161 void AliTPCcalibCosmic::Init(){
164 // Make performance histograms
167 // tracking performance bins
168 // 0 - delta of interest
169 // 1 - min (track0, track1) number of clusters
170 // 2 - R - vertex radius
172 // 4 - P2 - snp(phi) at inner wall of TPC
173 // 5 - P3 - tan(theta) at inner wall of TPC
174 // 6 - P4 - 1/pt mean
177 // 9 - is corss indicator
179 Double_t xminTrack[10], xmaxTrack[10];
181 TString axisName[10];
184 axisName[0] ="#Delta";
187 xminTrack[1] =80; xmaxTrack[1]=160;
188 axisName[1] ="N_{cl}";
191 xminTrack[2] =0; xmaxTrack[2]=90; //
192 axisName[2] ="dca_{r} (cm)";
195 xminTrack[3] =-250; xmaxTrack[3]=250; //
196 axisName[3] ="z (cm)";
199 xminTrack[4] =-0.8; xmaxTrack[4]=0.8; //
200 axisName[4] ="sin(#phi)";
203 xminTrack[5] =-1; xmaxTrack[5]=1; //
204 axisName[5] ="tan(#theta)";
207 xminTrack[6] =-2; xmaxTrack[6]=2; //
208 axisName[6] ="1/pt (1/GeV)";
211 xminTrack[7] =1; xmaxTrack[7]=1000; //
212 axisName[7] ="pt (GeV)";
215 xminTrack[8] =0; xmaxTrack[8]=TMath::Pi(); //
216 axisName[8] ="alpha";
219 xminTrack[9] =-0.1; xmaxTrack[9]=2.1; //
220 axisName[9] ="cross";
223 xminTrack[0] =-1; xmaxTrack[0]=1; //
224 fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", ndim, binsTrack,xminTrack, xmaxTrack);
225 xminTrack[0] =-5; xmaxTrack[0]=5; //
226 fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
229 xminTrack[0] =-1; xmaxTrack[0]=1; //
230 fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", ndim, binsTrack,xminTrack, xmaxTrack);
231 xminTrack[0] =-5; xmaxTrack[0]=5; //
232 fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
235 xminTrack[0] =-10; xmaxTrack[0]=10; //
236 fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", ndim, binsTrack,xminTrack, xmaxTrack);
237 xminTrack[0] =-5; xmaxTrack[0]=5; //
238 fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
241 xminTrack[0] =-10; xmaxTrack[0]=10; //
242 fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", ndim, binsTrack,xminTrack, xmaxTrack);
243 xminTrack[0] =-5; xmaxTrack[0]=5; //
244 fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
247 xminTrack[0] =-0.2; xmaxTrack[0]=0.2; //
248 fHistoDelta[4] = new THnSparseS("#Delta_{1/pt} (1/GeV)","#Delta_{1/pt} (1/GeV)", ndim, binsTrack,xminTrack, xmaxTrack);
249 xminTrack[0] =-5; xmaxTrack[0]=5; //
250 fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
254 xminTrack[0] =-0.5; xmaxTrack[0]=0.5; //
255 fHistoDelta[5] = new THnSparseS("#Delta_{pt}/p_{t}","#Delta_{pt}/p_{t}", ndim, binsTrack,xminTrack, xmaxTrack);
256 xminTrack[0] =-5; xmaxTrack[0]=5; //
257 fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
260 for (Int_t idedx=0;idedx<4;idedx++){
261 xminTrack[0] =0.5; xmaxTrack[0]=1.5; //
263 xminTrack[1] =10; xmaxTrack[1]=160;
265 fHistodEdxMax[idedx] = new THnSparseS(Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
266 Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
267 ndim, binsTrack,xminTrack, xmaxTrack);
268 fHistodEdxTot[idedx] = new THnSparseS(Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
269 Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
270 ndim, binsTrack,xminTrack, xmaxTrack);
275 for (Int_t ivar=0;ivar<6;ivar++){
276 for (Int_t ivar2=0;ivar2<ndim;ivar2++){
277 fHistoDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
278 fHistoDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
279 fHistoPull[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
280 fHistoPull[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
281 BinLogX(fHistoDelta[ivar],7);
282 BinLogX(fHistoPull[ivar],7);
284 fHistodEdxMax[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
285 fHistodEdxMax[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
286 fHistodEdxTot[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
287 fHistodEdxTot[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
288 BinLogX(fHistodEdxMax[ivar],7);
289 BinLogX(fHistodEdxTot[ivar],7);
296 void AliTPCcalibCosmic::Add(const AliTPCcalibCosmic* cosmic){
298 // merge the content of the cosmic componentnts
300 for (Int_t ivar=0; ivar<6;ivar++){
301 if (fHistoDelta[ivar] && cosmic->fHistoDelta[ivar]){
302 fHistoDelta[ivar]->Add(cosmic->fHistoDelta[ivar]);
304 if (fHistoPull[ivar] && cosmic->fHistoPull[ivar]){
305 fHistoPull[ivar]->Add(cosmic->fHistoPull[ivar]);
308 for (Int_t ivar=0; ivar<4;ivar++){
309 if (fHistodEdxMax[ivar] && cosmic->fHistodEdxMax[ivar]){
310 fHistodEdxMax[ivar]->Add(cosmic->fHistodEdxMax[ivar]);
312 if (fHistodEdxTot[ivar] && cosmic->fHistodEdxTot[ivar]){
313 fHistodEdxTot[ivar]->Add(cosmic->fHistodEdxTot[ivar]);
316 if (cosmic->fCosmicTree){
318 fCosmicTree = new TTree("pairs","pairs");
319 fCosmicTree->SetDirectory(0);
321 AliTPCcalibCosmic::AddTree(fCosmicTree,cosmic->fCosmicTree);
328 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
330 // Process of the ESD event - fill calibration components
333 Printf("ERROR: ESD not available");
339 // COSMIC not signed properly
340 // UInt_t specie = event->GetEventSpecie(); // select only cosmic events
341 //if (specie==AliRecoParam::kCosmic || specie==AliRecoParam::kCalib) {
346 FindCosmicPairs(event);
347 //const AliMultiplicity *multiplicity = event->GetMultiplicity();
348 // Int_t ntracklets = multiplicity->GetNumberOfTracklets();
349 //if (ntracklets>6) return; // filter out "normal" event with high multiplicity
350 //const TString &trigger = event->GetFiredTriggerClasses();
351 //if (trigger.Contains("C0OB0")==0) return;
354 FindPairs(event); // nearly everything takes place in find pairs...
356 if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
357 Int_t ntracks=event->GetNumberOfTracks();
358 fHistNTracks->Fill(ntracks);
363 void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, const AliExternalTrackParam *par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam */*inner1*/, AliTPCseed *seed0, AliTPCseed *seed1, const AliExternalTrackParam *param0Combined , Int_t cross){
365 // par0,par1 - parameter of tracks at DCA 0
366 // inner0,inner1 - parameter of tracks at the TPC entrance
367 // seed0, seed1 - detailed track information
368 // param0Combined - Use combined track parameters for binning
370 Int_t kMinCldEdx =20;
371 Int_t ncl0 = seed0->GetNumberOfClusters();
372 Int_t ncl1 = seed1->GetNumberOfClusters();
373 const Double_t kpullCut = 10;
379 Double_t radius0 = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]);
380 Double_t radius1 = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
381 inner0->GetXYZ(xyz0);
382 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
384 x[1] = TMath::Min(ncl0,ncl1);
385 x[2] = (radius0+radius1)*0.5;
386 x[3] = param0Combined->GetZ();
387 x[4] = inner0->GetSnp();
388 x[5] = param0Combined->GetTgl();
389 x[6] = param0Combined->GetSigned1Pt();
390 x[7] = param0Combined->Pt();
396 delta[0] = (par0->GetY()+par1->GetY());
397 delta[1] = (par0->GetZ()-par1->GetZ());
398 delta[2] = (par0->GetAlpha()-par1->GetAlpha()-TMath::Pi());
399 delta[3] = (par0->GetTgl()+par1->GetTgl());
400 delta[4] = (par0->GetParameter()[4]+par1->GetParameter()[4]);
401 delta[5] = (par0->Pt()-par1->Pt())/((par0->Pt()+par1->Pt())*0.5);
403 sigma[0] = TMath::Sqrt(par0->GetSigmaY2()+par1->GetSigmaY2());
404 sigma[1] = TMath::Sqrt(par0->GetSigmaZ2()+par1->GetSigmaZ2());
405 sigma[2] = TMath::Sqrt(par0->GetSigmaSnp2()+par1->GetSigmaSnp2());
406 sigma[3] = TMath::Sqrt(par0->GetSigmaTgl2()+par1->GetSigmaTgl2());
407 sigma[4] = TMath::Sqrt(par0->GetSigma1Pt2()+par1->GetSigma1Pt2());
408 sigma[5] = sigma[4]*((par0->Pt()+par1->Pt())*0.5);
411 for (Int_t ivar=0;ivar<6;ivar++){
412 if (sigma[ivar]==0) isOK=kFALSE;
413 x[0]= delta[ivar]/sigma[ivar];
414 if (TMath::Abs(x[0])>kpullCut) isOK = kFALSE;
418 if (isOK) for (Int_t ivar=0;ivar<6;ivar++){
419 x[0]= delta[ivar]; // Modifiation 10.10 use not normalized deltas
420 if (ivar==2 || ivar ==3) x[0]*=1000; // angles in radian
421 fHistoDelta[ivar]->Fill(x);
423 x[0]= delta[ivar]/sigma[ivar];
424 fHistoPull[ivar]->Fill(x);
429 // Fill dedx performance
431 for (Int_t ipad=0; ipad<4;ipad++){
437 if (ipad==0) row1=63;
438 if (ipad==1) {row0=63; row1=63+64;}
439 if (ipad==2) {row0=128;}
440 Int_t nclUp = TMath::Nint(seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
441 Int_t nclDown = TMath::Nint(seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
442 Int_t minCl = TMath::Min(nclUp,nclDown);
443 if (minCl<kMinCldEdx) continue;
446 Float_t dEdxTotUp = seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
447 Float_t dEdxTotDown = seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
448 Float_t dEdxMaxUp = seed0->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
449 Float_t dEdxMaxDown = seed1->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
451 if (dEdxTotDown<=0) continue;
452 if (dEdxMaxDown<=0) continue;
453 x[0]=dEdxTotUp/dEdxTotDown;
454 fHistodEdxTot[ipad]->Fill(x);
455 x[0]=dEdxMaxUp/dEdxMaxDown;
456 fHistodEdxMax[ipad]->Fill(x);
463 void AliTPCcalibCosmic::FindPairs(const AliESDEvent *event){
467 // Track0 is choosen in upper TPC part
468 // Track1 is choosen in lower TPC part
470 if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
471 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
472 Int_t ntracks=event->GetNumberOfTracks();
473 TObjArray tpcSeeds(ntracks);
474 if (ntracks==0) return;
475 Double_t vtxx[3]={0,0,0};
476 Double_t svtxx[3]={0.000001,0.000001,100.};
477 AliESDVertex vtx(vtxx,svtxx);
481 for (Int_t i=0;i<ntracks;++i) {
482 AliESDtrack *track = event->GetTrack(i);
483 fClusters->Fill(track->GetTPCNcls());
485 const AliExternalTrackParam * trackIn = track->GetInnerParam();
486 const AliExternalTrackParam * trackOut = track->GetOuterParam();
487 if (!trackIn) continue;
488 if (!trackOut) continue;
489 if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // filter laser
492 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
493 if (!friendTrack) continue;
494 TObject *calibObject;
495 AliTPCseed *seed = 0;
496 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
497 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
499 if (seed) tpcSeeds.AddAt(seed,i);
501 Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
502 if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
503 fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159));
505 if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159));
507 // if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159)>300) {
508 // //TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
509 // //if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
510 // // if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
517 if (ntracks<2) return;
521 for (Int_t i=0;i<ntracks;++i) {
522 AliESDtrack *track0 = event->GetTrack(i);
523 // track0 - choosen upper part
524 if (!track0) continue;
525 if (!track0->GetOuterParam()) continue;
526 if (track0->GetOuterParam()->GetAlpha()<0) continue;
528 track0->GetDirection(dir0);
529 for (Int_t j=0;j<ntracks;++j) {
531 AliESDtrack *track1 = event->GetTrack(j);
533 if (!track1) continue;
534 if (!track1->GetOuterParam()) continue;
535 if (track1->GetOuterParam()->GetAlpha()>0) continue;
538 track1->GetDirection(dir1);
540 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
541 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
542 if (! seed0) continue;
543 if (! seed1) continue;
544 Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159);
545 Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159);
547 Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63);
548 Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63);
550 Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159);
551 Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159);
553 Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
554 Float_t d0 = track0->GetLinearD(0,0);
555 Float_t d1 = track1->GetLinearD(0,0);
557 // conservative cuts - convergence to be guarantied
558 // applying before track propagation
559 if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0
560 if (dir>fCutMinDir) continue; // direction vector product
561 Float_t bz = AliTracker::GetBz();
562 Float_t dvertex0[2]; //distance to 0,0
563 Float_t dvertex1[2]; //distance to 0,0
564 track0->GetDZ(0,0,0,bz,dvertex0);
565 track1->GetDZ(0,0,0,bz,dvertex1);
566 if (TMath::Abs(dvertex0[1])>250) continue;
567 if (TMath::Abs(dvertex1[1])>250) continue;
571 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
572 AliExternalTrackParam param0(*track0);
573 AliExternalTrackParam param1(*track1);
575 // Propagate using Magnetic field and correct fo material budget
579 Double_t maxsnp=0.90;
580 AliTracker::PropagateTrackToBxByBz(¶m0,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign0);
581 AliTracker::PropagateTrackToBxByBz(¶m1,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign1);
583 // Propagate rest to the 0,0 DCA - z should be ignored
585 Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
586 Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
588 param0.GetDZ(0,0,0,bz,dvertex0);
589 param1.GetDZ(0,0,0,bz,dvertex1);
590 if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
592 Double_t xyz0[3];//,pxyz0[3];
593 Double_t xyz1[3];//,pxyz1[3];
596 Bool_t isPair = IsPair(¶m0,¶m1);
598 if (isPair) FillAcordeHist(track0);
599 if (isPair &¶m0.Pt()>1) {
600 const TString &trigger = event->GetFiredTriggerClasses();
601 UInt_t specie = event->GetEventSpecie();
602 printf("COSMIC ?\t%s\t%d\t%f\t%f\n", trigger.Data(),specie, param0.GetZ(), param1.GetZ());
605 // combined track params
607 AliExternalTrackParam *par0U=MakeCombinedTrack(¶m0,¶m1);
608 AliExternalTrackParam *par1U=MakeCombinedTrack(¶m1,¶m0);
612 TTreeSRedirector * cstream = GetDebugStreamer();
613 //printf("My stream=%p\n",(void*)cstream);
614 AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
615 AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
616 AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
617 AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
618 Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
619 Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
620 Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
621 Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
625 Int_t cross =0; // 0 no cross, 2 cross on both sides
626 if (isCrossI) cross+=1;
627 if (isCrossO) cross+=1;
628 FillHistoPerformance(¶m0, ¶m1, ip0, ip1, seed0, seed1,par0U, cross);
630 (*cstream) << "Track0" <<
631 "run="<<fRun<< // run number
632 "event="<<fEvent<< // event number
633 "time="<<fTime<< // time stamp of event
634 "trigger="<<fTrigger<< // trigger
635 "triggerClass="<<&fTriggerClass<< // trigger
636 "mag="<<fMagF<< // magnetic field
637 "dir="<<dir<< // direction
638 "OK="<<isPair<< // will be accepted
639 "b0="<<b0<< // propagate status
640 "b1="<<b1<< // propagate status
641 "crossI="<<isCrossI<< // cross inner
642 "crossO="<<isCrossO<< // cross outer
644 "Orig0.=" << track0 << // original track 0
645 "Orig1.=" << track1 << // original track 1
646 "Tr0.="<<¶m0<< // track propagated to the DCA 0,0
647 "Tr1.="<<¶m1<< // track propagated to the DCA 0,0
648 "Ip0.="<<ip0<< // inner param - upper
649 "Ip1.="<<ip1<< // inner param - lower
650 "Op0.="<<op0<< // outer param - upper
651 "Op1.="<<op1<< // outer param - lower
652 "Up0.="<<par0U<< // combined track 0
653 "Up1.="<<par1U<< // combined track 1
655 "v00="<<dvertex0[0]<< // distance using kalman
656 "v01="<<dvertex0[1]<< //
657 "v10="<<dvertex1[0]<< //
658 "v11="<<dvertex1[1]<< //
659 "d0="<<d0<< // linear distance to 0,0
660 "d1="<<d1<< // linear distance to 0,0
664 "x00="<<xyz0[0]<< // global position close to vertex
668 "x10="<<xyz1[0]<< // global position close to vertex
674 "dir00="<<dir0[0]<< // direction upper
678 "dir10="<<dir1[0]<< // direction lower
683 "Seed0.=" << seed0 << // original seed 0
684 "Seed1.=" << seed1 << // original seed 1
686 "dedx0="<<dedx0<< // dedx0 - all
687 "dedx1="<<dedx1<< // dedx1 - all
689 "dedx0I="<<dedx0I<< // dedx0 - inner ROC
690 "dedx1I="<<dedx1I<< // dedx1 - inner ROC
692 "dedx0O="<<dedx0O<< // dedx0 - outer ROC
693 "dedx1O="<<dedx1O<< // dedx1 - outer ROC
706 void AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
708 // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
709 if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
711 const Double_t acordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
712 const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
715 upperTrack->GetXYZ(r);
717 upperTrack->GetDirection(d);
719 z = r[2] + (d[2]/d[1])*(acordePlane - r[1]);
720 x = r[0] + (d[0]/d[1])*(acordePlane - r[1]);
723 x = r[0] + (d[0]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
724 z = r[2] + (d[2]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
727 x = r[0] + (d[0]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
728 z = r[2] + (d[2]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
731 fModules->Fill(z, x);
737 Long64_t AliTPCcalibCosmic::Merge(TCollection *const li) {
742 TIterator* iter = li->MakeIterator();
743 AliTPCcalibCosmic* cal = 0;
745 while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
746 if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
747 //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
751 fHistNTracks->Add(cal->GetHistNTracks());
752 fClusters->Add(cal-> GetHistClusters());
753 fModules->Add(cal->GetHistAcorde());
754 fHistPt->Add(cal->GetHistPt());
755 fDeDx->Add(cal->GetHistDeDx());
756 fDeDxMIP->Add(cal->GetHistMIP());
764 Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1) const{
768 // 0. Same direction - OPOSITE - cutDir +cutT
769 TCut cutDir("cutDir","dir<-0.99")
771 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
774 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
778 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
781 const Double_t *p0 = tr0->GetParameter();
782 const Double_t *p1 = tr1->GetParameter();
783 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
784 if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
785 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
787 Double_t d0[3], d1[3];
788 tr0->GetDirection(d0);
789 tr1->GetDirection(d1);
790 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
797 Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
799 // Calculate the MIP value - gaussian fit used
802 TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
803 funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
804 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
805 hist->Fit(funcDoubleGaus);
806 Double_t aMIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
808 delete funcDoubleGaus;
817 void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) {
819 // Not implemented yet
826 void AliTPCcalibCosmic::BinLogX(THnSparse *const h, Int_t axisDim) {
828 // Method for the correct logarithmic binning of histograms
830 TAxis *axis = h->GetAxis(axisDim);
831 int bins = axis->GetNbins();
833 Double_t from = axis->GetXmin();
834 Double_t to = axis->GetXmax();
835 Double_t *newBins = new Double_t[bins + 1];
838 Double_t factor = pow(to/from, 1./bins);
840 for (int i = 1; i <= bins; i++) {
841 newBins[i] = factor * newBins[i-1];
843 axis->Set(bins, newBins);
849 void AliTPCcalibCosmic::BinLogX(TH1 *const h) {
851 // Method for the correct logarithmic binning of histograms
853 TAxis *axis = h->GetXaxis();
854 int bins = axis->GetNbins();
856 Double_t from = axis->GetXmin();
857 Double_t to = axis->GetXmax();
858 Double_t *newBins = new Double_t[bins + 1];
861 Double_t factor = pow(to/from, 1./bins);
863 for (int i = 1; i <= bins; i++) {
864 newBins[i] = factor * newBins[i-1];
866 axis->Set(bins, newBins);
872 AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
874 // Make a atrack using the kalman update of track0 and track1
876 AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
877 par1R->Rotate(track0->GetAlpha());
878 par1R->PropagateTo(track0->GetX(),AliTracker::GetBz());
881 Double_t * param = (Double_t*)par1R->GetParameter();
882 Double_t * covar = (Double_t*)par1R->GetCovariance();
890 covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
891 //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
896 AliExternalTrackParam *AliTPCcalibCosmic::MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
898 // Make combined track
901 AliExternalTrackParam * par1T = MakeTrack(track0,track1);
902 AliExternalTrackParam * par0U = new AliExternalTrackParam(*track0);
904 UpdateTrack(*par0U,*par1T);
910 void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
912 // Update track 1 with track 2
916 TMatrixD vecXk(5,1); // X vector
917 TMatrixD covXk(5,5); // X covariance
918 TMatrixD matHk(5,5); // vector to mesurement
919 TMatrixD measR(5,5); // measurement error
920 TMatrixD vecZk(5,1); // measurement
922 TMatrixD vecYk(5,1); // Innovation or measurement residual
923 TMatrixD matHkT(5,5);
924 TMatrixD matSk(5,5); // Innovation (or residual) covariance
925 TMatrixD matKk(5,5); // Optimal Kalman gain
926 TMatrixD mat1(5,5); // update covariance matrix
927 TMatrixD covXk2(5,5); //
928 TMatrixD covOut(5,5);
930 Double_t *param1=(Double_t*) track1.GetParameter();
931 Double_t *covar1=(Double_t*) track1.GetCovariance();
932 Double_t *param2=(Double_t*) track2.GetParameter();
933 Double_t *covar2=(Double_t*) track2.GetCovariance();
935 // copy data to the matrix
936 for (Int_t ipar=0; ipar<5; ipar++){
937 for (Int_t jpar=0; jpar<5; jpar++){
938 covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
939 measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
943 vecXk(ipar,0) = param1[ipar];
944 vecZk(ipar,0) = param2[ipar];
953 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
954 matHkT=matHk.T(); matHk.T();
955 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
957 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
958 vecXk += matKk*vecYk; // updated vector
959 covXk2 = (mat1-(matKk*matHk));
960 covOut = covXk2*covXk;
964 // copy from matrix to parameters
977 for (Int_t ipar=0; ipar<5; ipar++){
978 param1[ipar]= vecXk(ipar,0) ;
979 for (Int_t jpar=0; jpar<5; jpar++){
980 covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
987 void AliTPCcalibCosmic::FindCosmicPairs(const AliESDEvent * event) {
989 // find cosmic pairs trigger by random trigger
992 AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD();
993 AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC();
994 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
995 const Double_t kMinPt=1;
996 const Double_t kMinPtMax=0.8;
997 const Double_t kMinNcl=50;
998 const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};
999 Int_t ntracks=event->GetNumberOfTracks();
1000 // Float_t dcaTPC[2]={0,0};
1001 // Float_t covTPC[3]={0,0,0};
1003 UInt_t specie = event->GetEventSpecie(); // skip laser events
1004 if (specie==AliRecoParam::kCalib) return;
1008 for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
1009 AliESDtrack *track0 = event->GetTrack(itrack0);
1010 if (!track0) continue;
1011 if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;
1013 if (TMath::Abs(AliTracker::GetBz())>1&&track0->Pt()<kMinPt) continue;
1014 if (track0->GetTPCncls()<kMinNcl) continue;
1015 if (TMath::Abs(track0->GetY())<kMaxDelta[0]) continue;
1016 if (track0->GetKinkIndex(0)>0) continue;
1017 const Double_t * par0=track0->GetParameter(); //track param at rhe DCA
1020 //track0->GetImpactParametersTPC(dcaTPC,covTPC);
1021 //if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
1022 //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
1023 // const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();
1024 for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
1025 AliESDtrack *track1 = event->GetTrack(itrack1);
1026 if (!track1) continue;
1027 if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;
1028 if (track1->GetKinkIndex(0)>0) continue;
1029 if (TMath::Abs(AliTracker::GetBz())>1&&track1->Pt()<kMinPt) continue;
1030 if (track1->GetTPCncls()<kMinNcl) continue;
1031 if (TMath::Abs(AliTracker::GetBz())>1&&TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
1032 if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;
1033 //track1->GetImpactParametersTPC(dcaTPC,covTPC);
1034 // if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
1035 //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
1037 const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
1039 Bool_t isPair=kTRUE;
1040 for (Int_t ipar=0; ipar<5; ipar++){
1041 if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
1042 if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
1044 if (!isPair) continue;
1045 if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
1046 //delta with correct sign
1048 TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"
1049 TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"
1050 TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"
1052 if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign
1053 if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
1054 if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
1055 if (!isPair) continue;
1056 TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1057 Int_t eventNumber = event->GetEventNumberInFile();
1058 Bool_t hasFriend=(esdFriend) ? (esdFriend->GetTrack(itrack0)!=0):0;
1059 Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4);
1060 printf("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS);
1063 // const AliExternalTrackParam * trackIn1 = track1->GetInnerParam();
1066 TTreeSRedirector * pcstream = GetDebugStreamer();
1067 Int_t ntracksSPD = vertexSPD->GetNContributors();
1068 Int_t ntracksTPC = vertexTPC->GetNContributors();
1070 AliESDfriendTrack *friendTrack0 = esdFriend->GetTrack(itrack0);
1071 if (!friendTrack0) continue;
1072 AliESDfriendTrack *friendTrack1 = esdFriend->GetTrack(itrack1);
1073 if (!friendTrack1) continue;
1074 TObject *calibObject;
1075 AliTPCseed *seed0 = 0;
1076 AliTPCseed *seed1 = 0;
1078 for (Int_t l=0;(calibObject=friendTrack0->GetCalibObject(l));++l) {
1079 if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
1081 for (Int_t l=0;(calibObject=friendTrack1->GetCalibObject(l));++l) {
1082 if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
1086 (*pcstream)<<"pairs"<<
1087 "run="<<fRun<< // run number
1088 "event="<<fEvent<< // event number
1089 "time="<<fTime<< // time stamp of event
1090 "trigger="<<fTrigger<< // trigger
1091 "triggerClass="<<&fTriggerClass<< // trigger
1092 "mag="<<fMagF<< // magnetic field
1094 "nSPD="<<ntracksSPD<<
1095 "nTPC="<<ntracksTPC<<
1096 "vSPD.="<<vertexSPD<< //primary vertex -SPD
1097 "vTPC.="<<vertexTPC<< //primary vertex -TPC
1098 "t0.="<<track0<< //track0
1099 "t1.="<<track1<< //track1
1100 "ft0.="<<friendTrack0<< //track0
1101 "ft1.="<<friendTrack1<< //track1
1102 "s0.="<<seed0<< //track0
1103 "s1.="<<seed1<< //track1
1107 fCosmicTree = new TTree("pairs","pairs");
1108 fCosmicTree->SetDirectory(0);
1110 if (fCosmicTree->GetEntries()==0){
1112 fCosmicTree->SetDirectory(0);
1113 fCosmicTree->Branch("t0.",&track0);
1114 fCosmicTree->Branch("t1.",&track1);
1115 fCosmicTree->Branch("ft0.",&friendTrack0);
1116 fCosmicTree->Branch("ft1.",&friendTrack1);
1118 fCosmicTree->SetBranchAddress("t0.",&track0);
1119 fCosmicTree->SetBranchAddress("t1.",&track1);
1120 fCosmicTree->SetBranchAddress("ft0.",&friendTrack0);
1121 fCosmicTree->SetBranchAddress("ft1.",&friendTrack1);
1123 fCosmicTree->Fill();
1129 void AliTPCcalibCosmic::Terminate(){
1131 // copy the cosmic tree to memory resident tree
1133 static Int_t counter=0;
1134 printf("AliTPCcalibCosmic::Terminate\t%d\n",counter);
1136 AliTPCcalibBase::Terminate();
1140 void AliTPCcalibCosmic::AddTree(TTree* treeOutput, TTree * treeInput){
1142 // Add the content of tree:
1143 // Notice automatic copy of tree in ROOT does not work for such complicated tree
1146 //if (TMath::Abs(fMagF)<0.1) return; // work around - otherwise crashes
1147 AliESDtrack *track0=new AliESDtrack;
1148 AliESDtrack *track1=new AliESDtrack;
1149 AliESDfriendTrack *ftrack0=new AliESDfriendTrack;
1150 AliESDfriendTrack *ftrack1=new AliESDfriendTrack;
1151 treeInput->SetBranchAddress("t0.",&track0);
1152 treeInput->SetBranchAddress("t1.",&track1);
1153 treeInput->SetBranchAddress("ft0.",&ftrack0);
1154 treeInput->SetBranchAddress("ft1.",&ftrack1);
1155 treeOutput->SetDirectory(0);
1157 Int_t entries= treeInput->GetEntries();
1158 Int_t step=1+Int_t(TMath::Log(1+treeOutput->GetEntries()/10.));
1159 for (Int_t i=0; i<entries; i+=step){
1160 treeInput->SetBranchAddress("t0.",&track0);
1161 treeInput->SetBranchAddress("t1.",&track1);
1162 treeInput->SetBranchAddress("ft0.",&ftrack0);
1163 treeInput->SetBranchAddress("ft1.",&ftrack1);
1164 treeInput->GetEntry(i);
1165 if (!track0) continue;
1166 if (!track1) continue;
1167 if (!ftrack0) continue;
1168 if (!ftrack1) continue;
1169 if (track0->GetTPCncls()<=0) continue;
1170 if (track1->GetTPCncls()<=0) continue;
1171 if (!track0->GetInnerParam()) continue;
1172 if (!track1->GetInnerParam()) continue;
1173 if (!track0->GetTPCInnerParam()) continue;
1174 if (!track1->GetTPCInnerParam()) continue;
1176 treeOutput->SetBranchAddress("t0.",&track0);
1177 treeOutput->SetBranchAddress("t1.",&track1);
1178 treeOutput->SetBranchAddress("ft0.",&ftrack0);
1179 treeOutput->SetBranchAddress("ft1.",&ftrack1);
1194 void AliTPCcalibCosmic::MakeFitTree(TTree * treeInput, TTreeSRedirector *pcstream, const TObjArray * corrArray, Int_t step, Int_t run){
1197 // refit the tracks with original points + corrected points for each correction
1199 // treeInput - tree with cosmic tracks
1200 // pcstream - debug output
1203 // Loop over pair of cosmic tracks:
1204 // 1. Find trigger offset between cosmic event and "physic" trigger
1205 // 2. Refit tracks with current transformation
1206 // 3. Refit tracks using additional "primitive" distortion on top
1207 // Best correction estimated as linear combination of corrections
1208 // minimizing the observed distortions
1209 // Observed distortions - matching in the y,z, snp, theta and 1/pt
1211 const Double_t kResetCov=20.;
1212 const Double_t kMaxDelta[5]={1,40,0.03,0.01,0.2};
1213 const Double_t kSigma=2.;
1214 const Double_t kMaxTime=1050;
1215 const Double_t kMaxSnp=0.8;
1216 Int_t ncorr=corrArray->GetEntries();
1217 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1218 AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
1219 AliGRPObject* grp = AliTPCcalibDB::Instance()->GetGRP(run);
1220 Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
1221 transform->SetCurrentRun(run);
1222 transform->SetCurrentTimeStamp(TMath::Nint(time));
1224 for (Int_t i=0;i<15;i++) covar[i]=0;
1225 covar[0]=kSigma*kSigma;
1226 covar[2]=kSigma*kSigma;
1227 covar[5]=kSigma*kSigma/Float_t(150*150);
1228 covar[9]=kSigma*kSigma/Float_t(150*150);
1230 Double_t *distortions = new Double_t[ncorr+1];
1232 AliESDtrack *track0=new AliESDtrack;
1233 AliESDtrack *track1=new AliESDtrack;
1234 AliESDfriendTrack *ftrack0=new AliESDfriendTrack;
1235 AliESDfriendTrack *ftrack1=new AliESDfriendTrack;
1236 treeInput->SetBranchAddress("t0.",&track0);
1237 treeInput->SetBranchAddress("t1.",&track1);
1238 treeInput->SetBranchAddress("ft0.",&ftrack0);
1239 treeInput->SetBranchAddress("ft1.",&ftrack1);
1240 Int_t entries= treeInput->GetEntries();
1241 for (Int_t i=0; i<entries; i+=step){
1242 treeInput->GetEntry(i);
1243 if (i%20==0) printf("%d\n",i);
1244 if (!ftrack0->GetTPCOut()) continue;
1245 if (!ftrack1->GetTPCOut()) continue;
1246 AliTPCseed *seed0=0;
1247 AliTPCseed *seed1=0;
1248 TObject *calibObject;
1249 for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
1250 if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
1252 for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
1253 if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
1255 if (!seed0) continue;
1256 if (!seed1) continue;
1257 if (TMath::Abs(seed0->GetSnp())>kMaxSnp) continue;
1258 if (TMath::Abs(seed1->GetSnp())>kMaxSnp) continue;
1261 Int_t nclA0=0, nclC0=0; // number of clusters
1262 Int_t nclA1=0, nclC1=0; // number of clusters
1263 Int_t ncl0=0,ncl1=0;
1264 Double_t rmin0=300, rmax0=-300; // variables to estimate the time0 - trigger offset
1265 Double_t rmin1=300, rmax1=-300;
1266 Double_t tmin0=2000, tmax0=-2000;
1267 Double_t tmin1=2000, tmax1=-2000;
1270 // calculate trigger offset usig "missing clusters"
1271 for (Int_t irow=0; irow<159; irow++){
1272 AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
1273 if (cluster0 &&cluster0->GetX()>10){
1274 if (cluster0->GetX()<rmin0) { rmin0=cluster0->GetX(); tmin0=cluster0->GetTimeBin();}
1275 if (cluster0->GetX()>rmax0) { rmax0=cluster0->GetX(); tmax0=cluster0->GetTimeBin();}
1277 if (cluster0->GetDetector()%36<18) nclA0++;
1278 if (cluster0->GetDetector()%36>=18) nclC0++;
1280 AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
1281 if (cluster1&&cluster1->GetX()>10){
1282 if (cluster1->GetX()<rmin1) { rmin1=cluster1->GetX(); tmin1=cluster1->GetTimeBin();}
1283 if (cluster1->GetX()>rmax1) { rmax1=cluster1->GetX(); tmax1=cluster1->GetTimeBin();}
1285 if (cluster1->GetDetector()%36<18) nclA1++;
1286 if (cluster1->GetDetector()%36>=18) nclC1++;
1289 Int_t cosmicType=0; // types of cosmic topology
1290 if ((nclA0>nclC0) && (nclA1>nclC1)) cosmicType=0; // AA side
1291 if ((nclA0<nclC0) && (nclA1<nclC1)) cosmicType=1; // CC side
1292 if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=2; // AC side
1293 if ((nclA0<nclC0) && (nclA1>nclC1)) cosmicType=3; // CA side
1294 //if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=6; // AC side out of time
1295 //if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=7; // CA side out of time
1297 Double_t deltaTime = 0; // correction for trigger not in time - equalizing the time arival
1298 if ((cosmicType>1)&&TMath::Abs(track1->GetZ()-track0->GetZ())>4){
1300 deltaTime=0.5*(track1->GetZ()-track0->GetZ())/param->GetZWidth();
1301 if (nclA0>nclC0) deltaTime*=-1; // if A side track
1304 TVectorD vectorDT(8);
1305 Int_t crossCounter=0;
1306 Double_t deltaTimeCross = AliTPCcalibCosmic::GetDeltaTime(rmin0, rmax0, rmin1, rmax1, tmin0, tmax0, tmin1, tmax1, TMath::Abs(track0->GetY()),vectorDT);
1307 Bool_t isOKTrigger=kTRUE;
1308 for (Int_t ic=0; ic<6;ic++) {
1309 if (TMath::Abs(vectorDT[ic])>0) {
1310 if (vectorDT[ic]+vectorDT[6]<0) isOKTrigger=kFALSE;
1311 if (vectorDT[ic]+vectorDT[7]>kMaxTime) isOKTrigger=kFALSE;
1317 Double_t deltaTimeCluster=deltaTime;
1318 if ((cosmicType==0 || cosmicType==1) && crossCounter>0){
1319 deltaTimeCluster=deltaTimeCross;
1322 if (nclA0*nclC0>0 || nclA1*nclC1>0) cosmicType+=16; // mixed A side C side - bad for visualization
1324 // Apply current transformation
1327 for (Int_t irow=0; irow<159; irow++){
1328 AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
1329 if (cluster0 &&cluster0->GetX()>10){
1330 Double_t x0[3]={cluster0->GetRow(),cluster0->GetPad(),cluster0->GetTimeBin()+deltaTimeCluster};
1331 Int_t index0[1]={cluster0->GetDetector()};
1332 transform->Transform(x0,index0,0,1);
1333 cluster0->SetX(x0[0]);
1334 cluster0->SetY(x0[1]);
1335 cluster0->SetZ(x0[2]);
1338 AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
1339 if (cluster1&&cluster1->GetX()>10){
1340 Double_t x1[3]={cluster1->GetRow(),cluster1->GetPad(),cluster1->GetTimeBin()+deltaTimeCluster};
1341 Int_t index1[1]={cluster1->GetDetector()};
1342 transform->Transform(x1,index1,0,1);
1343 cluster1->SetX(x1[0]);
1344 cluster1->SetY(x1[1]);
1345 cluster1->SetZ(x1[2]);
1350 Double_t alpha=track0->GetAlpha(); // rotation frame
1351 Double_t cos = TMath::Cos(alpha);
1352 Double_t sin = TMath::Sin(alpha);
1353 Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
1354 AliExternalTrackParam btrack0=*(ftrack0->GetTPCOut());
1355 AliExternalTrackParam btrack1=*(ftrack1->GetTPCOut());
1356 btrack0.Rotate(alpha);
1357 btrack1.Rotate(alpha);
1358 // change the sign for track 1
1359 Double_t* par1=(Double_t*)btrack0.GetParameter();
1362 btrack0.AddCovariance(covar);
1363 btrack1.AddCovariance(covar);
1364 btrack0.ResetCovariance(kResetCov);
1365 btrack1.ResetCovariance(kResetCov);
1368 TObjArray tracks0(ncorr+1);
1369 TObjArray tracks1(ncorr+1);
1371 Double_t dEdx0Tot=seed0->CookdEdxAnalytical(0.02,0.6,kTRUE);
1372 Double_t dEdx1Tot=seed1->CookdEdxAnalytical(0.02,0.6,kTRUE);
1373 Double_t dEdx0Max=seed0->CookdEdxAnalytical(0.02,0.6,kFALSE);
1374 Double_t dEdx1Max=seed1->CookdEdxAnalytical(0.02,0.6,kFALSE);
1375 //if (TMath::Abs((dEdx0Max+1)/(dEdx0Tot+1)-1.)>0.1) isOK=kFALSE;
1376 //if (TMath::Abs((dEdx1Max+1)/(dEdx1Tot+1)-1.)>0.1) isOK=kFALSE;
1378 for (Int_t icorr=-1; icorr<ncorr; icorr++){
1379 AliExternalTrackParam rtrack0=btrack0;
1380 AliExternalTrackParam rtrack1=btrack1;
1381 AliTPCCorrection *corr = 0;
1382 if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
1384 for (Int_t irow=159; irow>0; irow--){
1385 AliTPCclusterMI *cluster=seed0->GetClusterPointer(irow);
1386 if (!cluster) continue;
1388 Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1389 transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD); // transform to global
1390 Float_t r[3]={rD[0],rD[1],rD[2]};
1392 corr->DistortPoint(r, cluster->GetDetector());
1395 Double_t cov[3]={0.01,0.,0.01};
1396 Double_t lx =cos*r[0]+sin*r[1];
1397 Double_t ly =-sin*r[0]+cos*r[1];
1398 rD[1]=ly; rD[0]=lx; rD[2]=r[2]; //transform to track local
1399 if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, lx,mass,1.,kFALSE)) isOKT=kFALSE;;
1400 if (!rtrack0.Update(&rD[1],cov)) isOKT =kFALSE;
1401 if (icorr<0) ncl0++;
1404 for (Int_t irow=159; irow>0; irow--){
1405 AliTPCclusterMI *cluster=seed1->GetClusterPointer(irow);
1406 if (!cluster) continue;
1408 Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1409 transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD);
1410 Float_t r[3]={rD[0],rD[1],rD[2]};
1412 corr->DistortPoint(r, cluster->GetDetector());
1415 Double_t cov[3]={0.01,0.,0.01};
1416 Double_t lx =cos*r[0]+sin*r[1];
1417 Double_t ly =-sin*r[0]+cos*r[1];
1418 rD[1]=ly; rD[0]=lx; rD[2]=r[2];
1419 if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, lx,mass,1.,kFALSE)) isOKT=kFALSE;
1420 if (!rtrack1.Update(&rD[1],cov)) isOKT=kFALSE;
1421 if (icorr<0) ncl1++;
1423 if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,10.,kFALSE)) isOKT=kFALSE;
1424 if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,10.,kFALSE)) isOKT=kFALSE;
1425 if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,1.,kFALSE)) isOKT=kFALSE;
1426 if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,1.,kFALSE)) isOKT=kFALSE;
1427 const Double_t *param0=rtrack0.GetParameter();
1428 const Double_t *param1=rtrack1.GetParameter();
1429 for (Int_t ipar=0; ipar<4; ipar++){
1430 if (TMath::Abs(param1[ipar]-param0[ipar])>kMaxDelta[ipar]) isOK=kFALSE;
1432 tracks0.AddAt(rtrack0.Clone(), icorr+1);
1433 tracks1.AddAt(rtrack1.Clone(), icorr+1);
1434 AliExternalTrackParam out0=*(ftrack0->GetTPCOut());
1435 AliExternalTrackParam out1=*(ftrack1->GetTPCOut());
1436 Int_t nentries=TMath::Min(ncl0,ncl1);
1439 (*pcstream)<<"cosmic"<<
1440 "isOK="<<isOK<< // correct all propagation update and also residuals
1441 "isOKT="<<isOKT<< // correct all propagation update
1442 "isOKTrigger="<<isOKTrigger<< // correct? estimate of trigger offset
1446 "cross="<<crossCounter<<
1447 "vDT.="<<&vectorDT<<
1449 "dTime="<<deltaTime<< // delta time using the A-c side cross
1450 "dTimeCross="<<deltaTimeCross<< // delta time using missing clusters
1452 "dEdx0Max="<<dEdx0Max<<
1453 "dEdx0Tot="<<dEdx0Tot<<
1454 "dEdx1Max="<<dEdx1Max<<
1455 "dEdx1Tot="<<dEdx1Tot<<
1457 "track0.="<<track0<< // original track 0
1458 "track1.="<<track1<< // original track 1
1459 "out0.="<<&out0<< // outer track 0
1460 "out1.="<<&out1<< // outer track 1
1461 "rtrack0.="<<&rtrack0<< // refitted track with current transform
1462 "rtrack1.="<<&rtrack1<< //
1465 "entries="<<nentries<< // number of clusters
1472 Int_t nentries=TMath::Min(ncl0,ncl1);
1473 for (Int_t ipar=0; ipar<5; ipar++){
1474 for (Int_t icorr=-1; icorr<ncorr; icorr++){
1475 AliTPCCorrection *corr = 0;
1476 if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
1478 AliExternalTrackParam *param0=(AliExternalTrackParam *) tracks0.At(icorr+1);
1479 AliExternalTrackParam *param1=(AliExternalTrackParam *) tracks1.At(icorr+1);
1480 distortions[icorr+1]=param1->GetParameter()[ipar]-param0->GetParameter()[ipar];
1482 distortions[icorr+1]-=distortions[0];
1486 Double_t bz=AliTrackerBase::GetBz();
1488 param0->GetXYZ(gxyz);
1490 Double_t theta=param0->GetParameter()[3];
1491 Double_t phi = alpha;
1492 Double_t snp = track0->GetInnerParam()->GetSnp();
1493 Double_t mean= distortions[0];
1494 Int_t index = param0->GetIndex(ipar,ipar);
1495 Double_t rms=TMath::Sqrt(param1->GetCovariance()[index]+param1->GetCovariance()[index]);
1496 if (crossCounter<1) rms*=1;
1497 Double_t sector=9*phi/TMath::Pi();
1498 Double_t dsec = sector-TMath::Nint(sector+0.5);
1499 Double_t gx=gxyz[0],gy=gxyz[1],gz=gxyz[2];
1500 Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
1502 // Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5;
1503 Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5;
1505 (*pcstream)<<"fit"<< // dump valus for fit
1506 "run="<<run<< //run number
1507 "bz="<<bz<< // magnetic filed used
1508 "dtype="<<dtype<< // detector match type 20
1509 "ptype="<<ipar<< // parameter type
1510 "theta="<<theta<< // theta
1511 "phi="<<phi<< // phi
1512 "snp="<<snp<< // snp
1513 "mean="<<mean<< // mean dist value
1514 "rms="<<rms<< // rms
1518 "refX="<<refX<< // reference radius
1519 "gx="<<gx<< // global position
1520 "gy="<<gy<< // global position
1521 "gz="<<gz<< // global position
1522 "dRrec="<<dRrec<< // delta Radius in reconstruction
1524 "id="<<cosmicType<< //type of the cosmic used
1525 "entries="<<nentries;// number of entries in bin
1526 (*pcstream)<<"cosmicDebug"<<
1527 "p0.="<<param0<< // dump distorted track 0
1528 "p1.="<<param1; // dump distorted track 1
1531 (*pcstream)<<"fit"<<
1532 Form("%s=",corr->GetName())<<distortions[icorr+1]; // dump correction value
1533 (*pcstream)<<"cosmicDebug"<<
1534 Form("%s=",corr->GetName())<<distortions[icorr+1]<< // dump correction value
1535 Form("%sp0.=",corr->GetName())<<param0<< // dump distorted track 0
1536 Form("%sp1.=",corr->GetName())<<param1; // dump distorted track 1
1538 } //loop corrections
1539 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1540 (*pcstream)<<"cosmicDebug"<<"isOK="<<isOK<<"\n";
1541 } //loop over parameters
1544 delete [] distortions;
1549 Double_t AliTPCcalibCosmic::GetDeltaTime(Double_t rmin0, Double_t rmax0, Double_t rmin1, Double_t rmax1, Double_t tmin0, Double_t tmax0, Double_t tmin1, Double_t tmax1, Double_t dcaR, TVectorD &vectorDT)
1552 // Estimate trigger offset between random cosmic event and "physics" trigger
1553 // Efficiency about 50 % of cases:
1555 // 0. Tracks crossing A side and C side - no match in z - 30 % of cases
1556 // 1. Track only on one side and dissapear at small or at high radiuses - 50 % of cases
1557 // 1.a) rmax<Rc && tmax<Tcmax && tmax>tmin => deltaT=Tcmax-tmax
1558 // 1.b) rmin>Rcmin && tmin<Tcmax && tmin>tmax => deltaT=Tcmax-tmin
1559 // // case the z matching gives proper time
1560 // 1.c) rmax<Rc && tmax>Tcmin && tmax<tmin => deltaT=-tmax
1563 // TCut cutStop = "(min(rmax0,rmax1)<235||abs(rmin0-rmin1)>10)"; // tracks not registered for full time
1568 // 0-3 - occur - wrong correlation
1569 // 1-2 - occur - wrong correlation
1571 // 2-3 - occur - small number of outlyers -20%
1578 const Double_t kMaxRCut=235; // max radius
1579 const Double_t kMinRCut=TMath::Max(dcaR,90.); // min radius
1580 const Double_t kMaxDCut=30; // max distance for minimal radius
1581 const Double_t kMinTime=110;
1582 const Double_t kMaxTime=950;
1585 vectorDT[6]=TMath::Min(TMath::Min(tmin0,tmin1),TMath::Min(tmax0,tmax1));
1586 vectorDT[7]=TMath::Max(TMath::Max(tmin0,tmin1),TMath::Max(tmax0,tmax1));
1587 if (TMath::Min(rmax0,rmax1)<kMaxRCut){
1588 // max cross - deltaT>0
1589 if (rmax0<kMaxRCut && tmax0 <kMaxTime && tmax0>tmin0) vectorDT[0]=kMaxTime-tmax0; // disapear at CE
1590 if (rmax1<kMaxRCut && tmax1 <kMaxTime && tmax1>tmin1) vectorDT[1]=kMaxTime-tmax1; // disapear at CE
1591 // min cross - deltaT<0 - OK they are correlated
1592 if (rmax0<kMaxRCut && tmax0 >kMinTime && tmax0<tmin0) vectorDT[2]=-tmax0; // disapear at ROC
1593 if (rmax1<kMaxRCut && tmax1 >kMinTime && tmax1<tmin1) vectorDT[3]=-tmax1; // disapear at ROC
1595 if (rmin0> kMinRCut+kMaxDCut && tmin0 <kMaxTime && tmin0>tmax0) vectorDT[4]=kMaxTime-tmin0;
1596 if (rmin1> kMinRCut+kMaxDCut && tmin1 <kMaxTime && tmin1>tmax1) vectorDT[5]=kMaxTime-tmin1;
1598 for (Int_t i=0; i<6;i++) {
1599 if (TMath::Abs(vectorDT[i])>0) {
1601 if (vectorDT[i]+vectorDT[6]<0) isOK=kFALSE;
1602 if (vectorDT[i]+vectorDT[7]>kMaxTime) isOK=kFALSE;
1603 if (isOK) deltaT=vectorDT[i];