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