]> 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
Eff C++ warning removal
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibCosmic.cxx
CommitLineData
f7f33dec 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
54b76c13 16/*
108953e9 17 Comments to be written here:
54b76c13 18 1. What do we calibrate.
19 2. How to interpret results
20 3. Simple example
21 4. Analysis using debug streamers.
22
23
24
25 3.Simple example
26 // To make cosmic scan the user interaction neccessary
27 //
28 .x ~/UliStyle.C
29 gSystem->Load("libANALYSIS");
30 gSystem->Load("libTPCcalib");
31 TFile fcalib("CalibObjects.root");
32 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
33 AliTPCcalibCosmic * cosmic = ( AliTPCcalibCosmic *)array->FindObject("cosmicTPC");
34
35
36
37*/
38
39
40
f7f33dec 41#include "Riostream.h"
42#include "TChain.h"
43#include "TTree.h"
44#include "TH1F.h"
45#include "TH2F.h"
46#include "TList.h"
47#include "TMath.h"
48#include "TCanvas.h"
49#include "TFile.h"
54b76c13 50#include "TF1.h"
f7f33dec 51
54b76c13 52#include "AliTPCclusterMI.h"
f7f33dec 53#include "AliTPCseed.h"
54#include "AliESDVertex.h"
55#include "AliESDEvent.h"
56#include "AliESDfriend.h"
57#include "AliESDInputHandler.h"
54b76c13 58#include "AliAnalysisManager.h"
f7f33dec 59
60#include "AliTracker.h"
61#include "AliMagFMaps.h"
54b76c13 62#include "AliTPCCalROC.h"
f7f33dec 63
64#include "AliLog.h"
65
66#include "AliTPCcalibCosmic.h"
67
68#include "TTreeStream.h"
69#include "AliTPCTracklet.h"
70
71ClassImp(AliTPCcalibCosmic)
72
73
74AliTPCcalibCosmic::AliTPCcalibCosmic()
75 :AliTPCcalibBase(),
bca20570 76 fGainMap(0),
9b27d39b 77 fHistNTracks(0),
78 fClusters(0),
79 fModules(0),
80 fHistPt(0),
9b27d39b 81 fDeDx(0),
54b76c13 82 fDeDxMIP(0),
83 fMIPvalue(1),
9b27d39b 84 fCutMaxD(5), // maximal distance in rfi ditection
85 fCutTheta(0.03), // maximal distan theta
86 fCutMinDir(-0.99) // direction vector products
f7f33dec 87{
54b76c13 88 AliInfo("Default Constructor");
f7f33dec 89}
90
91
92AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
93 :AliTPCcalibBase(),
bca20570 94 fGainMap(0),
9b27d39b 95 fHistNTracks(0),
96 fClusters(0),
97 fModules(0),
98 fHistPt(0),
9b27d39b 99 fDeDx(0),
54b76c13 100 fDeDxMIP(0),
101 fMIPvalue(1),
9b27d39b 102 fCutMaxD(5), // maximal distance in rfi ditection
103 fCutTheta(0.03), // maximal distan theta
104 fCutMinDir(-0.99) // direction vector products
f7f33dec 105{
106 SetName(name);
107 SetTitle(title);
108 AliMagFMaps * field = new AliMagFMaps("dummy1", "dummy2",0,5,0);
109 AliTracker::SetFieldMap(field, kTRUE);
54b76c13 110
111 fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
112 fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
113 fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
114 fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
115 fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
f7f33dec 116 BinLogX(fDeDx);
54b76c13 117 fDeDxMIP = new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
118
9b27d39b 119 AliInfo("Non Default Constructor");
bca20570 120 //
f7f33dec 121}
122
123AliTPCcalibCosmic::~AliTPCcalibCosmic(){
124 //
125 //
126 //
127}
128
129
9b27d39b 130
131
132
f7f33dec 133void AliTPCcalibCosmic::Process(AliESDEvent *event) {
9b27d39b 134 //
135 //
136 //
f7f33dec 137 if (!event) {
138 Printf("ERROR: ESD not available");
139 return;
9b27d39b 140 }
f7f33dec 141 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
142 if (!ESDfriend) {
143 Printf("ERROR: ESDfriend not available");
144 return;
145 }
2acad464 146
108953e9 147
54b76c13 148 FindPairs(event); // nearly everything takes place in find pairs...
149
2acad464 150 if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
9b27d39b 151 Int_t ntracks=event->GetNumberOfTracks();
152 fHistNTracks->Fill(ntracks);
9b27d39b 153 if (ntracks==0) return;
9b27d39b 154
54b76c13 155}
9b27d39b 156
f7f33dec 157
158
54b76c13 159void AliTPCcalibCosmic::Analyze() {
160
161 fMIPvalue = CalculateMIPvalue(fDeDxMIP);
162
163 return;
164
165}
166
f7f33dec 167
f7f33dec 168
9b27d39b 169void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
170 //
171 // Find cosmic pairs
04087794 172 //
173 // Track0 is choosen in upper TPC part
174 // Track1 is choosen in lower TPC part
9b27d39b 175 //
2acad464 176 if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
9b27d39b 177 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
178 Int_t ntracks=event->GetNumberOfTracks();
9b27d39b 179 TObjArray tpcSeeds(ntracks);
180 if (ntracks==0) return;
181 Double_t vtxx[3]={0,0,0};
182 Double_t svtxx[3]={0.000001,0.000001,100.};
183 AliESDVertex vtx(vtxx,svtxx);
184 //
185 //track loop
186 //
54b76c13 187 for (Int_t i=0;i<ntracks;++i) {
188 AliESDtrack *track = event->GetTrack(i);
189 fClusters->Fill(track->GetTPCNcls());
190
191 const AliExternalTrackParam * trackIn = track->GetInnerParam();
192 const AliExternalTrackParam * trackOut = track->GetOuterParam();
193 if (!trackIn) continue;
194 if (!trackOut) continue;
9b27d39b 195
196 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
197 TObject *calibObject;
198 AliTPCseed *seed = 0;
199 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
200 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
201 }
202 if (seed) tpcSeeds.AddAt(seed,i);
54b76c13 203
204 Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
205 if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
206 fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
207 //
208 if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
209 //
210 if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap)>300) {
211 TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
212 if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
213 if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
214 }
215
216 }
217
9b27d39b 218 }
54b76c13 219
9b27d39b 220 if (ntracks<2) return;
221 //
222 // Find pairs
223 //
224 for (Int_t i=0;i<ntracks;++i) {
225 AliESDtrack *track0 = event->GetTrack(i);
04087794 226 // track0 - choosen upper part
227 if (!track0) continue;
228 if (!track0->GetOuterParam()) continue;
229 if (track0->GetOuterParam()->GetAlpha()<0) continue;
e9362f9d 230 Double_t dir0[3];
231 track0->GetDirection(dir0);
04087794 232 for (Int_t j=0;j<ntracks;++j) {
233 if (i==j) continue;
234 AliESDtrack *track1 = event->GetTrack(j);
235 //track 1 lower part
236 if (!track1) continue;
237 if (!track1->GetOuterParam()) continue;
238 if (track1->GetOuterParam()->GetAlpha()>0) continue;
239 //
e9362f9d 240 Double_t dir1[3];
241 track1->GetDirection(dir1);
54b76c13 242
9b27d39b 243 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
244 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
245 if (! seed0) continue;
246 if (! seed1) continue;
bca20570 247 Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
248 Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
249 //
250 Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
251 Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
252 //
253 Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
254 Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
255 //
e9362f9d 256 Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
9b27d39b 257 Float_t d0 = track0->GetLinearD(0,0);
258 Float_t d1 = track1->GetLinearD(0,0);
259 //
260 // conservative cuts - convergence to be guarantied
261 // applying before track propagation
262 if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0
263 if (dir>fCutMinDir) continue; // direction vector product
264 Float_t bz = AliTracker::GetBz();
265 Float_t dvertex0[2]; //distance to 0,0
266 Float_t dvertex1[2]; //distance to 0,0
267 track0->GetDZ(0,0,0,bz,dvertex0);
268 track1->GetDZ(0,0,0,bz,dvertex1);
269 if (TMath::Abs(dvertex0[1])>250) continue;
270 if (TMath::Abs(dvertex1[1])>250) continue;
271 //
272 //
273 //
274 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
275 AliExternalTrackParam param0(*track0);
276 AliExternalTrackParam param1(*track1);
277 //
278 // Propagate using Magnetic field and correct fo material budget
279 //
280 AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
281 AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
282 //
283 // Propagate rest to the 0,0 DCA - z should be ignored
284 //
285 Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
286 Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
287 //
288 param0.GetDZ(0,0,0,bz,dvertex0);
289 param1.GetDZ(0,0,0,bz,dvertex1);
290 //
291 Double_t xyz0[3];//,pxyz0[3];
292 Double_t xyz1[3];//,pxyz1[3];
293 param0.GetXYZ(xyz0);
294 param1.GetXYZ(xyz1);
295 Bool_t isPair = IsPair(&param0,&param1);
296 //
54b76c13 297 if (isPair) FillAcordeHist(track0);
298 //
9b27d39b 299 if (fStreamLevel>0){
300 TTreeSRedirector * cstream = GetDebugStreamer();
108953e9 301 //printf("My stream=%p\n",(void*)cstream);
e9362f9d 302 AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
303 AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
304 AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
305 AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
306 Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
307 Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
308 Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
309 Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
9b27d39b 310 if (cstream) {
311 (*cstream) << "Track0" <<
108953e9 312 "run="<<fRun<< // run number
313 "event="<<fEvent<< // event number
314 "time="<<fTime<< // time stamp of event
315 "trigger="<<fTrigger<< // trigger
316 "mag="<<fMagF<< // magnetic field
9b27d39b 317 "dir="<<dir<< // direction
54b76c13 318 "OK="<<isPair<< // will be accepted
9b27d39b 319 "b0="<<b0<< // propagate status
320 "b1="<<b1<< // propagate status
e9362f9d 321 "crossI="<<isCrossI<< // cross inner
322 "crossO="<<isCrossO<< // cross outer
323 //
9b27d39b 324 "Orig0.=" << track0 << // original track 0
325 "Orig1.=" << track1 << // original track 1
326 "Tr0.="<<&param0<< // track propagated to the DCA 0,0
327 "Tr1.="<<&param1<< // track propagated to the DCA 0,0
e9362f9d 328 "Ip0.="<<ip0<< // inner param - upper
329 "Ip1.="<<ip1<< // inner param - lower
330 "Op0.="<<op0<< // outer param - upper
331 "Op1.="<<op1<< // outer param - lower
332 //
9b27d39b 333 "v00="<<dvertex0[0]<< // distance using kalman
334 "v01="<<dvertex0[1]<< //
335 "v10="<<dvertex1[0]<< //
336 "v11="<<dvertex1[1]<< //
337 "d0="<<d0<< // linear distance to 0,0
338 "d1="<<d1<< // linear distance to 0,0
339 //
e9362f9d 340 //
341 //
342 "x00="<<xyz0[0]<< // global position close to vertex
9b27d39b 343 "x01="<<xyz0[1]<<
344 "x02="<<xyz0[2]<<
345 //
e9362f9d 346 "x10="<<xyz1[0]<< // global position close to vertex
9b27d39b 347 "x11="<<xyz1[1]<<
348 "x12="<<xyz1[2]<<
349 //
e9362f9d 350 "alpha0="<<alpha0<<
351 "alpha1="<<alpha1<<
352 "dir00="<<dir0[0]<< // direction upper
353 "dir01="<<dir0[1]<<
354 "dir02="<<dir0[2]<<
355 //
356 "dir10="<<dir1[0]<< // direction lower
357 "dir11="<<dir1[1]<<
358 "dir12="<<dir1[2]<<
359 //
360 //
54b76c13 361 "Seed0.=" << seed0 << // original seed 0
362 "Seed1.=" << seed1 << // original seed 1
bca20570 363 //
364 "dedx0="<<dedx0<< // dedx0 - all
365 "dedx1="<<dedx1<< // dedx1 - all
366 //
54b76c13 367 "dedx0I="<<dedx0I<< // dedx0 - inner ROC
368 "dedx1I="<<dedx1I<< // dedx1 - inner ROC
bca20570 369 //
54b76c13 370 "dedx0O="<<dedx0O<< // dedx0 - outer ROC
371 "dedx1O="<<dedx1O<< // dedx1 - outer ROC
9b27d39b 372 "\n";
373 }
374 }
375 }
376 }
377}
378
9b27d39b 379
bca20570 380
381
54b76c13 382void AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
383
384 // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
385 if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
386
387 const Double_t AcordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
388 const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
389
390 Double_t r[3];
391 upperTrack->GetXYZ(r);
392 Double_t d[3];
393 upperTrack->GetDirection(d);
394 Double_t x,z;
395 z = r[2] + (d[2]/d[1])*(AcordePlane - r[1]);
396 x = r[0] + (d[0]/d[1])*(AcordePlane - r[1]);
397
398 if (x > roof) {
399 x = r[0] + (d[0]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
400 z = r[2] + (d[2]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
401 }
402 if (x < -roof) {
403 x = r[0] + (d[0]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
404 z = r[2] + (d[2]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
405 }
406
407 fModules->Fill(z, x);
408
409}
410
411
412
bca20570 413Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
414
415 TIterator* iter = li->MakeIterator();
416 AliTPCcalibCosmic* cal = 0;
417
418 while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
419 if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
420 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
421 return -1;
422 }
423
54b76c13 424 fHistNTracks->Add(cal->GetHistNTracks());
425 fClusters->Add(cal-> GetHistClusters());
426 fModules->Add(cal->GetHistAcorde());
427 fHistPt->Add(cal->GetHistPt());
428 fDeDx->Add(cal->GetHistDeDx());
429 fDeDxMIP->Add(cal->GetHistMIP());
bca20570 430
431 }
432
433 return 0;
9b27d39b 434
f7f33dec 435}
436
54b76c13 437
9b27d39b 438Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
439 //
440 //
441 /*
442 // 0. Same direction - OPOSITE - cutDir +cutT
443 TCut cutDir("cutDir","dir<-0.99")
444 // 1.
445 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
446 //
447 // 2. The same rphi
448 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
449 //
450 //
451 //
452 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
453 // 1/Pt diff cut
454 */
455 const Double_t *p0 = tr0->GetParameter();
456 const Double_t *p1 = tr1->GetParameter();
457 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
458 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
459 Double_t d0[3], d1[3];
460 tr0->GetDirection(d0);
461 tr1->GetDirection(d1);
462 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
463 //
464 return kTRUE;
465}
54b76c13 466
467
468
469Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
470
471 TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
472 funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
473 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
474 hist->Fit(funcDoubleGaus);
475 Double_t MIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
476
477 delete funcDoubleGaus;
478
479 return MIPvalue;
480
481}
9b27d39b 482
483
f7f33dec 484
54b76c13 485
486void AliTPCcalibCosmic::CalculateBetheParams(TH2F *hist, Double_t * initialParam) {
487
488 return;
489
490}
491
492
f7f33dec 493void AliTPCcalibCosmic::BinLogX(TH1 *h) {
494
495 // Method for the correct logarithmic binning of histograms
496
497 TAxis *axis = h->GetXaxis();
498 int bins = axis->GetNbins();
499
500 Double_t from = axis->GetXmin();
501 Double_t to = axis->GetXmax();
502 Double_t *new_bins = new Double_t[bins + 1];
503
504 new_bins[0] = from;
505 Double_t factor = pow(to/from, 1./bins);
506
507 for (int i = 1; i <= bins; i++) {
508 new_bins[i] = factor * new_bins[i-1];
509 }
510 axis->Set(bins, new_bins);
511 delete new_bins;
512
513}
514
428b47be 515AliExternalTrackParam *AliTPCcalibCosmic::Invert(AliExternalTrackParam *input)
516{
517 //
518 // Invert paramerameter - not covariance yet
519 //
520 AliExternalTrackParam *output = new AliExternalTrackParam(*input);
521 Double_t * param = (Double_t*)output->GetParameter();
522 param[0]*=-1;
523 param[3]*=-1;
524 param[4]*=-1;
525 //
526 return output;
527}
7b18d067 528
e9362f9d 529/*
530
531void Init(){
532
533.x ~/UliStyle.C
534.x ~/rootlogon.C
535gSystem->Load("libSTAT.so");
536gSystem->Load("libANALYSIS");
537gSystem->Load("libTPCcalib");
538gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
539
540gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
541AliXRDPROOFtoolkit tool;
542TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
543chain->Lookup();
544
545TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01"); // OK
546TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2"); // OK
108953e9 547TCut cutP1("cutP1","abs(Tr0.fP[1]-Tr1.fP[1])<10"); // OK
e9362f9d 548TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
549TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
550TCut cutA=cutT+cutD+cutPt+cutN+cutP1;
551
552TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
553
554//
555chain->Draw(">>listEL",cutA,"entryList");
556TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
557chain->SetEntryList(elist);
558//
559chain->Draw(">>listV40Z100","abs(d0)<40&&abs(v01)<100","entryList");
560TEntryList *elistV40Z100 = (TEntryList*)gDirectory->Get("listV40Z100");
561chain->SetEntryList(elistV40Z100);
562
563//
564// Aliases
565//
566chain->SetAlias("side","(-1+(Tr0.fP[1]>0)*2)");
567chain->SetAlias("hpt","abs(Tr0.fP[4])<0.2");
108953e9 568chain->SetAlias("signy","(-1+(Tr0.fP[0]>0)*2)");
569
570chain->SetAlias("dy","Tr0.fP[0]+Tr1.fP[0]");
571chain->SetAlias("dz","Tr0.fP[1]-Tr1.fP[1]");
572chain->SetAlias("d1pt","Tr0.fP[4]+Tr1.fP[4]");
e9362f9d 573chain->SetAlias("dtheta","(Tr0.fP[3]+Tr1.fP[3])");
108953e9 574chain->SetAlias("dphi","(Tr0.fAlpha-Tr1.fAlpha-pi)");
575
e9362f9d 576chain->SetAlias("mtheta","(Tr0.fP[3]-Tr1.fP[3])*0.5")
577chain->SetAlias("sa","(sin(Tr0.fAlpha+0.))");
578chain->SetAlias("ca","(cos(Tr0.fAlpha+0.))");
579
108953e9 580
581
582chain->Draw("dy:sqrt(abs(Tr0.fP[4]))>>hisdyA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
583hisdyA->FitSlicesY();
584hisdyA_2->SetXTitle("#sqrt{1/p_{t}}");
585hisdyA_2->SetYTitle("#sigma_{y}(cm)");
586hisdyA_2->SetTitle("Cosmic - Y matching");
587hisdyA_2->SetMaximum(0.5);
588
589
590chain->Draw("dy:sqrt(abs(Tr0.fP[4]))>>hisdyC(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
591hisdyC->FitSlicesY();
592hisdyC_2->SetXTitle("#sqrt{1/p_{t}}");
593hisdyC_2->SetYTitle("#sigma_{y}(cm)");
594hisdyC_2->SetTitle("Cosmic - Y matching");
595hisdyC_2->SetMaximum(1);
596hisdyC_2->SetMinimum(0);
597hisdyC_2->SetMarkerStyle(22);
598hisdyA_2->SetMarkerStyle(21);
599hisdyC_2->SetMarkerSize(1.5);
600hisdzA_2->SetMarkerSize(1.5);
601hisdyC_2->Draw();
602hisdyA_2->Draw("same");
603gPad->SaveAs("~/Calibration/Cosmic/pic/ymatching.gif")
604
605chain->Draw("dz:sqrt(abs(Tr0.fP[4]))>>hisdzA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
606hisdzA->FitSlicesY();
607hisdzA_2->SetXTitle("#sqrt{1/p_{t}}");
608hisdzA_2->SetYTitle("#sigma_{z}(cm)");
609hisdzA_2->SetTitle("Cosmic - Z matching - A side ");
610hisdzA_2->SetMaximum(0.5);
611
612chain->Draw("dz:sqrt(abs(Tr0.fP[4]))>>hisdzC(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
613hisdzC->FitSlicesY();
614hisdzC_2->SetXTitle("#sqrt{1/p_{t}}");
615hisdzC_2->SetYTitle("#sigma_{z}(cm)");
616hisdzC_2->SetTitle("Cosmic - Z matching");
617hisdzC_2->SetMaximum(0.5);
618hisdzC_2->SetMarkerStyle(22);
619hisdzA_2->SetMarkerStyle(21);
620hisdzC_2->SetMarkerSize(1.5);
621hisdzA_2->SetMarkerSize(1.5);
622
623hisdzC_2->Draw();
624hisdzA_2->Draw("same");
625
626
627//
628// PICTURE 1/pt
629//
630chain->Draw("d1pt:sqrt(abs(Tr0.fP[4]))>>hisd1ptA(5,0,1,30,-0.1,0.1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
631hisd1ptA->FitSlicesY();
632hisd1ptA_2->SetXTitle("#sqrt{1/p_{t}}");
633hisd1ptA_2->SetYTitle("#sigma_{z}(cm)");
634hisd1ptA_2->SetTitle("Cosmic - Z matching - A side ");
635hisd1ptA_2->SetMaximum(0.5);
636
637chain->Draw("d1pt:sqrt(abs(Tr0.fP[4]))>>hisd1ptC(5,0,1,30,-0.1,0.1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
638hisd1ptC->FitSlicesY();
639hisd1ptC_2->SetXTitle("#sqrt{1/p_{t}}");
640hisd1ptC_2->SetYTitle("#sigma_{1/pt}(1/GeV)");
641hisd1ptC_2->SetTitle("Cosmic - 1/pt matching");
642hisd1ptC_2->SetMaximum(0.05);
643hisd1ptC_2->SetMarkerStyle(22);
644hisd1ptA_2->SetMarkerStyle(21);
645hisd1ptC_2->SetMarkerSize(1.5);
646hisd1ptA_2->SetMarkerSize(1.5);
647
648hisd1ptC_2->Draw();
649hisd1ptA_2->Draw("same");
650gPad->SaveAs("~/Calibration/Cosmic/pic/1ptmatching.gif")
651
652//
653// Theta
654//
655chain->Draw("dtheta:sqrt(abs(Tr0.fP[4]))>>hisdthetaA(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
656hisdthetaA->FitSlicesY();
657hisdthetaA_2->SetXTitle("#sqrt{1/p_{t}}");
658hisdthetaA_2->SetYTitle("#sigma_{#theta}(cm)");
659hisdthetaA_2->SetTitle("Cosmic - Z matching - A side ");
660hisdthetaA_2->SetMaximum(0.5);
661
662chain->Draw("dtheta:sqrt(abs(Tr0.fP[4]))>>hisdthetaC(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
663hisdthetaC->FitSlicesY();
664hisdthetaC_2->SetXTitle("#sqrt{1/p_{t}}");
665hisdthetaC_2->SetYTitle("#sigma_{#theta}(rad)");
666hisdthetaC_2->SetTitle("Cosmic - Theta matching");
667hisdthetaC_2->SetMaximum(0.01);
668hisdthetaC_2->SetMinimum(0.0);
669hisdthetaC_2->SetMarkerStyle(22);
670hisdthetaA_2->SetMarkerStyle(21);
671hisdthetaC_2->SetMarkerSize(1.5);
672hisdthetaA_2->SetMarkerSize(1.5);
673
674hisdthetaC_2->Draw();
675hisdthetaA_2->Draw("same");
676gPad->SaveAs("~/Calibration/Cosmic/pic/thetamatching.gif")
677//
678// Phi
679//
680chain->Draw("dphi:sqrt(abs(Tr0.fP[4]))>>hisdphiA(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
681hisdphiA->FitSlicesY();
682hisdphiA_2->SetXTitle("#sqrt{1/p_{t}}");
683hisdphiA_2->SetYTitle("#sigma_{#phi}(rad)");
684hisdphiA_2->SetTitle("Cosmic - Z matching - A side ");
685hisdphiA_2->SetMaximum(0.5);
686
687chain->Draw("dphi:sqrt(abs(Tr0.fP[4]))>>hisdphiC(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
688hisdphiC->FitSlicesY();
689hisdphiC_2->SetXTitle("#sqrt{1/p_{t}}");
690hisdphiC_2->SetYTitle("#sigma_{#phi}(rad)");
691hisdphiC_2->SetTitle("Cosmic - Phi matching");
692hisdphiC_2->SetMaximum(0.01);
693hisdphiC_2->SetMinimum(0.0);
694hisdphiC_2->SetMarkerStyle(22);
695hisdphiA_2->SetMarkerStyle(21);
696hisdphiC_2->SetMarkerSize(1.5);
697hisdphiA_2->SetMarkerSize(1.5);
698
699hisdphiC_2->Draw();
700hisdphiA_2->Draw("same");
701gPad->SaveAs("~/Calibration/Cosmic/pic/phimatching.gif")
702
703
704
e9362f9d 705}
706
707
708*/
709
710
711/*
712void MatchTheta(){
713
714TStatToolkit toolkit;
715Double_t chi2=0;
716Int_t npoints=0;
717TVectorD fitParamA0;
718TVectorD fitParamA1;
719TVectorD fitParamC0;
720TVectorD fitParamC1;
721TMatrixD covMatrix;
722
723
724TString fstring="";
725//
108953e9 726fstring+="mtheta++";
e9362f9d 727fstring+="ca++";
728fstring+="sa++";
729fstring+="ca*mtheta++";
730fstring+="sa*mtheta++";
108953e9 731//
e9362f9d 732fstring+="side++";
108953e9 733fstring+="side*mtheta++";
e9362f9d 734fstring+="side*ca++";
735fstring+="side*sa++";
736fstring+="side*ca*mtheta++";
737fstring+="side*sa*mtheta++";
738
739
108953e9 740TString *strTheta0 = toolkit.FitPlane(chain,"dtheta",fstring->Data(), "hpt&&!crossI&&!crossO", chi2,npoints,fitParamA0,covMatrix,0.8);
e9362f9d 741chain->SetAlias("dtheta0",strTheta0.Data());
742strTheta0->Tokenize("+")->Print();
743
744
745//fstring+="mtheta++";
746//fstring+="mtheta^2++";
747//fstring+="ca*mtheta^2++";
748//fstring+="sa*mtheta^2++";
749
750
751
752}
753
754*/
755
756
757
758
759/*
760 void PosCorrection()
761
762
763
764
765 TStatToolkit toolkit;
766 Double_t chi2=0;
767 Int_t npoints=0;
768 TVectorD fitParam;
769 TMatrixD covMatrix;
770
771 //Theta
772chain->SetAlias("dthe","(Tr0.fP[3]+Tr1.fP[3])");
773chain->SetAlias("sign","(-1+(Tr0.fP[1]>0)*2)");
774chain->SetAlias("di","(1-abs(Tr0.fP[1])/250)");
775//
776//
777TString strFit="";
778//
779strFit+="sign++"; //1
780strFit+="Tr0.fP[3]++"; //2
781//
782strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]>0)++"; //3
783strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]<0)++"; //4
784//
785strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]>0)++"; //5
786strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]<0)++"; //6
787//
788strFit+="sin(Tr0.fAlpha)*Tr0.fP[3]++"; //7
789strFit+="cos(Tr0.fAlpha)*Tr0.fP[3]++"; //8
790
791
792 //
793 TString * thetaParam = toolkit.FitPlane(chain,"dthe", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix,0.8,0,10000)
794
795 chain->SetAlias("corTheta",thetaParam->Data());
796 chain->Draw("dthe:Tr0.fP[1]","","",50000);
797
798
799
800*/
801
802
bca20570 803
54b76c13 804/*
bca20570 805
7b18d067 806void AliTPCcalibCosmic::dEdxCorrection(){
e9362f9d 807 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01"); // OK
808 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2"); // OK
809 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
810 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
811 TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
812 TCut cutA=cutT+cutD+cutPt+cutN+cutS;
7b18d067 813
5ab97d2b 814
815 .x ~/UliStyle.C
816 gSystem->Load("libANALYSIS");
817 gSystem->Load("libTPCcalib");
818 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
819 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
820 AliXRDPROOFtoolkit tool;
bca20570 821 TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
822 chain->Lookup();
e9362f9d 823
824 chain->Draw(">>listEL",cutA,"entryList");
825 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
826 chain->SetEntryList(elist);
7b18d067 827
828 .x ~/rootlogon.C
829 gSystem->Load("libSTAT.so");
e9362f9d 830 TStatToolkit toolkit;
7b18d067 831 Double_t chi2=0;
832 Int_t npoints=0;
833 TVectorD fitParam;
834 TMatrixD covMatrix;
835
836 chain->Draw("Tr0.fP[4]+Tr1.fP[4]","OK"+cutA);
837
838 TString strFit;
839 strFit+="(Tr0.fP[1]/250)++";
840 strFit+="(Tr0.fP[1]/250)^2++";
841 strFit+="(Tr0.fP[3])++";
842 strFit+="(Tr0.fP[3])^2++";
843
e9362f9d 844 TString * ptParam = TStatToolkit::FitPlane(chain,"Tr0.fP[4]+Tr1.fP[4]", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix)
bca20570 845
846
bca20570 847
7b18d067 848*/
849
850
428b47be 851/*
852gSystem->Load("libANALYSIS");
853gSystem->Load("libSTAT");
854gSystem->Load("libTPCcalib");
855
856TStatToolkit toolkit;
857Double_t chi2;
858TVectorD fitParam;
859TMatrixD covMatrix;
860Int_t npoints;
861//
862TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03"); // OK
863TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5"); // OK
864TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.2&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
865TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>110");
108953e9 866TCut cutA=cutT+cutD+cutPt+cutN;
428b47be 867
868
869
870TTree * chain = Track0;
871
872
873chain->SetAlias("norm","signalTot0.fElements[3]/signalTot1.fElements[3]");
874//
875chain->SetAlias("dr1","(signalTot0.fElements[1]/signalTot0.fElements[3])");
876chain->SetAlias("dr2","(signalTot0.fElements[2]/signalTot0.fElements[3])");
877chain->SetAlias("dr4","(signalTot0.fElements[4]/signalTot0.fElements[3])");
878chain->SetAlias("dr5","(signalTot0.fElements[5]/signalTot0.fElements[3])");
879
880TString fstring="";
881fstring+="dr1++";
882fstring+="dr2++";
883fstring+="dr4++";
884fstring+="dr5++";
885//
886fstring+="dr1*dr2++";
887fstring+="dr1*dr4++";
888fstring+="dr1*dr5++";
889fstring+="dr2*dr4++";
890fstring+="dr2*dr5++";
891fstring+="dr4*dr5++";
892
893
894
895TString *strqdedx = toolkit.FitPlane(chain,"norm",fstring->Data(), cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
896
897chain->SetAlias("corQT",strqdedx->Data());
898
899*/
900
901
108953e9 902/*
903 chain->SetProof(kTRUE);
904 chain->Draw("Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE):Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)",""+cutA,"",1000);
905
906
907chain->Draw("Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE)/Seed1.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE)>>his(100,0.5,1.5)","min(Orig0.fTPCncls,Orig1.fTPCncls)>130"+cutA,"",50000);
908
909*/
910
911