]> 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
Use debug stream only if requested
[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/*
17 Comments to be written here:
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
54b76c13 147 FindPairs(event); // nearly everything takes place in find pairs...
148
2acad464 149 if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
9b27d39b 150 Int_t ntracks=event->GetNumberOfTracks();
151 fHistNTracks->Fill(ntracks);
9b27d39b 152 if (ntracks==0) return;
9b27d39b 153
54b76c13 154}
9b27d39b 155
f7f33dec 156
157
54b76c13 158void AliTPCcalibCosmic::Analyze() {
159
160 fMIPvalue = CalculateMIPvalue(fDeDxMIP);
161
162 return;
163
164}
165
f7f33dec 166
f7f33dec 167
9b27d39b 168void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
169 //
170 // Find cosmic pairs
04087794 171 //
172 // Track0 is choosen in upper TPC part
173 // Track1 is choosen in lower TPC part
9b27d39b 174 //
2acad464 175 if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
9b27d39b 176 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
177 Int_t ntracks=event->GetNumberOfTracks();
9b27d39b 178 TObjArray tpcSeeds(ntracks);
179 if (ntracks==0) return;
180 Double_t vtxx[3]={0,0,0};
181 Double_t svtxx[3]={0.000001,0.000001,100.};
182 AliESDVertex vtx(vtxx,svtxx);
183 //
184 //track loop
185 //
54b76c13 186 for (Int_t i=0;i<ntracks;++i) {
187 AliESDtrack *track = event->GetTrack(i);
188 fClusters->Fill(track->GetTPCNcls());
189
190 const AliExternalTrackParam * trackIn = track->GetInnerParam();
191 const AliExternalTrackParam * trackOut = track->GetOuterParam();
192 if (!trackIn) continue;
193 if (!trackOut) continue;
9b27d39b 194
195 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
196 TObject *calibObject;
197 AliTPCseed *seed = 0;
198 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
199 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
200 }
201 if (seed) tpcSeeds.AddAt(seed,i);
54b76c13 202
203 Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
204 if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
205 fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
206 //
207 if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
208 //
209 if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap)>300) {
210 TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
211 if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
212 if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
213 }
214
215 }
216
9b27d39b 217 }
54b76c13 218
9b27d39b 219 if (ntracks<2) return;
220 //
221 // Find pairs
222 //
223 for (Int_t i=0;i<ntracks;++i) {
224 AliESDtrack *track0 = event->GetTrack(i);
04087794 225 // track0 - choosen upper part
226 if (!track0) continue;
227 if (!track0->GetOuterParam()) continue;
228 if (track0->GetOuterParam()->GetAlpha()<0) continue;
e9362f9d 229 Double_t dir0[3];
230 track0->GetDirection(dir0);
04087794 231 for (Int_t j=0;j<ntracks;++j) {
232 if (i==j) continue;
233 AliESDtrack *track1 = event->GetTrack(j);
234 //track 1 lower part
235 if (!track1) continue;
236 if (!track1->GetOuterParam()) continue;
237 if (track1->GetOuterParam()->GetAlpha()>0) continue;
238 //
e9362f9d 239 Double_t dir1[3];
240 track1->GetDirection(dir1);
54b76c13 241
9b27d39b 242 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
243 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
244 if (! seed0) continue;
245 if (! seed1) continue;
bca20570 246 Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
247 Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
248 //
249 Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
250 Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
251 //
252 Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
253 Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
254 //
e9362f9d 255 Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
9b27d39b 256 Float_t d0 = track0->GetLinearD(0,0);
257 Float_t d1 = track1->GetLinearD(0,0);
258 //
259 // conservative cuts - convergence to be guarantied
260 // applying before track propagation
261 if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0
262 if (dir>fCutMinDir) continue; // direction vector product
263 Float_t bz = AliTracker::GetBz();
264 Float_t dvertex0[2]; //distance to 0,0
265 Float_t dvertex1[2]; //distance to 0,0
266 track0->GetDZ(0,0,0,bz,dvertex0);
267 track1->GetDZ(0,0,0,bz,dvertex1);
268 if (TMath::Abs(dvertex0[1])>250) continue;
269 if (TMath::Abs(dvertex1[1])>250) continue;
270 //
271 //
272 //
273 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
274 AliExternalTrackParam param0(*track0);
275 AliExternalTrackParam param1(*track1);
276 //
277 // Propagate using Magnetic field and correct fo material budget
278 //
279 AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
280 AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
281 //
282 // Propagate rest to the 0,0 DCA - z should be ignored
283 //
284 Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
285 Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
286 //
287 param0.GetDZ(0,0,0,bz,dvertex0);
288 param1.GetDZ(0,0,0,bz,dvertex1);
289 //
290 Double_t xyz0[3];//,pxyz0[3];
291 Double_t xyz1[3];//,pxyz1[3];
292 param0.GetXYZ(xyz0);
293 param1.GetXYZ(xyz1);
294 Bool_t isPair = IsPair(&param0,&param1);
295 //
54b76c13 296 if (isPair) FillAcordeHist(track0);
297 //
9b27d39b 298 if (fStreamLevel>0){
299 TTreeSRedirector * cstream = GetDebugStreamer();
300 printf("My stream=%p\n",(void*)cstream);
e9362f9d 301 AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
302 AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
303 AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
304 AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
305 Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
306 Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
307 Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
308 Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
309 Int_t runNr = event->GetRunNumber();
310 UInt_t time = event->GetTimeStamp();
9b27d39b 311 if (cstream) {
312 (*cstream) << "Track0" <<
e9362f9d 313 "runNr="<<runNr<< // run number
314 "time="<<time<< // time stamp of event
9b27d39b 315 "dir="<<dir<< // direction
54b76c13 316 "OK="<<isPair<< // will be accepted
9b27d39b 317 "b0="<<b0<< // propagate status
318 "b1="<<b1<< // propagate status
e9362f9d 319 "crossI="<<isCrossI<< // cross inner
320 "crossO="<<isCrossO<< // cross outer
321 //
9b27d39b 322 "Orig0.=" << track0 << // original track 0
323 "Orig1.=" << track1 << // original track 1
324 "Tr0.="<<&param0<< // track propagated to the DCA 0,0
325 "Tr1.="<<&param1<< // track propagated to the DCA 0,0
e9362f9d 326 "Ip0.="<<ip0<< // inner param - upper
327 "Ip1.="<<ip1<< // inner param - lower
328 "Op0.="<<op0<< // outer param - upper
329 "Op1.="<<op1<< // outer param - lower
330 //
9b27d39b 331 "v00="<<dvertex0[0]<< // distance using kalman
332 "v01="<<dvertex0[1]<< //
333 "v10="<<dvertex1[0]<< //
334 "v11="<<dvertex1[1]<< //
335 "d0="<<d0<< // linear distance to 0,0
336 "d1="<<d1<< // linear distance to 0,0
337 //
e9362f9d 338 //
339 //
340 "x00="<<xyz0[0]<< // global position close to vertex
9b27d39b 341 "x01="<<xyz0[1]<<
342 "x02="<<xyz0[2]<<
343 //
e9362f9d 344 "x10="<<xyz1[0]<< // global position close to vertex
9b27d39b 345 "x11="<<xyz1[1]<<
346 "x12="<<xyz1[2]<<
347 //
e9362f9d 348 "alpha0="<<alpha0<<
349 "alpha1="<<alpha1<<
350 "dir00="<<dir0[0]<< // direction upper
351 "dir01="<<dir0[1]<<
352 "dir02="<<dir0[2]<<
353 //
354 "dir10="<<dir1[0]<< // direction lower
355 "dir11="<<dir1[1]<<
356 "dir12="<<dir1[2]<<
357 //
358 //
54b76c13 359 "Seed0.=" << seed0 << // original seed 0
360 "Seed1.=" << seed1 << // original seed 1
bca20570 361 //
362 "dedx0="<<dedx0<< // dedx0 - all
363 "dedx1="<<dedx1<< // dedx1 - all
364 //
54b76c13 365 "dedx0I="<<dedx0I<< // dedx0 - inner ROC
366 "dedx1I="<<dedx1I<< // dedx1 - inner ROC
bca20570 367 //
54b76c13 368 "dedx0O="<<dedx0O<< // dedx0 - outer ROC
369 "dedx1O="<<dedx1O<< // dedx1 - outer ROC
9b27d39b 370 "\n";
371 }
372 }
373 }
374 }
375}
376
9b27d39b 377
bca20570 378
379
54b76c13 380void AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
381
382 // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
383 if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
384
385 const Double_t AcordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
386 const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
387
388 Double_t r[3];
389 upperTrack->GetXYZ(r);
390 Double_t d[3];
391 upperTrack->GetDirection(d);
392 Double_t x,z;
393 z = r[2] + (d[2]/d[1])*(AcordePlane - r[1]);
394 x = r[0] + (d[0]/d[1])*(AcordePlane - r[1]);
395
396 if (x > roof) {
397 x = r[0] + (d[0]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
398 z = r[2] + (d[2]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
399 }
400 if (x < -roof) {
401 x = r[0] + (d[0]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
402 z = r[2] + (d[2]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
403 }
404
405 fModules->Fill(z, x);
406
407}
408
409
410
bca20570 411Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
412
413 TIterator* iter = li->MakeIterator();
414 AliTPCcalibCosmic* cal = 0;
415
416 while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
417 if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
418 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
419 return -1;
420 }
421
54b76c13 422 fHistNTracks->Add(cal->GetHistNTracks());
423 fClusters->Add(cal-> GetHistClusters());
424 fModules->Add(cal->GetHistAcorde());
425 fHistPt->Add(cal->GetHistPt());
426 fDeDx->Add(cal->GetHistDeDx());
427 fDeDxMIP->Add(cal->GetHistMIP());
bca20570 428
429 }
430
431 return 0;
9b27d39b 432
f7f33dec 433}
434
54b76c13 435
9b27d39b 436Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
437 //
438 //
439 /*
440 // 0. Same direction - OPOSITE - cutDir +cutT
441 TCut cutDir("cutDir","dir<-0.99")
442 // 1.
443 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
444 //
445 // 2. The same rphi
446 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
447 //
448 //
449 //
450 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
451 // 1/Pt diff cut
452 */
453 const Double_t *p0 = tr0->GetParameter();
454 const Double_t *p1 = tr1->GetParameter();
455 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
456 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
457 Double_t d0[3], d1[3];
458 tr0->GetDirection(d0);
459 tr1->GetDirection(d1);
460 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
461 //
462 return kTRUE;
463}
54b76c13 464
465
466
467Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
468
469 TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
470 funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
471 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
472 hist->Fit(funcDoubleGaus);
473 Double_t MIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
474
475 delete funcDoubleGaus;
476
477 return MIPvalue;
478
479}
9b27d39b 480
481
f7f33dec 482
54b76c13 483
484void AliTPCcalibCosmic::CalculateBetheParams(TH2F *hist, Double_t * initialParam) {
485
486 return;
487
488}
489
490
f7f33dec 491void AliTPCcalibCosmic::BinLogX(TH1 *h) {
492
493 // Method for the correct logarithmic binning of histograms
494
495 TAxis *axis = h->GetXaxis();
496 int bins = axis->GetNbins();
497
498 Double_t from = axis->GetXmin();
499 Double_t to = axis->GetXmax();
500 Double_t *new_bins = new Double_t[bins + 1];
501
502 new_bins[0] = from;
503 Double_t factor = pow(to/from, 1./bins);
504
505 for (int i = 1; i <= bins; i++) {
506 new_bins[i] = factor * new_bins[i-1];
507 }
508 axis->Set(bins, new_bins);
509 delete new_bins;
510
511}
512
428b47be 513AliExternalTrackParam *AliTPCcalibCosmic::Invert(AliExternalTrackParam *input)
514{
515 //
516 // Invert paramerameter - not covariance yet
517 //
518 AliExternalTrackParam *output = new AliExternalTrackParam(*input);
519 Double_t * param = (Double_t*)output->GetParameter();
520 param[0]*=-1;
521 param[3]*=-1;
522 param[4]*=-1;
523 //
524 return output;
525}
7b18d067 526
e9362f9d 527/*
528
529void Init(){
530
531.x ~/UliStyle.C
532.x ~/rootlogon.C
533gSystem->Load("libSTAT.so");
534gSystem->Load("libANALYSIS");
535gSystem->Load("libTPCcalib");
536gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
537
538gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
539AliXRDPROOFtoolkit tool;
540TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
541chain->Lookup();
542
543TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01"); // OK
544TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2"); // OK
545TCut cutP1("cutP1","abs(Tr0.fP[1]-Tr1.fP[1])<4"); // OK
546TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
547TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
548TCut cutA=cutT+cutD+cutPt+cutN+cutP1;
549
550TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
551
552//
553chain->Draw(">>listEL",cutA,"entryList");
554TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
555chain->SetEntryList(elist);
556//
557chain->Draw(">>listV40Z100","abs(d0)<40&&abs(v01)<100","entryList");
558TEntryList *elistV40Z100 = (TEntryList*)gDirectory->Get("listV40Z100");
559chain->SetEntryList(elistV40Z100);
560
561//
562// Aliases
563//
564chain->SetAlias("side","(-1+(Tr0.fP[1]>0)*2)");
565chain->SetAlias("hpt","abs(Tr0.fP[4])<0.2");
566chain->SetAlias("dtheta","(Tr0.fP[3]+Tr1.fP[3])");
567chain->SetAlias("mtheta","(Tr0.fP[3]-Tr1.fP[3])*0.5")
568chain->SetAlias("sa","(sin(Tr0.fAlpha+0.))");
569chain->SetAlias("ca","(cos(Tr0.fAlpha+0.))");
570
571}
572
573
574*/
575
576
577/*
578void MatchTheta(){
579
580TStatToolkit toolkit;
581Double_t chi2=0;
582Int_t npoints=0;
583TVectorD fitParamA0;
584TVectorD fitParamA1;
585TVectorD fitParamC0;
586TVectorD fitParamC1;
587TMatrixD covMatrix;
588
589
590TString fstring="";
591//
592fstring+="ca++";
593fstring+="sa++";
594fstring+="ca*mtheta++";
595fstring+="sa*mtheta++";
596fstring+="side++";
597fstring+="side*ca++";
598fstring+="side*sa++";
599fstring+="side*ca*mtheta++";
600fstring+="side*sa*mtheta++";
601
602
603TString *strTheta0 = toolkit.FitPlane(chain,"dtheta",fstring->Data(), "hpt"+cutS, chi2,npoints,fitParamA0,covMatrix,0.8);
604chain->SetAlias("dtheta0",strTheta0.Data());
605strTheta0->Tokenize("+")->Print();
606
607
608//fstring+="mtheta++";
609//fstring+="mtheta^2++";
610//fstring+="ca*mtheta^2++";
611//fstring+="sa*mtheta^2++";
612
613
614
615}
616
617*/
618
619
620
621
622/*
623 void PosCorrection()
624
625
626
627
628 TStatToolkit toolkit;
629 Double_t chi2=0;
630 Int_t npoints=0;
631 TVectorD fitParam;
632 TMatrixD covMatrix;
633
634 //Theta
635chain->SetAlias("dthe","(Tr0.fP[3]+Tr1.fP[3])");
636chain->SetAlias("sign","(-1+(Tr0.fP[1]>0)*2)");
637chain->SetAlias("di","(1-abs(Tr0.fP[1])/250)");
638//
639//
640TString strFit="";
641//
642strFit+="sign++"; //1
643strFit+="Tr0.fP[3]++"; //2
644//
645strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]>0)++"; //3
646strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]<0)++"; //4
647//
648strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]>0)++"; //5
649strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]<0)++"; //6
650//
651strFit+="sin(Tr0.fAlpha)*Tr0.fP[3]++"; //7
652strFit+="cos(Tr0.fAlpha)*Tr0.fP[3]++"; //8
653
654
655 //
656 TString * thetaParam = toolkit.FitPlane(chain,"dthe", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix,0.8,0,10000)
657
658 chain->SetAlias("corTheta",thetaParam->Data());
659 chain->Draw("dthe:Tr0.fP[1]","","",50000);
660
661
662
663*/
664
665
bca20570 666
54b76c13 667/*
bca20570 668
7b18d067 669void AliTPCcalibCosmic::dEdxCorrection(){
e9362f9d 670 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01"); // OK
671 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2"); // OK
672 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
673 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
674 TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
675 TCut cutA=cutT+cutD+cutPt+cutN+cutS;
7b18d067 676
5ab97d2b 677
678 .x ~/UliStyle.C
679 gSystem->Load("libANALYSIS");
680 gSystem->Load("libTPCcalib");
681 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
682 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
683 AliXRDPROOFtoolkit tool;
bca20570 684 TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
685 chain->Lookup();
e9362f9d 686
687 chain->Draw(">>listEL",cutA,"entryList");
688 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
689 chain->SetEntryList(elist);
7b18d067 690
691 .x ~/rootlogon.C
692 gSystem->Load("libSTAT.so");
e9362f9d 693 TStatToolkit toolkit;
7b18d067 694 Double_t chi2=0;
695 Int_t npoints=0;
696 TVectorD fitParam;
697 TMatrixD covMatrix;
698
699 chain->Draw("Tr0.fP[4]+Tr1.fP[4]","OK"+cutA);
700
701 TString strFit;
702 strFit+="(Tr0.fP[1]/250)++";
703 strFit+="(Tr0.fP[1]/250)^2++";
704 strFit+="(Tr0.fP[3])++";
705 strFit+="(Tr0.fP[3])^2++";
706
e9362f9d 707 TString * ptParam = TStatToolkit::FitPlane(chain,"Tr0.fP[4]+Tr1.fP[4]", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix)
bca20570 708
709
bca20570 710
7b18d067 711*/
712
713
428b47be 714/*
715gSystem->Load("libANALYSIS");
716gSystem->Load("libSTAT");
717gSystem->Load("libTPCcalib");
718
719TStatToolkit toolkit;
720Double_t chi2;
721TVectorD fitParam;
722TMatrixD covMatrix;
723Int_t npoints;
724//
725TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03"); // OK
726TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5"); // OK
727TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.2&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
728TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>110");
729TCut cutA="abs(norm-1)<0.3"+cutT+cutD+cutPt+cutN;
730
731
732
733TTree * chain = Track0;
734
735
736chain->SetAlias("norm","signalTot0.fElements[3]/signalTot1.fElements[3]");
737//
738chain->SetAlias("dr1","(signalTot0.fElements[1]/signalTot0.fElements[3])");
739chain->SetAlias("dr2","(signalTot0.fElements[2]/signalTot0.fElements[3])");
740chain->SetAlias("dr4","(signalTot0.fElements[4]/signalTot0.fElements[3])");
741chain->SetAlias("dr5","(signalTot0.fElements[5]/signalTot0.fElements[3])");
742
743TString fstring="";
744fstring+="dr1++";
745fstring+="dr2++";
746fstring+="dr4++";
747fstring+="dr5++";
748//
749fstring+="dr1*dr2++";
750fstring+="dr1*dr4++";
751fstring+="dr1*dr5++";
752fstring+="dr2*dr4++";
753fstring+="dr2*dr5++";
754fstring+="dr4*dr5++";
755
756
757
758TString *strqdedx = toolkit.FitPlane(chain,"norm",fstring->Data(), cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
759
760chain->SetAlias("corQT",strqdedx->Data());
761
762*/
763
764