Line data Source code
1 :
2 :
3 : /**************************************************************************
4 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 : * *
6 : * Author: The ALICE Off-line Project. *
7 : * Contributors are mentioned in the code where appropriate. *
8 : * *
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 : **************************************************************************/
17 :
18 : /*
19 : Comments to be written here:
20 : 1. What do we calibrate.
21 : 2. How to interpret results
22 : 3. Simple example
23 : 4. Analysis using debug streamers.
24 :
25 :
26 :
27 : 3.Simple example
28 : // To make cosmic scan the user interaction neccessary
29 : //
30 :
31 : */
32 :
33 :
34 :
35 : #include "Riostream.h"
36 : #include "TChain.h"
37 : #include "TTree.h"
38 : #include "TH1F.h"
39 : #include "TH2F.h"
40 : #include "TList.h"
41 : #include "TMath.h"
42 : #include "TCanvas.h"
43 : #include "TFile.h"
44 : #include "TF1.h"
45 : #include "THnSparse.h"
46 : #include "TDatabasePDG.h"
47 :
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"
55 :
56 : #include "AliTracker.h"
57 : #include "AliMagF.h"
58 : #include "AliTPCCalROC.h"
59 : #include "AliTPCParam.h"
60 : #include "AliLog.h"
61 :
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 : #include "AliTPCreco.h"
74 6 : ClassImp(AliTPCcalibCosmic)
75 :
76 :
77 : AliTPCcalibCosmic::AliTPCcalibCosmic()
78 0 : :AliTPCcalibBase(),
79 0 : fHistNTracks(0),
80 0 : fClusters(0),
81 0 : fModules(0),
82 0 : fHistPt(0),
83 0 : fDeDx(0),
84 0 : fDeDxMIP(0),
85 0 : fMIPvalue(1),
86 0 : fCutMaxD(5), // maximal distance in rfi ditection
87 0 : fCutMaxDz(40), // maximal distance in z ditection
88 0 : fCutTheta(0.03), // maximal distan theta
89 0 : fCutMinDir(-0.99), // direction vector products
90 0 : fCosmicTree(0) // tree with cosmic data
91 0 : {
92 : //
93 : // CONSTRUCTOR - SEE COMMENTS ABOVE
94 : //
95 0 : AliInfo("Default Constructor");
96 0 : for (Int_t ihis=0; ihis<6;ihis++){
97 0 : fHistoDelta[ihis]=0;
98 0 : fHistoPull[ihis]=0;
99 : }
100 0 : for (Int_t ihis=0; ihis<4;ihis++){
101 0 : fHistodEdxMax[ihis] =0;
102 0 : fHistodEdxTot[ihis] =0;
103 : }
104 0 : }
105 :
106 :
107 : AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
108 0 : :AliTPCcalibBase(),
109 0 : fHistNTracks(0),
110 0 : fClusters(0),
111 0 : fModules(0),
112 0 : fHistPt(0),
113 0 : fDeDx(0),
114 0 : fDeDxMIP(0),
115 0 : fMIPvalue(1),
116 0 : fCutMaxD(5), // maximal distance in rfi ditection
117 0 : fCutMaxDz(40), // maximal distance in z ditection
118 0 : fCutTheta(0.03), // maximal distan theta
119 0 : fCutMinDir(-0.99), // direction vector products
120 0 : fCosmicTree(0) // tree with cosmic data
121 0 : {
122 : //
123 : // cONSTRUCTOR - SEE COMENTS ABOVE
124 : //
125 0 : SetName(name);
126 0 : SetTitle(title);
127 :
128 0 : fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
129 0 : fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
130 0 : fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
131 0 : fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
132 0 : fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
133 0 : BinLogX(fDeDx);
134 0 : fDeDxMIP = new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
135 0 : Init();
136 0 : AliInfo("Non Default Constructor");
137 : //
138 0 : }
139 :
140 0 : AliTPCcalibCosmic::~AliTPCcalibCosmic(){
141 : //
142 : // destructor
143 : //
144 0 : for (Int_t ihis=0; ihis<6;ihis++){
145 0 : delete fHistoDelta[ihis];
146 0 : delete fHistoPull[ihis];
147 : }
148 0 : for (Int_t ihis=0; ihis<4;ihis++){
149 0 : delete fHistodEdxTot[ihis];
150 0 : delete fHistodEdxMax[ihis];
151 : }
152 :
153 0 : delete fHistNTracks; // histogram showing number of ESD tracks per event
154 0 : delete fClusters; // histogram showing the number of clusters per track
155 0 : delete fModules; // 2d histogram of tracks which are propagated to the ACORDE scintillator array
156 0 : delete fHistPt; // Pt histogram of reconstructed tracks
157 0 : delete fDeDx; // dEdx spectrum showing the different particle types
158 0 : delete fDeDxMIP; // TPC signal close to the MIP region of muons 0.4 < p < 0.45 GeV
159 0 : }
160 :
161 :
162 : void AliTPCcalibCosmic::Init(){
163 : //
164 : // init component
165 : // Make performance histograms
166 : //
167 :
168 : // tracking performance bins
169 : // 0 - delta of interest
170 : // 1 - min (track0, track1) number of clusters
171 : // 2 - R - vertex radius
172 : // 3 - P1 - mean z
173 : // 4 - P2 - snp(phi) at inner wall of TPC
174 : // 5 - P3 - tan(theta) at inner wall of TPC
175 : // 6 - P4 - 1/pt mean
176 : // 7 - pt - pt mean
177 : // 8 - alpha
178 : // 9 - is corss indicator
179 : Int_t ndim=10;
180 0 : Double_t xminTrack[10], xmaxTrack[10];
181 0 : Int_t binsTrack[10];
182 0 : TString axisName[10];
183 : //
184 0 : binsTrack[0] =100;
185 0 : axisName[0] ="#Delta";
186 : //
187 0 : binsTrack[1] =8;
188 0 : xminTrack[1] =80; xmaxTrack[1]=160;
189 0 : axisName[1] ="N_{cl}";
190 : //
191 0 : binsTrack[2] =10;
192 0 : xminTrack[2] =0; xmaxTrack[2]=90; //
193 0 : axisName[2] ="dca_{r} (cm)";
194 : //
195 0 : binsTrack[3] =25;
196 0 : xminTrack[3] =-250; xmaxTrack[3]=250; //
197 0 : axisName[3] ="z (cm)";
198 : //
199 0 : binsTrack[4] =10;
200 0 : xminTrack[4] =-0.8; xmaxTrack[4]=0.8; //
201 0 : axisName[4] ="sin(#phi)";
202 : //
203 0 : binsTrack[5] =10;
204 0 : xminTrack[5] =-1; xmaxTrack[5]=1; //
205 0 : axisName[5] ="tan(#theta)";
206 : //
207 0 : binsTrack[6] =40;
208 0 : xminTrack[6] =-2; xmaxTrack[6]=2; //
209 0 : axisName[6] ="1/pt (1/GeV)";
210 : //
211 0 : binsTrack[7] =50;
212 0 : xminTrack[7] =1; xmaxTrack[7]=1000; //
213 0 : axisName[7] ="pt (GeV)";
214 : //
215 0 : binsTrack[8] =18;
216 0 : xminTrack[8] =0; xmaxTrack[8]=TMath::Pi(); //
217 0 : axisName[8] ="alpha";
218 : //
219 0 : binsTrack[9] =3;
220 0 : xminTrack[9] =-0.1; xmaxTrack[9]=2.1; //
221 0 : axisName[9] ="cross";
222 : //
223 : // delta y
224 0 : xminTrack[0] =-1; xmaxTrack[0]=1; //
225 0 : fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", ndim, binsTrack,xminTrack, xmaxTrack);
226 0 : xminTrack[0] =-5; xmaxTrack[0]=5; //
227 0 : fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
228 : //
229 : // delta z
230 0 : xminTrack[0] =-1; xmaxTrack[0]=1; //
231 0 : fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", ndim, binsTrack,xminTrack, xmaxTrack);
232 0 : xminTrack[0] =-5; xmaxTrack[0]=5; //
233 0 : fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
234 : //
235 : // delta P2
236 0 : xminTrack[0] =-10; xmaxTrack[0]=10; //
237 0 : fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", ndim, binsTrack,xminTrack, xmaxTrack);
238 0 : xminTrack[0] =-5; xmaxTrack[0]=5; //
239 0 : fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
240 : //
241 : // delta P3
242 0 : xminTrack[0] =-10; xmaxTrack[0]=10; //
243 0 : fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", ndim, binsTrack,xminTrack, xmaxTrack);
244 0 : xminTrack[0] =-5; xmaxTrack[0]=5; //
245 0 : fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
246 : //
247 : // delta P4
248 0 : xminTrack[0] =-0.2; xmaxTrack[0]=0.2; //
249 0 : fHistoDelta[4] = new THnSparseS("#Delta_{1/pt} (1/GeV)","#Delta_{1/pt} (1/GeV)", ndim, binsTrack,xminTrack, xmaxTrack);
250 0 : xminTrack[0] =-5; xmaxTrack[0]=5; //
251 0 : fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
252 :
253 : //
254 : // delta Pt
255 0 : xminTrack[0] =-0.5; xmaxTrack[0]=0.5; //
256 0 : fHistoDelta[5] = new THnSparseS("#Delta_{pt}/p_{t}","#Delta_{pt}/p_{t}", ndim, binsTrack,xminTrack, xmaxTrack);
257 0 : xminTrack[0] =-5; xmaxTrack[0]=5; //
258 0 : fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
259 : //
260 :
261 0 : for (Int_t idedx=0;idedx<4;idedx++){
262 0 : xminTrack[0] =0.5; xmaxTrack[0]=1.5; //
263 0 : binsTrack[1] =40;
264 0 : xminTrack[1] =10; xmaxTrack[1]=160;
265 :
266 0 : fHistodEdxMax[idedx] = new THnSparseS(Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
267 0 : Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
268 : ndim, binsTrack,xminTrack, xmaxTrack);
269 0 : fHistodEdxTot[idedx] = new THnSparseS(Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
270 0 : Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
271 : ndim, binsTrack,xminTrack, xmaxTrack);
272 : }
273 :
274 :
275 :
276 0 : for (Int_t ivar=0;ivar<6;ivar++){
277 0 : for (Int_t ivar2=0;ivar2<ndim;ivar2++){
278 0 : fHistoDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
279 0 : fHistoDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
280 0 : fHistoPull[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
281 0 : fHistoPull[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
282 0 : BinLogX(fHistoDelta[ivar],7);
283 0 : BinLogX(fHistoPull[ivar],7);
284 0 : if (ivar<4){
285 0 : fHistodEdxMax[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
286 0 : fHistodEdxMax[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
287 0 : fHistodEdxTot[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
288 0 : fHistodEdxTot[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
289 0 : BinLogX(fHistodEdxMax[ivar],7);
290 0 : BinLogX(fHistodEdxTot[ivar],7);
291 : }
292 : }
293 : }
294 0 : }
295 :
296 :
297 : void AliTPCcalibCosmic::Add(const AliTPCcalibCosmic* cosmic){
298 : //
299 : // merge the content of the cosmic componentnts
300 : //
301 0 : for (Int_t ivar=0; ivar<6;ivar++){
302 0 : if (fHistoDelta[ivar] && cosmic->fHistoDelta[ivar]){
303 0 : fHistoDelta[ivar]->Add(cosmic->fHistoDelta[ivar]);
304 0 : }
305 0 : if (fHistoPull[ivar] && cosmic->fHistoPull[ivar]){
306 0 : fHistoPull[ivar]->Add(cosmic->fHistoPull[ivar]);
307 0 : }
308 : }
309 0 : for (Int_t ivar=0; ivar<4;ivar++){
310 0 : if (fHistodEdxMax[ivar] && cosmic->fHistodEdxMax[ivar]){
311 0 : fHistodEdxMax[ivar]->Add(cosmic->fHistodEdxMax[ivar]);
312 0 : }
313 0 : if (fHistodEdxTot[ivar] && cosmic->fHistodEdxTot[ivar]){
314 0 : fHistodEdxTot[ivar]->Add(cosmic->fHistodEdxTot[ivar]);
315 0 : }
316 : }
317 0 : if (cosmic->fCosmicTree){
318 0 : if (!fCosmicTree) {
319 0 : fCosmicTree = new TTree("pairs","pairs");
320 0 : fCosmicTree->SetDirectory(0);
321 0 : }
322 0 : AliTPCcalibCosmic::AddTree(fCosmicTree,cosmic->fCosmicTree);
323 0 : }
324 0 : }
325 :
326 :
327 :
328 :
329 : void AliTPCcalibCosmic::Process(AliESDEvent *event) {
330 : //
331 : // Process of the ESD event - fill calibration components
332 : //
333 0 : if (!event) {
334 0 : Printf("ERROR: ESD not available");
335 0 : return;
336 : }
337 :
338 : //
339 : //Int_t isOK=kTRUE;
340 : // COSMIC not signed properly
341 : // UInt_t specie = event->GetEventSpecie(); // select only cosmic events
342 : //if (specie==AliRecoParam::kCosmic || specie==AliRecoParam::kCalib) {
343 : // isOK = kTRUE;
344 : //}
345 : //if (!isOK) return;
346 : // Work around
347 0 : FindCosmicPairs(event);
348 : //const AliMultiplicity *multiplicity = event->GetMultiplicity();
349 : // Int_t ntracklets = multiplicity->GetNumberOfTracklets();
350 : //if (ntracklets>6) return; // filter out "normal" event with high multiplicity
351 : //const TString &trigger = event->GetFiredTriggerClasses();
352 : //if (trigger.Contains("C0OB0")==0) return;
353 :
354 :
355 0 : FindPairs(event); // nearly everything takes place in find pairs...
356 :
357 0 : if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
358 0 : Int_t ntracks=event->GetNumberOfTracks();
359 0 : fHistNTracks->Fill(ntracks);
360 :
361 0 : }
362 :
363 :
364 : 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 : //
366 : // par0,par1 - parameter of tracks at DCA 0
367 : // inner0,inner1 - parameter of tracks at the TPC entrance
368 : // seed0, seed1 - detailed track information
369 : // param0Combined - Use combined track parameters for binning
370 : //
371 : Int_t kMinCldEdx =20;
372 0 : Int_t ncl0 = seed0->GetNumberOfClusters();
373 0 : Int_t ncl1 = seed1->GetNumberOfClusters();
374 : const Double_t kpullCut = 10;
375 0 : Double_t x[10];
376 0 : Double_t xyz0[3];
377 0 : Double_t xyz1[3];
378 0 : par0->GetXYZ(xyz0);
379 0 : par1->GetXYZ(xyz1);
380 0 : Double_t radius0 = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]);
381 0 : Double_t radius1 = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
382 0 : inner0->GetXYZ(xyz0);
383 0 : Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
384 : // bin parameters
385 0 : x[1] = TMath::Min(ncl0,ncl1);
386 0 : x[2] = (radius0+radius1)*0.5;
387 0 : x[3] = param0Combined->GetZ();
388 0 : x[4] = inner0->GetSnp();
389 0 : x[5] = param0Combined->GetTgl();
390 0 : x[6] = param0Combined->GetSigned1Pt();
391 0 : x[7] = param0Combined->Pt();
392 0 : x[8] = alpha;
393 0 : x[9] = cross;
394 : // deltas
395 0 : Double_t delta[6];
396 0 : Double_t sigma[6];
397 0 : delta[0] = (par0->GetY()+par1->GetY());
398 0 : delta[1] = (par0->GetZ()-par1->GetZ());
399 0 : delta[2] = (par0->GetAlpha()-par1->GetAlpha()-TMath::Pi());
400 0 : delta[3] = (par0->GetTgl()+par1->GetTgl());
401 0 : delta[4] = (par0->GetParameter()[4]+par1->GetParameter()[4]);
402 0 : delta[5] = (par0->Pt()-par1->Pt())/((par0->Pt()+par1->Pt())*0.5);
403 : //
404 0 : sigma[0] = TMath::Sqrt(par0->GetSigmaY2()+par1->GetSigmaY2());
405 0 : sigma[1] = TMath::Sqrt(par0->GetSigmaZ2()+par1->GetSigmaZ2());
406 0 : sigma[2] = TMath::Sqrt(par0->GetSigmaSnp2()+par1->GetSigmaSnp2());
407 0 : sigma[3] = TMath::Sqrt(par0->GetSigmaTgl2()+par1->GetSigmaTgl2());
408 0 : sigma[4] = TMath::Sqrt(par0->GetSigma1Pt2()+par1->GetSigma1Pt2());
409 0 : sigma[5] = sigma[4]*((par0->Pt()+par1->Pt())*0.5);
410 : //
411 : Bool_t isOK = kTRUE;
412 0 : for (Int_t ivar=0;ivar<6;ivar++){
413 0 : if (sigma[ivar]==0) isOK=kFALSE;
414 0 : x[0]= delta[ivar]/sigma[ivar];
415 0 : if (TMath::Abs(x[0])>kpullCut) isOK = kFALSE;
416 : }
417 : //
418 :
419 0 : if (isOK) for (Int_t ivar=0;ivar<6;ivar++){
420 0 : x[0]= delta[ivar]; // Modifiation 10.10 use not normalized deltas
421 0 : if (ivar==2 || ivar ==3) x[0]*=1000; // angles in radian
422 0 : fHistoDelta[ivar]->Fill(x);
423 0 : if (sigma[ivar]>0){
424 0 : x[0]= delta[ivar]/sigma[ivar];
425 0 : fHistoPull[ivar]->Fill(x);
426 0 : }
427 0 : }
428 :
429 : //
430 : // Fill dedx performance
431 : //
432 0 : for (Int_t ipad=0; ipad<4;ipad++){
433 : //
434 : //
435 : //
436 : Int_t row0=0;
437 : Int_t row1=160;
438 0 : if (ipad==0) row1=63;
439 0 : if (ipad==1) {row0=63; row1=63+64;}
440 0 : if (ipad==2) {row0=128;}
441 0 : Int_t nclUp = TMath::Nint(seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
442 0 : Int_t nclDown = TMath::Nint(seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
443 0 : Int_t minCl = TMath::Min(nclUp,nclDown);
444 0 : if (minCl<kMinCldEdx) continue;
445 0 : x[1] = minCl;
446 : //
447 0 : Float_t dEdxTotUp = seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
448 0 : Float_t dEdxTotDown = seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
449 0 : Float_t dEdxMaxUp = seed0->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
450 0 : Float_t dEdxMaxDown = seed1->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
451 : //
452 0 : if (dEdxTotDown<=0) continue;
453 0 : if (dEdxMaxDown<=0) continue;
454 0 : x[0]=dEdxTotUp/dEdxTotDown;
455 0 : fHistodEdxTot[ipad]->Fill(x);
456 0 : x[0]=dEdxMaxUp/dEdxMaxDown;
457 0 : fHistodEdxMax[ipad]->Fill(x);
458 0 : }
459 :
460 :
461 :
462 0 : }
463 :
464 : void AliTPCcalibCosmic::FindPairs(const AliESDEvent *event){
465 : //
466 : // Find cosmic pairs
467 : //
468 : // Track0 is choosen in upper TPC part
469 : // Track1 is choosen in lower TPC part
470 : //
471 0 : if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
472 0 : AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
473 0 : Int_t ntracks=event->GetNumberOfTracks();
474 0 : TObjArray tpcSeeds(ntracks);
475 0 : if (ntracks==0) return;
476 0 : Double_t vtxx[3]={0,0,0};
477 0 : Double_t svtxx[3]={0.000001,0.000001,100.};
478 0 : AliESDVertex vtx(vtxx,svtxx);
479 : //
480 : //track loop
481 : //
482 0 : for (Int_t i=0;i<ntracks;++i) {
483 0 : AliESDtrack *track = event->GetTrack(i);
484 0 : fClusters->Fill(track->GetTPCNcls());
485 :
486 0 : const AliExternalTrackParam * trackIn = track->GetInnerParam();
487 0 : const AliExternalTrackParam * trackOut = track->GetOuterParam();
488 0 : if (!trackIn) continue;
489 0 : if (!trackOut) continue;
490 0 : if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // filter laser
491 :
492 :
493 0 : AliESDfriendTrack *friendTrack = (AliESDfriendTrack*)track->GetFriendTrack();
494 0 : if (!friendTrack) continue;
495 : TObject *calibObject;
496 : AliTPCseed *seed = 0;
497 0 : for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
498 0 : if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
499 : }
500 0 : if (seed) tpcSeeds.AddAt(seed,i);
501 :
502 0 : Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
503 0 : if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
504 0 : fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,kMaxRow));
505 : //
506 0 : if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,kMaxRow));
507 : //
508 : // if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,kMaxRow)>300) {
509 : // //TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
510 : // //if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
511 : // // if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
512 : // }
513 :
514 : }
515 :
516 0 : }
517 :
518 0 : if (ntracks<2) return;
519 : //
520 : // Find pairs
521 : //
522 0 : for (Int_t i=0;i<ntracks;++i) {
523 0 : AliESDtrack *track0 = event->GetTrack(i);
524 : // track0 - choosen upper part
525 0 : if (!track0) continue;
526 0 : if (!track0->GetOuterParam()) continue;
527 0 : if (track0->GetOuterParam()->GetAlpha()<0) continue;
528 0 : Double_t dir0[3];
529 0 : track0->GetDirection(dir0);
530 0 : for (Int_t j=0;j<ntracks;++j) {
531 0 : if (i==j) continue;
532 0 : AliESDtrack *track1 = event->GetTrack(j);
533 : //track 1 lower part
534 0 : if (!track1) continue;
535 0 : if (!track1->GetOuterParam()) continue;
536 0 : if (track1->GetOuterParam()->GetAlpha()>0) continue;
537 : //
538 0 : Double_t dir1[3];
539 0 : track1->GetDirection(dir1);
540 :
541 0 : AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
542 0 : AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
543 0 : if (! seed0) continue;
544 0 : if (! seed1) continue;
545 0 : Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,kMaxRow);
546 0 : Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,kMaxRow);
547 : //
548 0 : Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63);
549 0 : Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63);
550 : //
551 0 : Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,kMaxRow);
552 0 : Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,kMaxRow);
553 : //
554 0 : Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
555 0 : Float_t d0 = track0->GetLinearD(0,0);
556 0 : Float_t d1 = track1->GetLinearD(0,0);
557 : //
558 : // conservative cuts - convergence to be guarantied
559 : // applying before track propagation
560 0 : if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0
561 0 : if (dir>fCutMinDir) continue; // direction vector product
562 0 : Float_t bz = AliTracker::GetBz();
563 0 : Float_t dvertex0[2]; //distance to 0,0
564 0 : Float_t dvertex1[2]; //distance to 0,0
565 0 : track0->GetDZ(0,0,0,bz,dvertex0);
566 0 : track1->GetDZ(0,0,0,bz,dvertex1);
567 0 : if (TMath::Abs(dvertex0[1])>250) continue;
568 0 : if (TMath::Abs(dvertex1[1])>250) continue;
569 : //
570 : //
571 : //
572 0 : Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
573 0 : AliExternalTrackParam param0(*track0);
574 0 : AliExternalTrackParam param1(*track1);
575 : //
576 : // Propagate using Magnetic field and correct fo material budget
577 : //
578 : Double_t sign0=-1;
579 : Double_t sign1=1;
580 : Double_t maxsnp=0.90;
581 0 : AliTracker::PropagateTrackToBxByBz(¶m0,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign0);
582 0 : AliTracker::PropagateTrackToBxByBz(¶m1,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign1);
583 : //
584 : // Propagate rest to the 0,0 DCA - z should be ignored
585 : //
586 0 : Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
587 0 : Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
588 : //
589 0 : param0.GetDZ(0,0,0,bz,dvertex0);
590 0 : param1.GetDZ(0,0,0,bz,dvertex1);
591 0 : if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
592 : //
593 0 : Double_t xyz0[3];//,pxyz0[3];
594 0 : Double_t xyz1[3];//,pxyz1[3];
595 0 : param0.GetXYZ(xyz0);
596 0 : param1.GetXYZ(xyz1);
597 0 : Bool_t isPair = IsPair(¶m0,¶m1);
598 : //
599 0 : if (isPair) FillAcordeHist(track0);
600 0 : if (isPair &¶m0.Pt()>1) {
601 0 : const TString &trigger = event->GetFiredTriggerClasses();
602 0 : UInt_t specie = event->GetEventSpecie();
603 0 : printf("COSMIC ?\t%s\t%d\t%f\t%f\n", trigger.Data(),specie, param0.GetZ(), param1.GetZ());
604 0 : }
605 : //
606 : // combined track params
607 : //
608 0 : AliExternalTrackParam *par0U=MakeCombinedTrack(¶m0,¶m1);
609 0 : AliExternalTrackParam *par1U=MakeCombinedTrack(¶m1,¶m0);
610 :
611 : //
612 0 : if (fStreamLevel>0){
613 0 : TTreeSRedirector * cstream = GetDebugStreamer();
614 : //printf("My stream=%p\n",(void*)cstream);
615 0 : AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
616 0 : AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
617 0 : AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
618 0 : AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
619 0 : Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
620 0 : Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
621 0 : Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
622 0 : Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
623 : //
624 : //
625 : //
626 : Int_t cross =0; // 0 no cross, 2 cross on both sides
627 0 : if (isCrossI) cross+=1;
628 0 : if (isCrossO) cross+=1;
629 0 : FillHistoPerformance(¶m0, ¶m1, ip0, ip1, seed0, seed1,par0U, cross);
630 0 : if (cstream) {
631 0 : (*cstream) << "Track0" <<
632 0 : "run="<<fRun<< // run number
633 0 : "event="<<fEvent<< // event number
634 0 : "time="<<fTime<< // time stamp of event
635 0 : "trigger="<<fTrigger<< // trigger
636 0 : "triggerClass="<<&fTriggerClass<< // trigger
637 0 : "mag="<<fMagF<< // magnetic field
638 0 : "dir="<<dir<< // direction
639 0 : "OK="<<isPair<< // will be accepted
640 0 : "b0="<<b0<< // propagate status
641 0 : "b1="<<b1<< // propagate status
642 0 : "crossI="<<isCrossI<< // cross inner
643 0 : "crossO="<<isCrossO<< // cross outer
644 : //
645 0 : "Orig0.=" << track0 << // original track 0
646 0 : "Orig1.=" << track1 << // original track 1
647 0 : "Tr0.="<<¶m0<< // track propagated to the DCA 0,0
648 0 : "Tr1.="<<¶m1<< // track propagated to the DCA 0,0
649 0 : "Ip0.="<<ip0<< // inner param - upper
650 0 : "Ip1.="<<ip1<< // inner param - lower
651 0 : "Op0.="<<op0<< // outer param - upper
652 0 : "Op1.="<<op1<< // outer param - lower
653 0 : "Up0.="<<par0U<< // combined track 0
654 0 : "Up1.="<<par1U<< // combined track 1
655 : //
656 0 : "v00="<<dvertex0[0]<< // distance using kalman
657 0 : "v01="<<dvertex0[1]<< //
658 0 : "v10="<<dvertex1[0]<< //
659 0 : "v11="<<dvertex1[1]<< //
660 0 : "d0="<<d0<< // linear distance to 0,0
661 0 : "d1="<<d1<< // linear distance to 0,0
662 : //
663 : //
664 : //
665 0 : "x00="<<xyz0[0]<< // global position close to vertex
666 0 : "x01="<<xyz0[1]<<
667 0 : "x02="<<xyz0[2]<<
668 : //
669 0 : "x10="<<xyz1[0]<< // global position close to vertex
670 0 : "x11="<<xyz1[1]<<
671 0 : "x12="<<xyz1[2]<<
672 : //
673 0 : "alpha0="<<alpha0<<
674 0 : "alpha1="<<alpha1<<
675 0 : "dir00="<<dir0[0]<< // direction upper
676 0 : "dir01="<<dir0[1]<<
677 0 : "dir02="<<dir0[2]<<
678 : //
679 0 : "dir10="<<dir1[0]<< // direction lower
680 0 : "dir11="<<dir1[1]<<
681 0 : "dir12="<<dir1[2]<<
682 : //
683 : //
684 0 : "Seed0.=" << seed0 << // original seed 0
685 0 : "Seed1.=" << seed1 << // original seed 1
686 : //
687 0 : "dedx0="<<dedx0<< // dedx0 - all
688 0 : "dedx1="<<dedx1<< // dedx1 - all
689 : //
690 0 : "dedx0I="<<dedx0I<< // dedx0 - inner ROC
691 0 : "dedx1I="<<dedx1I<< // dedx1 - inner ROC
692 : //
693 0 : "dedx0O="<<dedx0O<< // dedx0 - outer ROC
694 0 : "dedx1O="<<dedx1O<< // dedx1 - outer ROC
695 : "\n";
696 : }
697 0 : }
698 0 : delete par0U;
699 0 : delete par1U;
700 0 : }
701 0 : }
702 0 : }
703 :
704 :
705 :
706 :
707 : void AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
708 :
709 : // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
710 0 : if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
711 :
712 : const Double_t acordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
713 : const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
714 :
715 0 : Double_t r[3];
716 0 : upperTrack->GetXYZ(r);
717 0 : Double_t d[3];
718 0 : upperTrack->GetDirection(d);
719 : Double_t x,z;
720 0 : z = r[2] + (d[2]/d[1])*(acordePlane - r[1]);
721 0 : x = r[0] + (d[0]/d[1])*(acordePlane - r[1]);
722 :
723 0 : if (x > roof) {
724 0 : x = r[0] + (d[0]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
725 0 : z = r[2] + (d[2]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
726 0 : }
727 0 : if (x < -roof) {
728 0 : x = r[0] + (d[0]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
729 0 : z = r[2] + (d[2]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
730 0 : }
731 :
732 0 : fModules->Fill(z, x);
733 :
734 0 : }
735 :
736 :
737 :
738 : Long64_t AliTPCcalibCosmic::Merge(TCollection *const li) {
739 : //
740 : // component merging
741 : //
742 :
743 0 : TIterator* iter = li->MakeIterator();
744 : AliTPCcalibCosmic* cal = 0;
745 :
746 0 : while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
747 0 : if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
748 : //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
749 0 : return -1;
750 : }
751 :
752 0 : fHistNTracks->Add(cal->GetHistNTracks());
753 0 : fClusters->Add(cal-> GetHistClusters());
754 0 : fModules->Add(cal->GetHistAcorde());
755 0 : fHistPt->Add(cal->GetHistPt());
756 0 : fDeDx->Add(cal->GetHistDeDx());
757 0 : fDeDxMIP->Add(cal->GetHistMIP());
758 0 : Add(cal);
759 : }
760 0 : return 0;
761 :
762 0 : }
763 :
764 :
765 : Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1) const{
766 : //
767 : //
768 : /*
769 : // 0. Same direction - OPOSITE - cutDir +cutT
770 : TCut cutDir("cutDir","dir<-0.99")
771 : // 1.
772 : TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
773 : //
774 : // 2. The same rphi
775 : TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
776 : //
777 : //
778 : //
779 : TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
780 : // 1/Pt diff cut
781 : */
782 0 : const Double_t *p0 = tr0->GetParameter();
783 0 : const Double_t *p1 = tr1->GetParameter();
784 0 : if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
785 0 : if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
786 0 : if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
787 :
788 0 : Double_t d0[3], d1[3];
789 0 : tr0->GetDirection(d0);
790 0 : tr1->GetDirection(d1);
791 0 : if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
792 : //
793 0 : return kTRUE;
794 0 : }
795 :
796 :
797 :
798 : Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
799 : //
800 : // Calculate the MIP value - gaussian fit used
801 : //
802 :
803 0 : TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
804 0 : funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
805 0 : hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
806 0 : hist->Fit(funcDoubleGaus);
807 0 : Double_t aMIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
808 :
809 0 : delete funcDoubleGaus;
810 :
811 0 : return aMIPvalue;
812 :
813 0 : }
814 :
815 :
816 :
817 :
818 : void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) {
819 : //
820 : // Not implemented yet
821 : //
822 0 : return;
823 :
824 : }
825 :
826 :
827 : void AliTPCcalibCosmic::BinLogX(THnSparse *const h, Int_t axisDim) {
828 :
829 : // Method for the correct logarithmic binning of histograms
830 :
831 0 : TAxis *axis = h->GetAxis(axisDim);
832 0 : int bins = axis->GetNbins();
833 :
834 0 : Double_t from = axis->GetXmin();
835 0 : Double_t to = axis->GetXmax();
836 0 : Double_t *newBins = new Double_t[bins + 1];
837 :
838 0 : newBins[0] = from;
839 0 : Double_t factor = pow(to/from, 1./bins);
840 :
841 0 : for (int i = 1; i <= bins; i++) {
842 0 : newBins[i] = factor * newBins[i-1];
843 : }
844 0 : axis->Set(bins, newBins);
845 0 : delete [] newBins;
846 :
847 0 : }
848 :
849 :
850 : void AliTPCcalibCosmic::BinLogX(TH1 *const h) {
851 :
852 : // Method for the correct logarithmic binning of histograms
853 :
854 0 : TAxis *axis = h->GetXaxis();
855 0 : int bins = axis->GetNbins();
856 :
857 0 : Double_t from = axis->GetXmin();
858 0 : Double_t to = axis->GetXmax();
859 0 : Double_t *newBins = new Double_t[bins + 1];
860 :
861 0 : newBins[0] = from;
862 0 : Double_t factor = pow(to/from, 1./bins);
863 :
864 0 : for (int i = 1; i <= bins; i++) {
865 0 : newBins[i] = factor * newBins[i-1];
866 : }
867 0 : axis->Set(bins, newBins);
868 0 : delete [] newBins;
869 :
870 0 : }
871 :
872 :
873 : AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
874 : //
875 : // Make a atrack using the kalman update of track0 and track1
876 : //
877 0 : AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
878 0 : par1R->Rotate(track0->GetAlpha());
879 0 : par1R->PropagateTo(track0->GetX(),AliTracker::GetBz());
880 : //
881 : //
882 0 : Double_t * param = (Double_t*)par1R->GetParameter();
883 0 : Double_t * covar = (Double_t*)par1R->GetCovariance();
884 :
885 0 : param[0]*=1; //OK
886 0 : param[1]*=1; //OK
887 0 : param[2]*=1; //?
888 0 : param[3]*=-1; //OK
889 0 : param[4]*=-1; //OK
890 : //
891 0 : covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
892 : //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
893 0 : covar[13]*=-1.;
894 0 : return par1R;
895 0 : }
896 :
897 : AliExternalTrackParam *AliTPCcalibCosmic::MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
898 : //
899 : // Make combined track
900 : //
901 : //
902 0 : AliExternalTrackParam * par1T = MakeTrack(track0,track1);
903 0 : AliExternalTrackParam * par0U = new AliExternalTrackParam(*track0);
904 : //
905 0 : UpdateTrack(*par0U,*par1T);
906 0 : delete par1T;
907 0 : return par0U;
908 0 : }
909 :
910 :
911 : void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
912 : //
913 : // Update track 1 with track 2
914 : //
915 : //
916 : //
917 0 : TMatrixD vecXk(5,1); // X vector
918 0 : TMatrixD covXk(5,5); // X covariance
919 0 : TMatrixD matHk(5,5); // vector to mesurement
920 0 : TMatrixD measR(5,5); // measurement error
921 0 : TMatrixD vecZk(5,1); // measurement
922 : //
923 0 : TMatrixD vecYk(5,1); // Innovation or measurement residual
924 0 : TMatrixD matHkT(5,5);
925 0 : TMatrixD matSk(5,5); // Innovation (or residual) covariance
926 0 : TMatrixD matKk(5,5); // Optimal Kalman gain
927 0 : TMatrixD mat1(5,5); // update covariance matrix
928 0 : TMatrixD covXk2(5,5); //
929 0 : TMatrixD covOut(5,5);
930 : //
931 0 : Double_t *param1=(Double_t*) track1.GetParameter();
932 0 : Double_t *covar1=(Double_t*) track1.GetCovariance();
933 0 : Double_t *param2=(Double_t*) track2.GetParameter();
934 0 : Double_t *covar2=(Double_t*) track2.GetCovariance();
935 : //
936 : // copy data to the matrix
937 0 : for (Int_t ipar=0; ipar<5; ipar++){
938 0 : for (Int_t jpar=0; jpar<5; jpar++){
939 0 : covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
940 0 : measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
941 0 : matHk(ipar,jpar)=0;
942 0 : mat1(ipar,jpar)=0;
943 : }
944 0 : vecXk(ipar,0) = param1[ipar];
945 0 : vecZk(ipar,0) = param2[ipar];
946 0 : matHk(ipar,ipar)=1;
947 0 : mat1(ipar,ipar)=0;
948 : }
949 : //
950 : //
951 : //
952 : //
953 : //
954 0 : vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
955 0 : matHkT=matHk.T(); matHk.T();
956 0 : matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
957 0 : matSk.Invert();
958 0 : matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
959 0 : vecXk += matKk*vecYk; // updated vector
960 0 : covXk2 = (mat1-(matKk*matHk));
961 0 : covOut = covXk2*covXk;
962 : //
963 : //
964 : //
965 : // copy from matrix to parameters
966 : if (0) {
967 : vecXk.Print();
968 : vecZk.Print();
969 : //
970 : measR.Print();
971 : covXk.Print();
972 : covOut.Print();
973 : //
974 : track1.Print();
975 : track2.Print();
976 : }
977 :
978 0 : for (Int_t ipar=0; ipar<5; ipar++){
979 0 : param1[ipar]= vecXk(ipar,0) ;
980 0 : for (Int_t jpar=0; jpar<5; jpar++){
981 0 : covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
982 : }
983 : }
984 0 : }
985 :
986 :
987 :
988 : void AliTPCcalibCosmic::FindCosmicPairs(const AliESDEvent * event) {
989 : //
990 : // find cosmic pairs trigger by random trigger
991 : //
992 : //
993 0 : AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD();
994 0 : AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC();
995 0 : AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
996 : const Double_t kMinPt=1;
997 : const Double_t kMinPtMax=0.8;
998 : const Double_t kMinNcl=50;
999 : const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};
1000 0 : Int_t ntracks=event->GetNumberOfTracks();
1001 : // Float_t dcaTPC[2]={0,0};
1002 : // Float_t covTPC[3]={0,0,0};
1003 :
1004 0 : UInt_t specie = event->GetEventSpecie(); // skip laser events
1005 0 : if (specie==AliRecoParam::kCalib) return;
1006 :
1007 :
1008 :
1009 0 : for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
1010 0 : AliESDtrack *track0 = event->GetTrack(itrack0);
1011 0 : if (!track0) continue;
1012 0 : if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;
1013 :
1014 0 : if (TMath::Abs(AliTracker::GetBz())>1&&track0->Pt()<kMinPt) continue;
1015 0 : if (track0->GetTPCncls()<kMinNcl) continue;
1016 0 : if (TMath::Abs(track0->GetY())<kMaxDelta[0]) continue;
1017 0 : if (track0->GetKinkIndex(0)>0) continue;
1018 0 : const Double_t * par0=track0->GetParameter(); //track param at rhe DCA
1019 : //rm primaries
1020 : //
1021 : //track0->GetImpactParametersTPC(dcaTPC,covTPC);
1022 : //if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
1023 : //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
1024 : // const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();
1025 0 : for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
1026 0 : AliESDtrack *track1 = event->GetTrack(itrack1);
1027 0 : if (!track1) continue;
1028 0 : if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;
1029 0 : if (track1->GetKinkIndex(0)>0) continue;
1030 0 : if (TMath::Abs(AliTracker::GetBz())>1&&track1->Pt()<kMinPt) continue;
1031 0 : if (track1->GetTPCncls()<kMinNcl) continue;
1032 0 : if (TMath::Abs(AliTracker::GetBz())>1&&TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
1033 0 : if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;
1034 : //track1->GetImpactParametersTPC(dcaTPC,covTPC);
1035 : // if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
1036 : //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
1037 : //
1038 0 : const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
1039 : //
1040 : Bool_t isPair=kTRUE;
1041 0 : for (Int_t ipar=0; ipar<5; ipar++){
1042 0 : if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
1043 0 : if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
1044 : }
1045 0 : if (!isPair) continue;
1046 0 : if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
1047 : //delta with correct sign
1048 : /*
1049 : TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"
1050 : TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"
1051 : TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"
1052 : */
1053 0 : if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign
1054 0 : if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
1055 0 : if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
1056 0 : if (!isPair) continue;
1057 0 : TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1058 0 : Int_t eventNumber = event->GetEventNumberInFile();
1059 0 : Bool_t hasFriend=(esdFriend) ? track0->GetFriendTrack() : 0;//( esdFriend->GetTrack(itrack0)!=0):0;
1060 0 : Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4);
1061 0 : printf("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS);
1062 :
1063 :
1064 : // const AliExternalTrackParam * trackIn1 = track1->GetInnerParam();
1065 : //
1066 : //
1067 0 : TTreeSRedirector * pcstream = GetDebugStreamer();
1068 0 : Int_t ntracksSPD = vertexSPD->GetNContributors();
1069 0 : Int_t ntracksTPC = vertexTPC->GetNContributors();
1070 : //
1071 0 : AliESDfriendTrack *friendTrack0 = (AliESDfriendTrack*)track0->GetFriendTrack();//esdFriend->GetTrack(itrack0);
1072 0 : if (!friendTrack0) continue;
1073 0 : AliESDfriendTrack *friendTrack1 = (AliESDfriendTrack*)track1->GetFriendTrack();//esdFriend->GetTrack(itrack1);
1074 0 : if (!friendTrack1) continue;
1075 : TObject *calibObject;
1076 : AliTPCseed *seed0 = 0;
1077 : AliTPCseed *seed1 = 0;
1078 : //
1079 0 : for (Int_t l=0;(calibObject=friendTrack0->GetCalibObject(l));++l) {
1080 0 : if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
1081 : }
1082 0 : for (Int_t l=0;(calibObject=friendTrack1->GetCalibObject(l));++l) {
1083 0 : if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
1084 : }
1085 : //
1086 0 : if (pcstream){
1087 0 : (*pcstream)<<"pairs"<<
1088 0 : "run="<<fRun<< // run number
1089 0 : "event="<<fEvent<< // event number
1090 0 : "time="<<fTime<< // time stamp of event
1091 0 : "trigger="<<fTrigger<< // trigger
1092 0 : "triggerClass="<<&fTriggerClass<< // trigger
1093 0 : "mag="<<fMagF<< // magnetic field
1094 : //
1095 0 : "nSPD="<<ntracksSPD<<
1096 0 : "nTPC="<<ntracksTPC<<
1097 0 : "vSPD.="<<vertexSPD<< //primary vertex -SPD
1098 0 : "vTPC.="<<vertexTPC<< //primary vertex -TPC
1099 0 : "t0.="<<track0<< //track0
1100 0 : "t1.="<<track1<< //track1
1101 0 : "ft0.="<<friendTrack0<< //track0
1102 0 : "ft1.="<<friendTrack1<< //track1
1103 0 : "s0.="<<seed0<< //track0
1104 0 : "s1.="<<seed1<< //track1
1105 : "\n";
1106 : }
1107 0 : if (!fCosmicTree) {
1108 0 : fCosmicTree = new TTree("pairs","pairs");
1109 0 : fCosmicTree->SetDirectory(0);
1110 : }
1111 0 : if (fCosmicTree->GetEntries()==0){
1112 : //
1113 0 : fCosmicTree->SetDirectory(0);
1114 0 : fCosmicTree->Branch("t0.",&track0);
1115 0 : fCosmicTree->Branch("t1.",&track1);
1116 0 : fCosmicTree->Branch("ft0.",&friendTrack0);
1117 0 : fCosmicTree->Branch("ft1.",&friendTrack1);
1118 : }else{
1119 0 : fCosmicTree->SetBranchAddress("t0.",&track0);
1120 0 : fCosmicTree->SetBranchAddress("t1.",&track1);
1121 0 : fCosmicTree->SetBranchAddress("ft0.",&friendTrack0);
1122 0 : fCosmicTree->SetBranchAddress("ft1.",&friendTrack1);
1123 : }
1124 0 : fCosmicTree->Fill();
1125 0 : }
1126 0 : }
1127 0 : }
1128 :
1129 :
1130 : void AliTPCcalibCosmic::Terminate(){
1131 : //
1132 : // copy the cosmic tree to memory resident tree
1133 : //
1134 : static Int_t counter=0;
1135 0 : printf("AliTPCcalibCosmic::Terminate\t%d\n",counter);
1136 0 : counter++;
1137 0 : AliTPCcalibBase::Terminate();
1138 0 : }
1139 :
1140 :
1141 : void AliTPCcalibCosmic::AddTree(TTree* treeOutput, TTree * treeInput){
1142 : //
1143 : // Add the content of tree:
1144 : // Notice automatic copy of tree in ROOT does not work for such complicated tree
1145 : //
1146 0 : return;
1147 : //if (TMath::Abs(fMagF)<0.1) return; // work around - otherwise crashes
1148 : AliESDtrack *track0=new AliESDtrack;
1149 : AliESDtrack *track1=new AliESDtrack;
1150 : AliESDfriendTrack *ftrack0=new AliESDfriendTrack;
1151 : AliESDfriendTrack *ftrack1=new AliESDfriendTrack;
1152 : treeInput->SetBranchAddress("t0.",&track0);
1153 : treeInput->SetBranchAddress("t1.",&track1);
1154 : treeInput->SetBranchAddress("ft0.",&ftrack0);
1155 : treeInput->SetBranchAddress("ft1.",&ftrack1);
1156 : treeOutput->SetDirectory(0);
1157 : //
1158 : Int_t entries= treeInput->GetEntries();
1159 : Int_t step=1+Int_t(TMath::Log(1+treeOutput->GetEntries()/10.));
1160 : for (Int_t i=0; i<entries; i+=step){
1161 : treeInput->SetBranchAddress("t0.",&track0);
1162 : treeInput->SetBranchAddress("t1.",&track1);
1163 : treeInput->SetBranchAddress("ft0.",&ftrack0);
1164 : treeInput->SetBranchAddress("ft1.",&ftrack1);
1165 : treeInput->GetEntry(i);
1166 : if (!track0) continue;
1167 : if (!track1) continue;
1168 : if (!ftrack0) continue;
1169 : if (!ftrack1) continue;
1170 : if (track0->GetTPCncls()<=0) continue;
1171 : if (track1->GetTPCncls()<=0) continue;
1172 : if (!track0->GetInnerParam()) continue;
1173 : if (!track1->GetInnerParam()) continue;
1174 : if (!track0->GetTPCInnerParam()) continue;
1175 : if (!track1->GetTPCInnerParam()) continue;
1176 : //track0
1177 : treeOutput->SetBranchAddress("t0.",&track0);
1178 : treeOutput->SetBranchAddress("t1.",&track1);
1179 : treeOutput->SetBranchAddress("ft0.",&ftrack0);
1180 : treeOutput->SetBranchAddress("ft1.",&ftrack1);
1181 : treeOutput->Fill();
1182 : delete track0;
1183 : delete track1;
1184 : delete ftrack0;
1185 : delete ftrack1;
1186 : track0=0;
1187 : track1=0;
1188 : ftrack0=0;
1189 : ftrack1=0;
1190 : }
1191 : }
1192 :
1193 :
1194 :
1195 : void AliTPCcalibCosmic::MakeFitTree(TTree * treeInput, TTreeSRedirector *pcstream, const TObjArray * corrArray, Int_t step, Int_t run){
1196 : //
1197 : // Make fit tree
1198 : // refit the tracks with original points + corrected points for each correction
1199 : // Input:
1200 : // treeInput - tree with cosmic tracks
1201 : // pcstream - debug output
1202 :
1203 : // Algorithm:
1204 : // Loop over pair of cosmic tracks:
1205 : // 1. Find trigger offset between cosmic event and "physic" trigger
1206 : // 2. Refit tracks with current transformation
1207 : // 3. Refit tracks using additional "primitive" distortion on top
1208 : // Best correction estimated as linear combination of corrections
1209 : // minimizing the observed distortions
1210 : // Observed distortions - matching in the y,z, snp, theta and 1/pt
1211 : //
1212 : const Double_t kResetCov=20.;
1213 : const Double_t kMaxDelta[5]={1,40,0.03,0.01,0.2};
1214 : const Double_t kSigma=2.;
1215 : const Double_t kMaxTime=1050;
1216 : const Double_t kMaxSnp=0.8;
1217 0 : Int_t ncorr=corrArray->GetEntries();
1218 0 : AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1219 0 : AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
1220 0 : AliGRPObject* grp = AliTPCcalibDB::Instance()->GetGRP(run);
1221 0 : Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
1222 0 : transform->SetCurrentRun(run);
1223 0 : transform->SetCurrentTimeStamp(TMath::Nint(time));
1224 0 : Double_t covar[15];
1225 0 : for (Int_t i=0;i<15;i++) covar[i]=0;
1226 0 : covar[0]=kSigma*kSigma;
1227 0 : covar[2]=kSigma*kSigma;
1228 0 : covar[5]=kSigma*kSigma/Float_t(150*150);
1229 0 : covar[9]=kSigma*kSigma/Float_t(150*150);
1230 0 : covar[14]=0.2*0.2;
1231 0 : Double_t *distortions = new Double_t[ncorr+1];
1232 :
1233 0 : AliESDtrack *track0=new AliESDtrack;
1234 0 : AliESDtrack *track1=new AliESDtrack;
1235 0 : AliESDfriendTrack *ftrack0=new AliESDfriendTrack;
1236 0 : AliESDfriendTrack *ftrack1=new AliESDfriendTrack;
1237 0 : treeInput->SetBranchAddress("t0.",&track0);
1238 0 : treeInput->SetBranchAddress("t1.",&track1);
1239 0 : treeInput->SetBranchAddress("ft0.",&ftrack0);
1240 0 : treeInput->SetBranchAddress("ft1.",&ftrack1);
1241 0 : Int_t entries= treeInput->GetEntries();
1242 0 : for (Int_t i=0; i<entries; i+=step){
1243 0 : treeInput->GetEntry(i);
1244 0 : if (i%20==0) printf("%d\n",i);
1245 0 : if (!ftrack0->GetTPCOut()) continue;
1246 0 : if (!ftrack1->GetTPCOut()) continue;
1247 : AliTPCseed *seed0=0;
1248 : AliTPCseed *seed1=0;
1249 : TObject *calibObject;
1250 0 : for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
1251 0 : if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
1252 : }
1253 0 : for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
1254 0 : if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
1255 : }
1256 0 : if (!seed0) continue;
1257 0 : if (!seed1) continue;
1258 0 : if (TMath::Abs(seed0->GetSnp())>kMaxSnp) continue;
1259 0 : if (TMath::Abs(seed1->GetSnp())>kMaxSnp) continue;
1260 : //
1261 : //
1262 : Int_t nclA0=0, nclC0=0; // number of clusters
1263 : Int_t nclA1=0, nclC1=0; // number of clusters
1264 0 : Int_t ncl0=0,ncl1=0;
1265 : Double_t rmin0=300, rmax0=-300; // variables to estimate the time0 - trigger offset
1266 : Double_t rmin1=300, rmax1=-300;
1267 : Double_t tmin0=2000, tmax0=-2000;
1268 : Double_t tmin1=2000, tmax1=-2000;
1269 : //
1270 : //
1271 : // calculate trigger offset usig "missing clusters"
1272 0 : for (Int_t irow=0; irow<kMaxRow; irow++){
1273 0 : AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
1274 0 : if (cluster0 &&cluster0->GetX()>10){
1275 0 : if (cluster0->GetX()<rmin0) { rmin0=cluster0->GetX(); tmin0=cluster0->GetTimeBin();}
1276 0 : if (cluster0->GetX()>rmax0) { rmax0=cluster0->GetX(); tmax0=cluster0->GetTimeBin();}
1277 0 : ncl0++;
1278 0 : if (cluster0->GetDetector()%36<18) nclA0++;
1279 0 : if (cluster0->GetDetector()%36>=18) nclC0++;
1280 : }
1281 0 : AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
1282 0 : if (cluster1&&cluster1->GetX()>10){
1283 0 : if (cluster1->GetX()<rmin1) { rmin1=cluster1->GetX(); tmin1=cluster1->GetTimeBin();}
1284 0 : if (cluster1->GetX()>rmax1) { rmax1=cluster1->GetX(); tmax1=cluster1->GetTimeBin();}
1285 0 : ncl1++;
1286 0 : if (cluster1->GetDetector()%36<18) nclA1++;
1287 0 : if (cluster1->GetDetector()%36>=18) nclC1++;
1288 : }
1289 : }
1290 0 : Int_t cosmicType=0; // types of cosmic topology
1291 0 : if ((nclA0>nclC0) && (nclA1>nclC1)) cosmicType=0; // AA side
1292 0 : if ((nclA0<nclC0) && (nclA1<nclC1)) cosmicType=1; // CC side
1293 0 : if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=2; // AC side
1294 0 : if ((nclA0<nclC0) && (nclA1>nclC1)) cosmicType=3; // CA side
1295 : //if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=6; // AC side out of time
1296 : //if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=7; // CA side out of time
1297 : //
1298 0 : Double_t deltaTime = 0; // correction for trigger not in time - equalizing the time arival
1299 0 : if ((cosmicType>1)&&TMath::Abs(track1->GetZ()-track0->GetZ())>4){
1300 0 : cosmicType+=4;
1301 0 : deltaTime=0.5*(track1->GetZ()-track0->GetZ())/param->GetZWidth();
1302 0 : if (nclA0>nclC0) deltaTime*=-1; // if A side track
1303 : }
1304 : //
1305 0 : TVectorD vectorDT(8);
1306 0 : Int_t crossCounter=0;
1307 0 : Double_t deltaTimeCross = AliTPCcalibCosmic::GetDeltaTime(rmin0, rmax0, rmin1, rmax1, tmin0, tmax0, tmin1, tmax1, TMath::Abs(track0->GetY()),vectorDT);
1308 0 : Bool_t isOKTrigger=kTRUE;
1309 0 : for (Int_t ic=0; ic<6;ic++) {
1310 0 : if (TMath::Abs(vectorDT[ic])>0) {
1311 0 : if (vectorDT[ic]+vectorDT[6]<0) isOKTrigger=kFALSE;
1312 0 : if (vectorDT[ic]+vectorDT[7]>kMaxTime) isOKTrigger=kFALSE;
1313 0 : if (isOKTrigger){
1314 0 : crossCounter++;
1315 0 : }
1316 : }
1317 : }
1318 0 : Double_t deltaTimeCluster=deltaTime;
1319 0 : if ((cosmicType==0 || cosmicType==1) && crossCounter>0){
1320 0 : deltaTimeCluster=deltaTimeCross;
1321 0 : cosmicType+=8;
1322 0 : }
1323 0 : if (nclA0*nclC0>0 || nclA1*nclC1>0) cosmicType+=16; // mixed A side C side - bad for visualization
1324 : //
1325 : // Apply current transformation
1326 : //
1327 : //
1328 0 : for (Int_t irow=0; irow<kMaxRow; irow++){
1329 0 : AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
1330 0 : if (cluster0 &&cluster0->GetX()>10){
1331 0 : Double_t x0[3]={ static_cast<Double_t>(cluster0->GetRow()),cluster0->GetPad(),cluster0->GetTimeBin()+deltaTimeCluster};
1332 0 : Int_t index0[1]={cluster0->GetDetector()};
1333 0 : transform->Transform(x0,index0,0,1);
1334 0 : cluster0->SetX(x0[0]);
1335 0 : cluster0->SetY(x0[1]);
1336 0 : cluster0->SetZ(x0[2]);
1337 : //
1338 0 : }
1339 0 : AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
1340 0 : if (cluster1&&cluster1->GetX()>10){
1341 0 : Double_t x1[3]={ static_cast<Double_t>(cluster1->GetRow()),cluster1->GetPad(),cluster1->GetTimeBin()+deltaTimeCluster};
1342 0 : Int_t index1[1]={cluster1->GetDetector()};
1343 0 : transform->Transform(x1,index1,0,1);
1344 0 : cluster1->SetX(x1[0]);
1345 0 : cluster1->SetY(x1[1]);
1346 0 : cluster1->SetZ(x1[2]);
1347 0 : }
1348 : }
1349 : //
1350 : //
1351 0 : Double_t alpha=track0->GetAlpha(); // rotation frame
1352 0 : Double_t cos = TMath::Cos(alpha);
1353 0 : Double_t sin = TMath::Sin(alpha);
1354 0 : Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
1355 0 : AliExternalTrackParam btrack0=*(ftrack0->GetTPCOut());
1356 0 : AliExternalTrackParam btrack1=*(ftrack1->GetTPCOut());
1357 0 : btrack0.Rotate(alpha);
1358 0 : btrack1.Rotate(alpha);
1359 : // change the sign for track 1
1360 0 : Double_t* par1=(Double_t*)btrack0.GetParameter();
1361 0 : par1[3]*=-1;
1362 0 : par1[4]*=-1;
1363 0 : btrack0.AddCovariance(covar);
1364 0 : btrack1.AddCovariance(covar);
1365 0 : btrack0.ResetCovariance(kResetCov);
1366 0 : btrack1.ResetCovariance(kResetCov);
1367 0 : Bool_t isOK=kTRUE;
1368 0 : Bool_t isOKT=kTRUE;
1369 0 : TObjArray tracks0(ncorr+1);
1370 0 : TObjArray tracks1(ncorr+1);
1371 : //
1372 0 : Double_t dEdx0Tot=seed0->CookdEdxAnalytical(0.02,0.6,kTRUE);
1373 0 : Double_t dEdx1Tot=seed1->CookdEdxAnalytical(0.02,0.6,kTRUE);
1374 0 : Double_t dEdx0Max=seed0->CookdEdxAnalytical(0.02,0.6,kFALSE);
1375 0 : Double_t dEdx1Max=seed1->CookdEdxAnalytical(0.02,0.6,kFALSE);
1376 : //if (TMath::Abs((dEdx0Max+1)/(dEdx0Tot+1)-1.)>0.1) isOK=kFALSE;
1377 : //if (TMath::Abs((dEdx1Max+1)/(dEdx1Tot+1)-1.)>0.1) isOK=kFALSE;
1378 0 : ncl0=0; ncl1=0;
1379 0 : for (Int_t icorr=-1; icorr<ncorr; icorr++){
1380 0 : AliExternalTrackParam rtrack0=btrack0;
1381 0 : AliExternalTrackParam rtrack1=btrack1;
1382 : AliTPCCorrection *corr = 0;
1383 0 : if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
1384 : //
1385 0 : for (Int_t irow=kMaxRow; irow--;){
1386 0 : AliTPCclusterMI *cluster=seed0->GetClusterPointer(irow);
1387 0 : if (!cluster) continue;
1388 0 : if (!isOKT) break;
1389 0 : Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1390 0 : transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD); // transform to global
1391 0 : Float_t r[3]={static_cast<Float_t>(rD[0]),static_cast<Float_t>(rD[1]),static_cast<Float_t>(rD[2])};
1392 0 : if (corr){
1393 0 : corr->DistortPoint(r, cluster->GetDetector());
1394 : }
1395 : //
1396 0 : Double_t cov[3]={0.01,0.,0.01};
1397 0 : Double_t lx =cos*r[0]+sin*r[1];
1398 0 : Double_t ly =-sin*r[0]+cos*r[1];
1399 0 : rD[1]=ly; rD[0]=lx; rD[2]=r[2]; //transform to track local
1400 0 : if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, lx,mass,1.,kFALSE)) isOKT=kFALSE;;
1401 0 : if (!rtrack0.Update(&rD[1],cov)) isOKT =kFALSE;
1402 0 : if (icorr<0) ncl0++;
1403 0 : }
1404 : //
1405 0 : for (Int_t irow=kMaxRow; irow--;){
1406 0 : AliTPCclusterMI *cluster=seed1->GetClusterPointer(irow);
1407 0 : if (!cluster) continue;
1408 0 : if (!isOKT) break;
1409 0 : Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1410 0 : transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD);
1411 0 : Float_t r[3]={static_cast<Float_t>(rD[0]),static_cast<Float_t>(rD[1]),static_cast<Float_t>(rD[2])};
1412 0 : if (corr){
1413 0 : corr->DistortPoint(r, cluster->GetDetector());
1414 : }
1415 : //
1416 0 : Double_t cov[3]={0.01,0.,0.01};
1417 0 : Double_t lx =cos*r[0]+sin*r[1];
1418 0 : Double_t ly =-sin*r[0]+cos*r[1];
1419 0 : rD[1]=ly; rD[0]=lx; rD[2]=r[2];
1420 0 : if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, lx,mass,1.,kFALSE)) isOKT=kFALSE;
1421 0 : if (!rtrack1.Update(&rD[1],cov)) isOKT=kFALSE;
1422 0 : if (icorr<0) ncl1++;
1423 0 : }
1424 0 : if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,10.,kFALSE)) isOKT=kFALSE;
1425 0 : if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,10.,kFALSE)) isOKT=kFALSE;
1426 0 : if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,1.,kFALSE)) isOKT=kFALSE;
1427 0 : if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,1.,kFALSE)) isOKT=kFALSE;
1428 0 : const Double_t *param0=rtrack0.GetParameter();
1429 0 : const Double_t *param1=rtrack1.GetParameter();
1430 0 : for (Int_t ipar=0; ipar<4; ipar++){
1431 0 : if (TMath::Abs(param1[ipar]-param0[ipar])>kMaxDelta[ipar]) isOK=kFALSE;
1432 : }
1433 0 : tracks0.AddAt(rtrack0.Clone(), icorr+1);
1434 0 : tracks1.AddAt(rtrack1.Clone(), icorr+1);
1435 0 : AliExternalTrackParam out0=*(ftrack0->GetTPCOut());
1436 0 : AliExternalTrackParam out1=*(ftrack1->GetTPCOut());
1437 0 : Int_t nentries=TMath::Min(ncl0,ncl1);
1438 :
1439 0 : if (icorr<0) {
1440 0 : (*pcstream)<<"cosmic"<<
1441 0 : "isOK="<<isOK<< // correct all propagation update and also residuals
1442 0 : "isOKT="<<isOKT<< // correct all propagation update
1443 0 : "isOKTrigger="<<isOKTrigger<< // correct? estimate of trigger offset
1444 0 : "id="<<cosmicType<<
1445 : //
1446 : //
1447 0 : "cross="<<crossCounter<<
1448 0 : "vDT.="<<&vectorDT<<
1449 : //
1450 0 : "dTime="<<deltaTime<< // delta time using the A-c side cross
1451 0 : "dTimeCross="<<deltaTimeCross<< // delta time using missing clusters
1452 : //
1453 0 : "dEdx0Max="<<dEdx0Max<<
1454 0 : "dEdx0Tot="<<dEdx0Tot<<
1455 0 : "dEdx1Max="<<dEdx1Max<<
1456 0 : "dEdx1Tot="<<dEdx1Tot<<
1457 : //
1458 0 : "track0.="<<track0<< // original track 0
1459 0 : "track1.="<<track1<< // original track 1
1460 0 : "out0.="<<&out0<< // outer track 0
1461 0 : "out1.="<<&out1<< // outer track 1
1462 0 : "rtrack0.="<<&rtrack0<< // refitted track with current transform
1463 0 : "rtrack1.="<<&rtrack1<< //
1464 0 : "ncl0="<<ncl0<<
1465 0 : "ncl1="<<ncl1<<
1466 0 : "entries="<<nentries<< // number of clusters
1467 : "\n";
1468 : }
1469 0 : }
1470 : //
1471 :
1472 0 : if (isOK){
1473 0 : Int_t nentries=TMath::Min(ncl0,ncl1);
1474 0 : for (Int_t ipar=0; ipar<5; ipar++){
1475 0 : for (Int_t icorr=-1; icorr<ncorr; icorr++){
1476 : AliTPCCorrection *corr = 0;
1477 0 : if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
1478 : //
1479 0 : AliExternalTrackParam *param0=(AliExternalTrackParam *) tracks0.At(icorr+1);
1480 0 : AliExternalTrackParam *param1=(AliExternalTrackParam *) tracks1.At(icorr+1);
1481 0 : distortions[icorr+1]=param1->GetParameter()[ipar]-param0->GetParameter()[ipar];
1482 0 : if (icorr>=0){
1483 0 : distortions[icorr+1]-=distortions[0];
1484 0 : }
1485 : //
1486 0 : if (icorr<0){
1487 0 : Double_t bz=AliTrackerBase::GetBz();
1488 0 : Double_t gxyz[3];
1489 0 : param0->GetXYZ(gxyz);
1490 0 : Int_t dtype=20;
1491 0 : Double_t theta=param0->GetParameter()[3];
1492 0 : Double_t phi = alpha;
1493 0 : Double_t snp = track0->GetInnerParam()->GetSnp();
1494 0 : Double_t mean= distortions[0];
1495 0 : Int_t index = param0->GetIndex(ipar,ipar);
1496 0 : Double_t rms=TMath::Sqrt(param1->GetCovariance()[index]+param1->GetCovariance()[index]);
1497 0 : if (crossCounter<1) rms*=1;
1498 0 : Double_t sector=9*phi/TMath::Pi();
1499 0 : Double_t dsec = sector-TMath::Nint(sector+0.5);
1500 0 : Double_t gx=gxyz[0],gy=gxyz[1],gz=gxyz[2];
1501 0 : Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
1502 0 : Double_t dRrec=0;
1503 : // Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5;
1504 0 : Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5;
1505 :
1506 0 : (*pcstream)<<"fit"<< // dump valus for fit
1507 0 : "run="<<run<< //run number
1508 0 : "bz="<<bz<< // magnetic filed used
1509 0 : "dtype="<<dtype<< // detector match type 20
1510 0 : "ptype="<<ipar<< // parameter type
1511 0 : "theta="<<theta<< // theta
1512 0 : "phi="<<phi<< // phi
1513 0 : "snp="<<snp<< // snp
1514 0 : "mean="<<mean<< // mean dist value
1515 0 : "rms="<<rms<< // rms
1516 0 : "sector="<<sector<<
1517 0 : "dsec="<<dsec<<
1518 : //
1519 0 : "refX="<<refX<< // reference radius
1520 0 : "gx="<<gx<< // global position
1521 0 : "gy="<<gy<< // global position
1522 0 : "gz="<<gz<< // global position
1523 0 : "dRrec="<<dRrec<< // delta Radius in reconstruction
1524 0 : "pt="<<pt<< //1/pt
1525 0 : "id="<<cosmicType<< //type of the cosmic used
1526 0 : "entries="<<nentries;// number of entries in bin
1527 0 : (*pcstream)<<"cosmicDebug"<<
1528 0 : "p0.="<<param0<< // dump distorted track 0
1529 0 : "p1.="<<param1; // dump distorted track 1
1530 0 : }
1531 0 : if (icorr>=0){
1532 0 : (*pcstream)<<"fit"<<
1533 0 : Form("%s=",corr->GetName())<<distortions[icorr+1]; // dump correction value
1534 0 : (*pcstream)<<"cosmicDebug"<<
1535 0 : Form("%s=",corr->GetName())<<distortions[icorr+1]<< // dump correction value
1536 0 : Form("%sp0.=",corr->GetName())<<param0<< // dump distorted track 0
1537 0 : Form("%sp1.=",corr->GetName())<<param1; // dump distorted track 1
1538 : }
1539 : } //loop corrections
1540 0 : (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1541 0 : (*pcstream)<<"cosmicDebug"<<"isOK="<<isOK<<"\n";
1542 : } //loop over parameters
1543 0 : } // dump results
1544 0 : }//loop tracks
1545 0 : delete [] distortions;
1546 0 : }
1547 :
1548 :
1549 :
1550 : 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)
1551 : {
1552 : //
1553 : // Estimate trigger offset between random cosmic event and "physics" trigger
1554 : // Efficiency about 50 % of cases:
1555 : // Cases:
1556 : // 0. Tracks crossing A side and C side - no match in z - 30 % of cases
1557 : // 1. Track only on one side and dissapear at small or at high radiuses - 50 % of cases
1558 : // 1.a) rmax<Rc && tmax<Tcmax && tmax>tmin => deltaT=Tcmax-tmax
1559 : // 1.b) rmin>Rcmin && tmin<Tcmax && tmin>tmax => deltaT=Tcmax-tmin
1560 : // // case the z matching gives proper time
1561 : // 1.c) rmax<Rc && tmax>Tcmin && tmax<tmin => deltaT=-tmax
1562 : //
1563 : // check algorithm:
1564 : // TCut cutStop = "(min(rmax0,rmax1)<235||abs(rmin0-rmin1)>10)"; // tracks not registered for full time
1565 :
1566 : // Combinations:
1567 : // 0-1 - forbidden
1568 : // 0-2 - forbidden
1569 : // 0-3 - occur - wrong correlation
1570 : // 1-2 - occur - wrong correlation
1571 : // 1-3 - forbidden
1572 : // 2-3 - occur - small number of outlyers -20%
1573 : // Frequency:
1574 : // 0 - 106
1575 : // 1 - 265
1576 : // 2 - 206
1577 : // 3 - 367
1578 : //
1579 : const Double_t kMaxRCut=235; // max radius
1580 0 : const Double_t kMinRCut=TMath::Max(dcaR,90.); // min radius
1581 : const Double_t kMaxDCut=30; // max distance for minimal radius
1582 : const Double_t kMinTime=110;
1583 : const Double_t kMaxTime=950;
1584 : Double_t deltaT=0;
1585 : Int_t counter=0;
1586 0 : vectorDT[6]=TMath::Min(TMath::Min(tmin0,tmin1),TMath::Min(tmax0,tmax1));
1587 0 : vectorDT[7]=TMath::Max(TMath::Max(tmin0,tmin1),TMath::Max(tmax0,tmax1));
1588 0 : if (TMath::Min(rmax0,rmax1)<kMaxRCut){
1589 : // max cross - deltaT>0
1590 0 : if (rmax0<kMaxRCut && tmax0 <kMaxTime && tmax0>tmin0) vectorDT[0]=kMaxTime-tmax0; // disapear at CE
1591 0 : if (rmax1<kMaxRCut && tmax1 <kMaxTime && tmax1>tmin1) vectorDT[1]=kMaxTime-tmax1; // disapear at CE
1592 : // min cross - deltaT<0 - OK they are correlated
1593 0 : if (rmax0<kMaxRCut && tmax0 >kMinTime && tmax0<tmin0) vectorDT[2]=-tmax0; // disapear at ROC
1594 0 : if (rmax1<kMaxRCut && tmax1 >kMinTime && tmax1<tmin1) vectorDT[3]=-tmax1; // disapear at ROC
1595 : }
1596 0 : if (rmin0> kMinRCut+kMaxDCut && tmin0 <kMaxTime && tmin0>tmax0) vectorDT[4]=kMaxTime-tmin0;
1597 0 : if (rmin1> kMinRCut+kMaxDCut && tmin1 <kMaxTime && tmin1>tmax1) vectorDT[5]=kMaxTime-tmin1;
1598 : Bool_t isOK=kTRUE;
1599 0 : for (Int_t i=0; i<6;i++) {
1600 0 : if (TMath::Abs(vectorDT[i])>0) {
1601 0 : counter++;
1602 0 : if (vectorDT[i]+vectorDT[6]<0) isOK=kFALSE;
1603 0 : if (vectorDT[i]+vectorDT[7]>kMaxTime) isOK=kFALSE;
1604 0 : if (isOK) deltaT=vectorDT[i];
1605 : }
1606 : }
1607 0 : return deltaT;
1608 : }
1609 :
|
