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