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