Adding simple example to load default debug streamer
[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;
9b27d39b 229 Double_t d1[3];
230 track0->GetDirection(d1);
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 //
239 Double_t d2[3];
240 track1->GetDirection(d2);
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 //
9b27d39b 255 Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
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);
301 if (cstream) {
302 (*cstream) << "Track0" <<
303 "dir="<<dir<< // direction
54b76c13 304 "OK="<<isPair<< // will be accepted
9b27d39b 305 "b0="<<b0<< // propagate status
306 "b1="<<b1<< // propagate status
307 "Orig0.=" << track0 << // original track 0
308 "Orig1.=" << track1 << // original track 1
309 "Tr0.="<<&param0<< // track propagated to the DCA 0,0
310 "Tr1.="<<&param1<< // track propagated to the DCA 0,0
311 "v00="<<dvertex0[0]<< // distance using kalman
312 "v01="<<dvertex0[1]<< //
313 "v10="<<dvertex1[0]<< //
314 "v11="<<dvertex1[1]<< //
315 "d0="<<d0<< // linear distance to 0,0
316 "d1="<<d1<< // linear distance to 0,0
317 //
318 "x00="<<xyz0[0]<<
319 "x01="<<xyz0[1]<<
320 "x02="<<xyz0[2]<<
321 //
322 "x10="<<xyz1[0]<<
323 "x11="<<xyz1[1]<<
324 "x12="<<xyz1[2]<<
325 //
54b76c13 326 "Seed0.=" << seed0 << // original seed 0
327 "Seed1.=" << seed1 << // original seed 1
bca20570 328 //
329 "dedx0="<<dedx0<< // dedx0 - all
330 "dedx1="<<dedx1<< // dedx1 - all
331 //
54b76c13 332 "dedx0I="<<dedx0I<< // dedx0 - inner ROC
333 "dedx1I="<<dedx1I<< // dedx1 - inner ROC
bca20570 334 //
54b76c13 335 "dedx0O="<<dedx0O<< // dedx0 - outer ROC
336 "dedx1O="<<dedx1O<< // dedx1 - outer ROC
9b27d39b 337 "\n";
338 }
339 }
340 }
341 }
342}
343
9b27d39b 344
bca20570 345
346
54b76c13 347void AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
348
349 // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
350 if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
351
352 const Double_t AcordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
353 const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
354
355 Double_t r[3];
356 upperTrack->GetXYZ(r);
357 Double_t d[3];
358 upperTrack->GetDirection(d);
359 Double_t x,z;
360 z = r[2] + (d[2]/d[1])*(AcordePlane - r[1]);
361 x = r[0] + (d[0]/d[1])*(AcordePlane - r[1]);
362
363 if (x > roof) {
364 x = r[0] + (d[0]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
365 z = r[2] + (d[2]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
366 }
367 if (x < -roof) {
368 x = r[0] + (d[0]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
369 z = r[2] + (d[2]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
370 }
371
372 fModules->Fill(z, x);
373
374}
375
376
377
bca20570 378Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
379
380 TIterator* iter = li->MakeIterator();
381 AliTPCcalibCosmic* cal = 0;
382
383 while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
384 if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
385 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
386 return -1;
387 }
388
54b76c13 389 fHistNTracks->Add(cal->GetHistNTracks());
390 fClusters->Add(cal-> GetHistClusters());
391 fModules->Add(cal->GetHistAcorde());
392 fHistPt->Add(cal->GetHistPt());
393 fDeDx->Add(cal->GetHistDeDx());
394 fDeDxMIP->Add(cal->GetHistMIP());
bca20570 395
396 }
397
398 return 0;
9b27d39b 399
f7f33dec 400}
401
54b76c13 402
9b27d39b 403Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
404 //
405 //
406 /*
407 // 0. Same direction - OPOSITE - cutDir +cutT
408 TCut cutDir("cutDir","dir<-0.99")
409 // 1.
410 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
411 //
412 // 2. The same rphi
413 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
414 //
415 //
416 //
417 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
418 // 1/Pt diff cut
419 */
420 const Double_t *p0 = tr0->GetParameter();
421 const Double_t *p1 = tr1->GetParameter();
422 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
423 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
424 Double_t d0[3], d1[3];
425 tr0->GetDirection(d0);
426 tr1->GetDirection(d1);
427 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
428 //
429 return kTRUE;
430}
54b76c13 431
432
433
434Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
435
436 TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
437 funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
438 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
439 hist->Fit(funcDoubleGaus);
440 Double_t MIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
441
442 delete funcDoubleGaus;
443
444 return MIPvalue;
445
446}
9b27d39b 447
448
f7f33dec 449
54b76c13 450
451void AliTPCcalibCosmic::CalculateBetheParams(TH2F *hist, Double_t * initialParam) {
452
453 return;
454
455}
456
457
f7f33dec 458void AliTPCcalibCosmic::BinLogX(TH1 *h) {
459
460 // Method for the correct logarithmic binning of histograms
461
462 TAxis *axis = h->GetXaxis();
463 int bins = axis->GetNbins();
464
465 Double_t from = axis->GetXmin();
466 Double_t to = axis->GetXmax();
467 Double_t *new_bins = new Double_t[bins + 1];
468
469 new_bins[0] = from;
470 Double_t factor = pow(to/from, 1./bins);
471
472 for (int i = 1; i <= bins; i++) {
473 new_bins[i] = factor * new_bins[i-1];
474 }
475 axis->Set(bins, new_bins);
476 delete new_bins;
477
478}
479
428b47be 480AliExternalTrackParam *AliTPCcalibCosmic::Invert(AliExternalTrackParam *input)
481{
482 //
483 // Invert paramerameter - not covariance yet
484 //
485 AliExternalTrackParam *output = new AliExternalTrackParam(*input);
486 Double_t * param = (Double_t*)output->GetParameter();
487 param[0]*=-1;
488 param[3]*=-1;
489 param[4]*=-1;
490 //
491 return output;
492}
7b18d067 493
bca20570 494
54b76c13 495/*
bca20570 496
7b18d067 497void AliTPCcalibCosmic::dEdxCorrection(){
2acad464 498 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03"); // OK
499 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5"); // OK
7b18d067 500 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.2&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
2acad464 501 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>70");
7b18d067 502 TCut cutA=cutT+cutD+cutPt+cutN;
503
5ab97d2b 504
505 .x ~/UliStyle.C
506 gSystem->Load("libANALYSIS");
507 gSystem->Load("libTPCcalib");
508 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
509 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
510 AliXRDPROOFtoolkit tool;
bca20570 511 TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
512 chain->Lookup();
7b18d067 513
514 .x ~/rootlogon.C
515 gSystem->Load("libSTAT.so");
516
517 Double_t chi2=0;
518 Int_t npoints=0;
519 TVectorD fitParam;
520 TMatrixD covMatrix;
521
522 chain->Draw("Tr0.fP[4]+Tr1.fP[4]","OK"+cutA);
523
524 TString strFit;
525 strFit+="(Tr0.fP[1]/250)++";
526 strFit+="(Tr0.fP[1]/250)^2++";
527 strFit+="(Tr0.fP[3])++";
528 strFit+="(Tr0.fP[3])^2++";
529
530 TString * ptParam = TStatToolkit::FitPlane(chain,"Tr0.fP[4]+Tr1.fP[4]", strFit.Data(),cutA, chi2,npoints,fitParam,covMatrix)
531
532strFit+="(Tr0.fP[1]/250)++";
533strFit+="(Tr0.fP[1]/250)^2++";
534strFit+="(Tr0.fP[3])++";
535strFit+="(Tr0.fP[3])^2++";
536strFit+="(Tr0.fP[1]/250)^2*Tr0.fP[3]++";
537strFit+="(Tr0.fP[1]/250)^2*Tr0.fP[3]^2++";
538//
539
540strFit+="sign(Tr0.fP[1])++"
541strFit+="sign(Tr0.fP[1])*(1-abs(Tr0.fP[1]/250))"
542
543TString * thetaParam = TStatToolkit::FitPlane(chain,"Tr0.fP[3]+Tr1.fP[3]", strFit.Data(),cutA, chi2,npoints,fitParam,covMatrix)
bca20570 544
545
546
547(-0.009263+(Tr0.fP[1]/250)*(-0.009693)+(Tr0.fP[1]/250)^2*(0.009062)+(Tr0.fP[3])*(0.002256)+(Tr0.fP[3])^2*(-0.052775)+(Tr0.fP[1]/250)^2*Tr0.fP[3]*(-0.020627)+(Tr0.fP[1]/250)^2*Tr0.fP[3]^2*(0.049030))
548
7b18d067 549*/
550
551
428b47be 552/*
553gSystem->Load("libANALYSIS");
554gSystem->Load("libSTAT");
555gSystem->Load("libTPCcalib");
556
557TStatToolkit toolkit;
558Double_t chi2;
559TVectorD fitParam;
560TMatrixD covMatrix;
561Int_t npoints;
562//
563TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03"); // OK
564TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5"); // OK
565TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.2&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
566TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>110");
567TCut cutA="abs(norm-1)<0.3"+cutT+cutD+cutPt+cutN;
568
569
570
571TTree * chain = Track0;
572
573
574chain->SetAlias("norm","signalTot0.fElements[3]/signalTot1.fElements[3]");
575//
576chain->SetAlias("dr1","(signalTot0.fElements[1]/signalTot0.fElements[3])");
577chain->SetAlias("dr2","(signalTot0.fElements[2]/signalTot0.fElements[3])");
578chain->SetAlias("dr4","(signalTot0.fElements[4]/signalTot0.fElements[3])");
579chain->SetAlias("dr5","(signalTot0.fElements[5]/signalTot0.fElements[3])");
580
581TString fstring="";
582fstring+="dr1++";
583fstring+="dr2++";
584fstring+="dr4++";
585fstring+="dr5++";
586//
587fstring+="dr1*dr2++";
588fstring+="dr1*dr4++";
589fstring+="dr1*dr5++";
590fstring+="dr2*dr4++";
591fstring+="dr2*dr5++";
592fstring+="dr4*dr5++";
593
594
595
596TString *strqdedx = toolkit.FitPlane(chain,"norm",fstring->Data(), cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
597
598chain->SetAlias("corQT",strqdedx->Data());
599
600*/
601
602
603
604
605
606