1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 Comments to be written here:
18 1. What do we calibrate.
19 2. How to interpret results
21 4. Analysis using debug streamers.
26 // To make cosmic scan the user interaction neccessary
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");
41 #include "Riostream.h"
52 #include "AliTPCclusterMI.h"
53 #include "AliTPCseed.h"
54 #include "AliESDVertex.h"
55 #include "AliESDEvent.h"
56 #include "AliESDfriend.h"
57 #include "AliESDInputHandler.h"
58 #include "AliAnalysisManager.h"
60 #include "AliTracker.h"
62 #include "AliTPCCalROC.h"
66 #include "AliTPCcalibCosmic.h"
67 #include "AliExternalComparison.h"
68 #include "TTreeStream.h"
69 #include "AliTPCTracklet.h"
71 ClassImp(AliTPCcalibCosmic)
74 AliTPCcalibCosmic::AliTPCcalibCosmic()
84 fCutMaxD(5), // maximal distance in rfi ditection
85 fCutMaxDz(40), // maximal distance in z ditection
86 fCutTheta(0.03), // maximal distan theta
87 fCutMinDir(-0.99) // direction vector products
89 AliInfo("Default Constructor");
93 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
103 fCutMaxD(5), // maximal distance in rfi ditection
104 fCutMaxDz(40), // maximal distance in z ditection
105 fCutTheta(0.03), // maximal distan theta
106 fCutMinDir(-0.99) // direction vector products
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);
117 fDeDxMIP = new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
119 AliInfo("Non Default Constructor");
123 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
133 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
138 Printf("ERROR: ESD not available");
141 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
143 Printf("ERROR: ESDfriend not available");
148 FindPairs(event); // nearly everything takes place in find pairs...
150 if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
151 Int_t ntracks=event->GetNumberOfTracks();
152 fHistNTracks->Fill(ntracks);
153 if (ntracks==0) return;
159 void AliTPCcalibCosmic::Analyze() {
161 fMIPvalue = CalculateMIPvalue(fDeDxMIP);
169 void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
173 // Track0 is choosen in upper TPC part
174 // Track1 is choosen in lower TPC part
176 if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
177 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
178 Int_t ntracks=event->GetNumberOfTracks();
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);
187 for (Int_t i=0;i<ntracks;++i) {
188 AliESDtrack *track = event->GetTrack(i);
189 fClusters->Fill(track->GetTPCNcls());
191 const AliExternalTrackParam * trackIn = track->GetInnerParam();
192 const AliExternalTrackParam * trackOut = track->GetOuterParam();
193 if (!trackIn) continue;
194 if (!trackOut) continue;
195 if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // filter laser
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;
204 if (seed) tpcSeeds.AddAt(seed,i);
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));
210 if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
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;
222 if (ntracks<2) return;
226 for (Int_t i=0;i<ntracks;++i) {
227 AliESDtrack *track0 = event->GetTrack(i);
228 // track0 - choosen upper part
229 if (!track0) continue;
230 if (!track0->GetOuterParam()) continue;
231 if (track0->GetOuterParam()->GetAlpha()<0) continue;
233 track0->GetDirection(dir0);
234 for (Int_t j=0;j<ntracks;++j) {
236 AliESDtrack *track1 = event->GetTrack(j);
238 if (!track1) continue;
239 if (!track1->GetOuterParam()) continue;
240 if (track1->GetOuterParam()->GetAlpha()>0) continue;
243 track1->GetDirection(dir1);
245 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
246 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
247 if (! seed0) continue;
248 if (! seed1) continue;
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);
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);
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);
258 Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
259 Float_t d0 = track0->GetLinearD(0,0);
260 Float_t d1 = track1->GetLinearD(0,0);
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;
276 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
277 AliExternalTrackParam param0(*track0);
278 AliExternalTrackParam param1(*track1);
280 // Propagate using Magnetic field and correct fo material budget
282 AliTracker::PropagateTrackTo(¶m0,dmax+1,0.0005,3,kTRUE);
283 AliTracker::PropagateTrackTo(¶m1,dmax+1,0.0005,3,kTRUE);
285 // Propagate rest to the 0,0 DCA - z should be ignored
287 Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
288 Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
290 param0.GetDZ(0,0,0,bz,dvertex0);
291 param1.GetDZ(0,0,0,bz,dvertex1);
292 if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
294 Double_t xyz0[3];//,pxyz0[3];
295 Double_t xyz1[3];//,pxyz1[3];
298 Bool_t isPair = IsPair(¶m0,¶m1);
300 if (isPair) FillAcordeHist(track0);
303 TTreeSRedirector * cstream = GetDebugStreamer();
304 //printf("My stream=%p\n",(void*)cstream);
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]);
314 (*cstream) << "Track0" <<
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
320 "dir="<<dir<< // direction
321 "OK="<<isPair<< // will be accepted
322 "b0="<<b0<< // propagate status
323 "b1="<<b1<< // propagate status
324 "crossI="<<isCrossI<< // cross inner
325 "crossO="<<isCrossO<< // cross outer
327 "Orig0.=" << track0 << // original track 0
328 "Orig1.=" << track1 << // original track 1
329 "Tr0.="<<¶m0<< // track propagated to the DCA 0,0
330 "Tr1.="<<¶m1<< // track propagated to the DCA 0,0
331 "Ip0.="<<ip0<< // inner param - upper
332 "Ip1.="<<ip1<< // inner param - lower
333 "Op0.="<<op0<< // outer param - upper
334 "Op1.="<<op1<< // outer param - lower
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
345 "x00="<<xyz0[0]<< // global position close to vertex
349 "x10="<<xyz1[0]<< // global position close to vertex
355 "dir00="<<dir0[0]<< // direction upper
359 "dir10="<<dir1[0]<< // direction lower
364 "Seed0.=" << seed0 << // original seed 0
365 "Seed1.=" << seed1 << // original seed 1
367 "dedx0="<<dedx0<< // dedx0 - all
368 "dedx1="<<dedx1<< // dedx1 - all
370 "dedx0I="<<dedx0I<< // dedx0 - inner ROC
371 "dedx1I="<<dedx1I<< // dedx1 - inner ROC
373 "dedx0O="<<dedx0O<< // dedx0 - outer ROC
374 "dedx1O="<<dedx1O<< // dedx1 - outer ROC
385 void AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
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;
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
394 upperTrack->GetXYZ(r);
396 upperTrack->GetDirection(d);
398 z = r[2] + (d[2]/d[1])*(AcordePlane - r[1]);
399 x = r[0] + (d[0]/d[1])*(AcordePlane - r[1]);
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]);
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]);
410 fModules->Fill(z, x);
416 Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
418 TIterator* iter = li->MakeIterator();
419 AliTPCcalibCosmic* cal = 0;
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());
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());
441 Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
445 // 0. Same direction - OPOSITE - cutDir +cutT
446 TCut cutDir("cutDir","dir<-0.99")
448 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
451 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
455 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
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;
461 if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
462 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
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;
474 Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
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));
482 delete funcDoubleGaus;
491 void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) {
493 // Not implemented yet
500 void AliTPCcalibCosmic::BinLogX(TH1 *h) {
502 // Method for the correct logarithmic binning of histograms
504 TAxis *axis = h->GetXaxis();
505 int bins = axis->GetNbins();
507 Double_t from = axis->GetXmin();
508 Double_t to = axis->GetXmax();
509 Double_t *new_bins = new Double_t[bins + 1];
512 Double_t factor = pow(to/from, 1./bins);
514 for (int i = 1; i <= bins; i++) {
515 new_bins[i] = factor * new_bins[i-1];
517 axis->Set(bins, new_bins);
522 AliExternalTrackParam *AliTPCcalibCosmic::Invert(AliExternalTrackParam *input)
525 // Invert paramerameter - not covariance yet
527 AliExternalTrackParam *output = new AliExternalTrackParam(*input);
528 Double_t * param = (Double_t*)output->GetParameter();
536 AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
540 AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
541 par1R->Rotate(track0->GetAlpha());
544 Double_t * param = (Double_t*)par1R->GetParameter();
545 Double_t * covar = (Double_t*)par1R->GetCovariance();
552 covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
553 //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
555 par1R->PropagateTo(track0->GetX(),0); // bz shold be set -
557 // printf("Print param\n");
564 void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
566 // Update track 1 with track 2
570 TMatrixD vecXk(5,1); // X vector
571 TMatrixD covXk(5,5); // X covariance
572 TMatrixD matHk(5,5); // vector to mesurement
573 TMatrixD measR(5,5); // measurement error
574 TMatrixD vecZk(5,1); // measurement
576 TMatrixD vecYk(5,1); // Innovation or measurement residual
577 TMatrixD matHkT(5,5);
578 TMatrixD matSk(5,5); // Innovation (or residual) covariance
579 TMatrixD matKk(5,5); // Optimal Kalman gain
580 TMatrixD mat1(5,5); // update covariance matrix
581 TMatrixD covXk2(5,5); //
582 TMatrixD covOut(5,5);
584 Double_t *param1=(Double_t*) track1.GetParameter();
585 Double_t *covar1=(Double_t*) track1.GetCovariance();
586 Double_t *param2=(Double_t*) track2.GetParameter();
587 Double_t *covar2=(Double_t*) track2.GetCovariance();
589 // copy data to the matrix
590 for (Int_t ipar=0; ipar<5; ipar++){
591 vecXk(ipar,0) = param1[ipar];
592 vecZk(ipar,0) = param2[ipar];
593 for (Int_t jpar=0; jpar<5; jpar++){
594 covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
595 measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
602 matHk(0,0)=1; matHk(1,1)= 1; matHk(2,2)= 1;
603 matHk(3,3)= 1; matHk(4,4)= 1; // vector to measurement
605 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
606 matHkT=matHk.T(); matHk.T();
607 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
609 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
610 vecXk += matKk*vecYk; // updated vector
611 mat1(0,0)=1; mat1(1,1)=1; mat1(2,2)=1; mat1(3,3)=1; mat1(4,4)=1;
612 covXk2 = (mat1-(matKk*matHk));
613 covOut = covXk2*covXk;
617 // copy from matrix to parameters
630 for (Int_t ipar=0; ipar<5; ipar++){
631 param1[ipar]= vecXk(ipar,0) ;
632 for (Int_t jpar=0; jpar<5; jpar++){
633 covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
638 void AliTPCcalibCosmic::ProcessTree(TTree * chainTracklet, AliExternalComparison *comp){
640 // Process the debug streamer tree
641 // Possible to modify selection criteria
643 TTreeSRedirector * cstream = new TTreeSRedirector("cosmicdump.root");
644 //AliTPCcalibCosmic *cosmic = this;
646 AliExternalTrackParam * tr0 = 0;
647 AliExternalTrackParam * tr1 = 0;
650 Int_t entries=chainTracklet->GetEntries();
651 for (Int_t i=0; i< entries; i++){
652 chainTracklet->GetBranch("Tr0.")->SetAddress(&tr0);
653 chainTracklet->GetBranch("Tr1.")->SetAddress(&tr1);
654 chainTracklet->GetEntry(i);
657 if (tr0->GetY()==0) continue;
658 if (tr1->GetY()==0) continue;
660 AliExternalTrackParam par0(*tr0);
661 AliExternalTrackParam par1(*tr1);
662 AliExternalTrackParam par1R(*tr1);
663 par1R.Rotate(par1.GetAlpha()+TMath::Pi());
664 AliExternalTrackParam *par1T = MakeTrack(tr0,tr1);
666 printf("%d\t%d\t\n",i, npoints);
670 AliExternalTrackParam par0U=par0;
671 AliExternalTrackParam par1U=*par1T;
673 UpdateTrack(par0U,*par1T);
674 UpdateTrack(par1U,par0);
677 if (i%100==0) printf("%d\t%d\tt\n",i, npoints);
678 Bool_t accept = comp->AcceptPair(&par0,par1T);
680 if (1||fStreamLevel>0){
681 (*cstream)<<"Tracklet"<<
683 "tr0.="<<&par0<< //original track up
684 "tr1.="<<&par1<< //original track down
685 "tr1R.="<<&par1R<< //track1 rotated to 0 frame
686 "tr1T.="<<par1T<< //track1 transformed to the track 0 frame
688 "tr0U.="<<&par0U<< //track 0 updated with track 1
689 "tr1U.="<<&par1U<< //track 1 updated with track 0
695 if (comp) comp->Process(&par0,par1T);
715 gSystem->Load("libSTAT.so");
716 gSystem->Load("libANALYSIS");
717 gSystem->Load("libTPCcalib");
718 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
720 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
721 AliXRDPROOFtoolkit tool;
722 TChain * chainCosmic = tool.MakeChain("cosmic.txt","Track0",0,1000000);
723 chainCosmic->Lookup();
725 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01"); // OK
726 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2"); // OK
727 TCut cutP1("cutP1","abs(Tr0.fP[1]-Tr1.fP[1])<20"); // OK
728 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
729 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>50");
730 TCut cutM("cutM","abs(mag)>0.01");
731 TCut cutA=cutT+cutD+cutPt+cutN+cutP1+"trigger!=16";
733 TCut cuthpt("abs(Tr0.fP[4])+abs(Tr1.fP[4])<0.2");
734 TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
737 chainCosmic->Draw(">>listELP",cutA,"entryList");
738 //TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
739 //TEntryList *elist = (TEntryList*)gProof->GetOutputList()->At(1);
740 chainCosmic->SetEntryList(elist);
742 chainCosmic->Draw(">>listV40Z100","abs(d0)<40&&abs(v01)<100","entryList");
743 TEntryList *elistV40Z100 = (TEntryList*)gDirectory->Get("listV40Z100");
744 chainCosmic->SetEntryList(elistV40Z100);
749 chainCosmic->SetAlias("side","(-1+(Tr0.fP[1]>0)*2)");
750 chainCosmic->SetAlias("hpt","abs(Tr0.fP[4])<0.2");
751 chainCosmic->SetAlias("signy","(-1+(Tr0.fP[0]>0)*2)");
753 chainCosmic->SetAlias("dy","Tr0.fP[0]+Tr1.fP[0]");
754 chainCosmic->SetAlias("dz","Tr0.fP[1]-Tr1.fP[1]");
755 chainCosmic->SetAlias("d1pt","Tr0.fP[4]+Tr1.fP[4]");
756 chainCosmic->SetAlias("dtheta","(Tr0.fP[3]+Tr1.fP[3])");
757 chainCosmic->SetAlias("dphi","(Tr0.fAlpha-Tr1.fAlpha-pi)");
759 chainCosmic->SetAlias("mtheta","(Tr0.fP[3]-Tr1.fP[3])*0.5")
760 chainCosmic->SetAlias("sa","(sin(Tr0.fAlpha+0.))");
761 chainCosmic->SetAlias("ca","(cos(Tr0.fAlpha+0.))");
765 chainCosmic->Draw("dy:sqrt(abs(Tr0.fP[4]))>>hisdyA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
766 hisdyA->FitSlicesY();
767 hisdyA_2->SetXTitle("#sqrt{1/p_{t}}");
768 hisdyA_2->SetYTitle("#sigma_{y}(cm)");
769 hisdyA_2->SetTitle("Cosmic - Y matching");
770 hisdyA_2->SetMaximum(0.5);
773 chainCosmic->Draw("dy:sqrt(abs(Tr0.fP[4]))>>hisdyC(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
774 hisdyC->FitSlicesY();
775 hisdyC_2->SetXTitle("#sqrt{1/p_{t}}");
776 hisdyC_2->SetYTitle("#sigma_{y}(cm)");
777 hisdyC_2->SetTitle("Cosmic - Y matching");
778 hisdyC_2->SetMaximum(1);
779 hisdyC_2->SetMinimum(0);
780 hisdyC_2->SetMarkerStyle(22);
781 hisdyA_2->SetMarkerStyle(21);
782 hisdyC_2->SetMarkerSize(1.5);
783 hisdyA_2->SetMarkerSize(1.5);
785 hisdyA_2->Draw("same");
786 gPad->SaveAs("~/Calibration/Cosmic/pic/ymatching.gif")
788 chainCosmic->Draw("dz:sqrt(abs(Tr0.fP[4]))>>hisdzA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
789 hisdzA->FitSlicesY();
790 hisdzA_2->SetXTitle("#sqrt{1/p_{t}}");
791 hisdzA_2->SetYTitle("#sigma_{z}(cm)");
792 hisdzA_2->SetTitle("Cosmic - Z matching - A side ");
793 hisdzA_2->SetMaximum(0.5);
795 chainCosmic->Draw("dz:sqrt(abs(Tr0.fP[4]))>>hisdzC(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
796 hisdzC->FitSlicesY();
797 hisdzC_2->SetXTitle("#sqrt{1/p_{t}}");
798 hisdzC_2->SetYTitle("#sigma_{z}(cm)");
799 hisdzC_2->SetTitle("Cosmic - Z matching");
800 hisdzC_2->SetMaximum(0.5);
801 hisdzC_2->SetMarkerStyle(22);
802 hisdzA_2->SetMarkerStyle(21);
803 hisdzC_2->SetMarkerSize(1.5);
804 hisdzA_2->SetMarkerSize(1.5);
807 hisdzA_2->Draw("same");
813 chainCosmic->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"+cutM);
814 hisd1ptA->FitSlicesY();
815 hisd1ptA_2->SetXTitle("#sqrt{1/p_{t}}");
816 hisd1ptA_2->SetYTitle("#sigma_{z}(cm)");
817 hisd1ptA_2->SetTitle("Cosmic - Z matching - A side ");
818 hisd1ptA_2->SetMaximum(0.5);
820 chainCosmic->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"+cutM);
821 hisd1ptC->FitSlicesY();
822 hisd1ptC_2->SetXTitle("#sqrt{1/p_{t}}");
823 hisd1ptC_2->SetYTitle("#sigma_{1/pt}(1/GeV)");
824 hisd1ptC_2->SetTitle("Cosmic - 1/pt matching");
825 hisd1ptC_2->SetMaximum(0.05);
826 hisd1ptC_2->SetMarkerStyle(22);
827 hisd1ptA_2->SetMarkerStyle(21);
828 hisd1ptC_2->SetMarkerSize(1.5);
829 hisd1ptA_2->SetMarkerSize(1.5);
832 hisd1ptA_2->Draw("same");
833 gPad->SaveAs("~/Calibration/Cosmic/pic/1ptmatching.gif")
838 chainCosmic->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");
839 hisdthetaA->FitSlicesY();
840 hisdthetaA_2->SetXTitle("#sqrt{1/p_{t}}");
841 hisdthetaA_2->SetYTitle("#sigma_{#theta}(cm)");
842 hisdthetaA_2->SetTitle("Cosmic - Z matching - A side ");
843 hisdthetaA_2->SetMaximum(0.5);
845 chainCosmic->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");
846 hisdthetaC->FitSlicesY();
847 hisdthetaC_2->SetXTitle("#sqrt{1/p_{t}}");
848 hisdthetaC_2->SetYTitle("#sigma_{#theta}(rad)");
849 hisdthetaC_2->SetTitle("Cosmic - Theta matching");
850 hisdthetaC_2->SetMaximum(0.01);
851 hisdthetaC_2->SetMinimum(0.0);
852 hisdthetaC_2->SetMarkerStyle(22);
853 hisdthetaA_2->SetMarkerStyle(21);
854 hisdthetaC_2->SetMarkerSize(1.5);
855 hisdthetaA_2->SetMarkerSize(1.5);
857 hisdthetaC_2->Draw();
858 hisdthetaA_2->Draw("same");
859 gPad->SaveAs("~/Calibration/Cosmic/pic/thetamatching.gif")
863 chainCosmic->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");
864 hisdphiA->FitSlicesY();
865 hisdphiA_2->SetXTitle("#sqrt{1/p_{t}}");
866 hisdphiA_2->SetYTitle("#sigma_{#phi}(rad)");
867 hisdphiA_2->SetTitle("Cosmic - Z matching - A side ");
868 hisdphiA_2->SetMaximum(0.5);
870 chainCosmic->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");
871 hisdphiC->FitSlicesY();
872 hisdphiC_2->SetXTitle("#sqrt{1/p_{t}}");
873 hisdphiC_2->SetYTitle("#sigma_{#phi}(rad)");
874 hisdphiC_2->SetTitle("Cosmic - Phi matching");
875 hisdphiC_2->SetMaximum(0.01);
876 hisdphiC_2->SetMinimum(0.0);
877 hisdphiC_2->SetMarkerStyle(22);
878 hisdphiA_2->SetMarkerStyle(21);
879 hisdphiC_2->SetMarkerSize(1.5);
880 hisdphiA_2->SetMarkerSize(1.5);
883 hisdphiA_2->Draw("same");
884 gPad->SaveAs("~/Calibration/Cosmic/pic/phimatching.gif")
897 TStatToolkit toolkit;
912 fstring+="ca*mtheta++";
913 fstring+="sa*mtheta++";
916 fstring+="side*mtheta++";
917 fstring+="side*ca++";
918 fstring+="side*sa++";
919 fstring+="side*ca*mtheta++";
920 fstring+="side*sa*mtheta++";
923 TString *strTheta0 = toolkit.FitPlane(chain,"dtheta",fstring->Data(), "hpt&&!crossI&&!crossO", chi2,npoints,fitParamA0,covMatrix,0.8);
924 chainCosmic->SetAlias("dtheta0",strTheta0.Data());
925 strTheta0->Tokenize("+")->Print();
928 //fstring+="mtheta++";
929 //fstring+="mtheta^2++";
930 //fstring+="ca*mtheta^2++";
931 //fstring+="sa*mtheta^2++";
948 TStatToolkit toolkit;
955 chainCosmic->SetAlias("dthe","(Tr0.fP[3]+Tr1.fP[3])");
956 chainCosmic->SetAlias("sign","(-1+(Tr0.fP[1]>0)*2)");
957 chainCosmic->SetAlias("di","(1-abs(Tr0.fP[1])/250)");
962 strFit+="sign++"; //1
963 strFit+="Tr0.fP[3]++"; //2
965 strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]>0)++"; //3
966 strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]<0)++"; //4
968 strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]>0)++"; //5
969 strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]<0)++"; //6
971 strFit+="sin(Tr0.fAlpha)*Tr0.fP[3]++"; //7
972 strFit+="cos(Tr0.fAlpha)*Tr0.fP[3]++"; //8
976 TString * thetaParam = toolkit.FitPlane(chain,"dthe", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix,0.8,0,10000)
978 chainCosmic->SetAlias("corTheta",thetaParam->Data());
979 chainCosmic->Draw("dthe:Tr0.fP[1]","","",50000);
989 void AliTPCcalibCosmic::dEdxCorrection(){
990 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01"); // OK
991 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2"); // OK
992 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
993 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
994 TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
995 TCut cutA=cutT+cutD+cutPt+cutN+cutS;
999 gSystem->Load("libANALYSIS");
1000 gSystem->Load("libTPCcalib");
1001 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
1002 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
1003 AliXRDPROOFtoolkit tool;
1004 TChain * chainCosmic = tool.MakeChain("cosmic.txt","Track0",0,1000000);
1005 chainCosmic->Lookup();
1007 chainCosmic->Draw(">>listEL",cutA,"entryList");
1008 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
1009 chainCosmic->SetEntryList(elist);
1012 gSystem->Load("libSTAT.so");
1013 TStatToolkit toolkit;
1019 chainCosmic->Draw("Tr0.fP[4]+Tr1.fP[4]","OK"+cutA);
1022 strFit+="(Tr0.fP[1]/250)++";
1023 strFit+="(Tr0.fP[1]/250)^2++";
1024 strFit+="(Tr0.fP[3])++";
1025 strFit+="(Tr0.fP[3])^2++";
1027 TString * ptParam = TStatToolkit::FitPlane(chain,"Tr0.fP[4]+Tr1.fP[4]", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix)
1035 gSystem->Load("libANALYSIS");
1036 gSystem->Load("libSTAT");
1037 gSystem->Load("libTPCcalib");
1039 TStatToolkit toolkit;
1045 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03"); // OK
1046 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5"); // OK
1047 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.2&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
1048 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>110");
1049 TCut cutA=cutT+cutD+cutPt+cutN;
1053 TTree * chainCosmic = Track0;
1056 chainCosmic->SetAlias("norm","signalTot0.fElements[3]/signalTot1.fElements[3]");
1058 chainCosmic->SetAlias("dr1","(signalTot0.fElements[1]/signalTot0.fElements[3])");
1059 chainCosmic->SetAlias("dr2","(signalTot0.fElements[2]/signalTot0.fElements[3])");
1060 chainCosmic->SetAlias("dr4","(signalTot0.fElements[4]/signalTot0.fElements[3])");
1061 chainCosmic->SetAlias("dr5","(signalTot0.fElements[5]/signalTot0.fElements[3])");
1069 fstring+="dr1*dr2++";
1070 fstring+="dr1*dr4++";
1071 fstring+="dr1*dr5++";
1072 fstring+="dr2*dr4++";
1073 fstring+="dr2*dr5++";
1074 fstring+="dr4*dr5++";
1078 TString *strqdedx = toolkit.FitPlane(chain,"norm",fstring->Data(), cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
1080 chainCosmic->SetAlias("corQT",strqdedx->Data());
1086 chainCosmic->SetProof(kTRUE);
1087 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);
1090 chainCosmic->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);
1096 chainCosmic->Draw("Tr0.fP[1]-Tr1.fP[1]:sqrt(abs(Tr0.fP[4]))>>hisdzA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0&&abs(mag)>0.1"+cutA);
1098 TGraph *grdzA = (TGraph*)gProof->GetOutputList()->At(1)->Clone();