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");
38 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
39 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
40 AliXRDPROOFtoolkit tool;
41 TChain * chainCosmic = tool.MakeChainRandom("cosmic.txt","Track0",0,10000);
42 TChain * chainBudget = tool.MakeChainRandom("cosmic.txt","material",0,1000);
43 TCut cutptV="abs(1/pt0V-1/pt1V)<0.1";
44 TCut cutptI="abs(1/pt0In-1/pt1In)<0.5";
45 TCut cutncl="nclmin>120";
46 TCut cutDz="abs(p0.fP[1])<50";
47 TCut cutDr="abs(p0.fP[0])<50";
49 chainBudget->Draw(">>listB",cutptV+cutptI+cutncl+cutDr+cutDz,"entryList");
50 TEntryList *elistB = (TEntryList*)gDirectory->Get("listB");
51 chainBudget->SetEntryList(elistB);
53 chainBudget->SetAlias("dptrel","(pt0V-pt1V)/((pt0V+pt1V)*0.5)");
54 chainBudget->SetAlias("dptInrel","(pt0In-pt1In)/((pt0In+pt1In)*0.5)");
55 chainBudget->SetAlias("ptcorr","(pt0In-pt0V)/(pt0V)+(pt1V-pt1In)/(pt1In)");
60 #include "Riostream.h"
70 #include "THnSparse.h"
72 #include "AliTPCclusterMI.h"
73 #include "AliTPCseed.h"
74 #include "AliESDVertex.h"
75 #include "AliESDEvent.h"
76 #include "AliESDfriend.h"
77 #include "AliESDInputHandler.h"
78 #include "AliAnalysisManager.h"
80 #include "AliTracker.h"
82 #include "AliTPCCalROC.h"
86 #include "AliTPCcalibCosmic.h"
87 #include "TTreeStream.h"
88 #include "AliTPCTracklet.h"
89 #include "AliESDcosmic.h"
92 ClassImp(AliTPCcalibCosmic)
95 AliTPCcalibCosmic::AliTPCcalibCosmic()
104 fCutMaxD(5), // maximal distance in rfi ditection
105 fCutMaxDz(40), // maximal distance in z ditection
106 fCutTheta(0.03), // maximal distan theta
107 fCutMinDir(-0.99) // direction vector products
109 AliInfo("Default Constructor");
110 for (Int_t ihis=0; ihis<6;ihis++){
114 for (Int_t ihis=0; ihis<4;ihis++){
115 fHistodEdxMax[ihis] =0;
116 fHistodEdxTot[ihis] =0;
121 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
130 fCutMaxD(5), // maximal distance in rfi ditection
131 fCutMaxDz(40), // maximal distance in z ditection
132 fCutTheta(0.03), // maximal distan theta
133 fCutMinDir(-0.99) // direction vector products
138 fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
139 fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
140 fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
141 fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
142 fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
144 fDeDxMIP = new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
146 AliInfo("Non Default Constructor");
150 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
154 for (Int_t ihis=0; ihis<6;ihis++){
155 delete fHistoDelta[ihis];
156 delete fHistoPull[ihis];
158 for (Int_t ihis=0; ihis<4;ihis++){
159 delete fHistodEdxTot[ihis];
160 delete fHistodEdxMax[ihis];
163 delete fHistNTracks; // histogram showing number of ESD tracks per event
164 delete fClusters; // histogram showing the number of clusters per track
165 delete fModules; // 2d histogram of tracks which are propagated to the ACORDE scintillator array
166 delete fHistPt; // Pt histogram of reconstructed tracks
167 delete fDeDx; // dEdx spectrum showing the different particle types
168 delete fDeDxMIP; // TPC signal close to the MIP region of muons 0.4 < p < 0.45 GeV
172 void AliTPCcalibCosmic::Init(){
175 // Make performance histograms
178 // tracking performance bins
179 // 0 - delta of interest
180 // 1 - min (track0, track1) number of clusters
181 // 2 - R - vertex radius
183 // 4 - P2 - snp(phi) at inner wall of TPC
184 // 5 - P3 - tan(theta) at inner wall of TPC
185 // 6 - P4 - 1/pt mean
189 Double_t xminTrack[9], xmaxTrack[9];
194 axisName[0] ="#Delta";
197 xminTrack[1] =80; xmaxTrack[1]=160;
198 axisName[1] ="N_{cl}";
201 xminTrack[2] =0; xmaxTrack[2]=90; //
202 axisName[2] ="dca_{r} (cm)";
205 xminTrack[3] =-250; xmaxTrack[3]=250; //
206 axisName[3] ="z (cm)";
209 xminTrack[4] =-0.8; xmaxTrack[4]=0.8; //
210 axisName[4] ="sin(#phi)";
213 xminTrack[5] =-1; xmaxTrack[5]=1; //
214 axisName[5] ="tan(#theta)";
217 xminTrack[6] =-2; xmaxTrack[6]=2; //
218 axisName[6] ="1/pt (1/GeV)";
221 xminTrack[7] =0.2; xmaxTrack[7]=50; //
222 axisName[7] ="pt (GeV)";
225 xminTrack[8] =0; xmaxTrack[8]=TMath::Pi(); //
226 axisName[8] ="alpha";
229 xminTrack[0] =-1; xmaxTrack[0]=1; //
230 fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
231 xminTrack[0] =-5; xmaxTrack[0]=5; //
232 fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
235 xminTrack[0] =-1; xmaxTrack[0]=1; //
236 fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
237 xminTrack[0] =-5; xmaxTrack[0]=5; //
238 fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
241 xminTrack[0] =-10; xmaxTrack[0]=10; //
242 fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", 9, binsTrack,xminTrack, xmaxTrack);
243 xminTrack[0] =-5; xmaxTrack[0]=5; //
244 fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
247 xminTrack[0] =-10; xmaxTrack[0]=10; //
248 fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", 9, binsTrack,xminTrack, xmaxTrack);
249 xminTrack[0] =-5; xmaxTrack[0]=5; //
250 fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
253 xminTrack[0] =-0.2; xmaxTrack[0]=0.2; //
254 fHistoDelta[4] = new THnSparseS("#Delta_{1/pt} (1/GeV)","#Delta_{1/pt} (1/GeV)", 9, binsTrack,xminTrack, xmaxTrack);
255 xminTrack[0] =-5; xmaxTrack[0]=5; //
256 fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
260 xminTrack[0] =-0.5; xmaxTrack[0]=0.5; //
261 fHistoDelta[5] = new THnSparseS("#Delta_{pt}/p_{t}","#Delta_{pt}/p_{t}", 9, binsTrack,xminTrack, xmaxTrack);
262 xminTrack[0] =-5; xmaxTrack[0]=5; //
263 fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
266 for (Int_t idedx=0;idedx<4;idedx++){
267 xminTrack[0] =0.5; xmaxTrack[0]=1.5; //
269 xminTrack[1] =10; xmaxTrack[1]=160;
271 fHistodEdxMax[idedx] = new THnSparseS(Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
272 Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
273 9, binsTrack,xminTrack, xmaxTrack);
274 fHistodEdxTot[idedx] = new THnSparseS(Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
275 Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
276 9, binsTrack,xminTrack, xmaxTrack);
281 for (Int_t ivar=0;ivar<6;ivar++){
282 for (Int_t ivar2=0;ivar2<9;ivar2++){
283 fHistoDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
284 fHistoDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
285 fHistoPull[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
286 fHistoPull[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
287 BinLogX(fHistoDelta[ivar],7);
288 BinLogX(fHistoPull[ivar],7);
290 fHistodEdxMax[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
291 fHistodEdxMax[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
292 fHistodEdxTot[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
293 fHistodEdxTot[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
294 BinLogX(fHistodEdxMax[ivar],7);
295 BinLogX(fHistodEdxTot[ivar],7);
302 void AliTPCcalibCosmic::Add(const AliTPCcalibCosmic* cosmic){
306 for (Int_t ivar=0; ivar<6;ivar++){
307 if (fHistoDelta[ivar] && cosmic->fHistoDelta[ivar]){
308 fHistoDelta[ivar]->Add(cosmic->fHistoDelta[ivar]);
310 if (fHistoPull[ivar] && cosmic->fHistoPull[ivar]){
311 fHistoPull[ivar]->Add(cosmic->fHistoPull[ivar]);
314 for (Int_t ivar=0; ivar<4;ivar++){
315 if (fHistodEdxMax[ivar] && cosmic->fHistodEdxMax[ivar]){
316 fHistodEdxMax[ivar]->Add(cosmic->fHistodEdxMax[ivar]);
318 if (fHistodEdxTot[ivar] && cosmic->fHistodEdxTot[ivar]){
319 fHistodEdxTot[ivar]->Add(cosmic->fHistodEdxTot[ivar]);
327 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
332 Printf("ERROR: ESD not available");
335 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
337 Printf("ERROR: esdFriend not available");
342 FindPairs(event); // nearly everything takes place in find pairs...
344 if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
345 Int_t ntracks=event->GetNumberOfTracks();
346 fHistNTracks->Fill(ntracks);
347 if (ntracks==0) return;
348 // AliESDcosmic cosmicESD;
349 // TTreeSRedirector * cstream = GetDebugStreamer();
350 // cosmicESD.SetDebugStreamer(cstream);
351 // cosmicESD.ProcessEvent(event);
352 // if (cstream) cosmicESD.DumpToTree();
358 void AliTPCcalibCosmic::FillHistoPerformance(AliExternalTrackParam *par0, AliExternalTrackParam *par1, AliExternalTrackParam *inner0, AliExternalTrackParam */*inner1*/, AliTPCseed *seed0, AliTPCseed *seed1, AliExternalTrackParam *param0Combined ){
360 // par0,par1 - parameter of tracks at DCA 0
361 // inner0,inner1 - parameter of tracks at the TPC entrance
362 // seed0, seed1 - detailed track information
363 // param0Combined - Use combined track parameters for binning
365 Int_t kMinCldEdx =20;
366 Int_t ncl0 = seed0->GetNumberOfClusters();
367 Int_t ncl1 = seed1->GetNumberOfClusters();
369 const Double_t kpullCut = 10;
375 Double_t radius0 = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]);
376 Double_t radius1 = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
377 inner0->GetXYZ(xyz0);
378 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
380 x[1] = TMath::Min(ncl0,ncl1);
381 x[2] = (radius0+radius1)*0.5;
382 x[3] = param0Combined->GetZ();
383 x[4] = inner0->GetSnp();
384 x[5] = param0Combined->GetTgl();
385 x[6] = param0Combined->GetSigned1Pt();
386 x[7] = param0Combined->Pt();
391 delta[0] = (par0->GetY()+par1->GetY());
392 delta[1] = (par0->GetZ()-par1->GetZ());
393 delta[2] = (par0->GetAlpha()-par1->GetAlpha()-TMath::Pi());
394 delta[3] = (par0->GetTgl()+par1->GetTgl());
395 delta[4] = (par0->GetParameter()[4]+par1->GetParameter()[4]);
396 delta[5] = (par0->Pt()-par1->Pt())/((par0->Pt()+par1->Pt())*0.5);
398 sigma[0] = TMath::Sqrt(par0->GetSigmaY2()+par1->GetSigmaY2());
399 sigma[1] = TMath::Sqrt(par0->GetSigmaZ2()+par1->GetSigmaZ2());
400 sigma[2] = TMath::Sqrt(par0->GetSigmaSnp2()+par1->GetSigmaSnp2());
401 sigma[3] = TMath::Sqrt(par0->GetSigmaTgl2()+par1->GetSigmaTgl2());
402 sigma[4] = TMath::Sqrt(par0->GetSigma1Pt2()+par1->GetSigma1Pt2());
403 sigma[5] = sigma[4]*((par0->Pt()+par1->Pt())*0.5);
406 for (Int_t ivar=0;ivar<6;ivar++){
407 if (sigma[ivar]==0) isOK=kFALSE;
408 x[0]= delta[ivar]/sigma[ivar];
409 if (TMath::Abs(x[0])>kpullCut) isOK = kFALSE;
413 if (isOK) for (Int_t ivar=0;ivar<6;ivar++){
414 x[0]= delta[ivar]/TMath::Sqrt(2);
415 if (ivar==2 || ivar ==3) x[0]*=1000;
416 fHistoDelta[ivar]->Fill(x);
418 x[0]= delta[ivar]/sigma[ivar];
419 fHistoPull[ivar]->Fill(x);
424 // Fill dedx performance
426 for (Int_t ipad=0; ipad<4;ipad++){
432 if (ipad==0) row1=63;
433 if (ipad==1) {row0=63; row1=63+64;}
434 if (ipad==2) {row0=128;}
435 Int_t nclUp = TMath::Nint(seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
436 Int_t nclDown = TMath::Nint(seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
437 Int_t minCl = TMath::Min(nclUp,nclDown);
438 if (minCl<kMinCldEdx) continue;
441 Float_t dEdxTotUp = seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
442 Float_t dEdxTotDown = seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
443 Float_t dEdxMaxUp = seed0->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
444 Float_t dEdxMaxDown = seed1->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
446 if (dEdxTotDown<=0) continue;
447 if (dEdxMaxDown<=0) continue;
448 x[0]=dEdxTotUp/dEdxTotDown;
449 fHistodEdxTot[ipad]->Fill(x);
450 x[0]=dEdxMaxUp/dEdxMaxDown;
451 fHistodEdxMax[ipad]->Fill(x);
458 void AliTPCcalibCosmic::MaterialBudgetDump(AliExternalTrackParam *par0, AliExternalTrackParam *par1, AliExternalTrackParam *inner0, AliExternalTrackParam *inner1, AliTPCseed *seed0, AliTPCseed *seed1, AliExternalTrackParam *param0Combined, AliExternalTrackParam *param1Combined){
460 // matrial budget AOD dump
462 // par0,par1 - parameter of tracks at DCA 0
463 // inner0,inner1 - parameter of tracks at the TPC entrance
464 // seed0, seed1 - detailed track information
465 // param0Combined - Use combined track parameters for binning
467 Double_t p0In = inner0->GetP();
468 Double_t p1In = inner1->GetP();
469 Double_t p0V = par0->GetP();
470 Double_t p1V = par1->GetP();
472 Double_t pt0In = inner0->Pt();
473 Double_t pt1In = inner1->Pt();
474 Double_t pt0V = par0->Pt();
475 Double_t pt1V = par1->Pt();
476 Int_t ncl0 = seed0->GetNumberOfClusters();
477 Int_t ncl1 = seed1->GetNumberOfClusters();
478 Int_t nclmin=TMath::Min(ncl0,ncl1);
479 Double_t sign = (param0Combined->GetSigned1Pt()>0) ? 1:-1.;
481 TTreeSRedirector * pcstream = GetDebugStreamer();
483 (*pcstream)<<"material"<<
484 "run="<<fRun<< // run number
485 "event="<<fEvent<< // event number
486 "time="<<fTime<< // time stamp of event
487 "trigger="<<fTrigger<< // trigger
488 "triggerClass="<<&fTriggerClass<< // trigger
489 "mag="<<fMagF<< // magnetic field
490 "sign="<<sign<< // sign of the track
506 "up0.="<<param0Combined<<
507 "up1.="<<param1Combined<<
513 void AliTPCcalibCosmic::Analyze() {
515 fMIPvalue = CalculateMIPvalue(fDeDxMIP);
523 void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
527 // Track0 is choosen in upper TPC part
528 // Track1 is choosen in lower TPC part
530 if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
531 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
532 Int_t ntracks=event->GetNumberOfTracks();
533 TObjArray tpcSeeds(ntracks);
534 if (ntracks==0) return;
535 Double_t vtxx[3]={0,0,0};
536 Double_t svtxx[3]={0.000001,0.000001,100.};
537 AliESDVertex vtx(vtxx,svtxx);
541 for (Int_t i=0;i<ntracks;++i) {
542 AliESDtrack *track = event->GetTrack(i);
543 fClusters->Fill(track->GetTPCNcls());
545 const AliExternalTrackParam * trackIn = track->GetInnerParam();
546 const AliExternalTrackParam * trackOut = track->GetOuterParam();
547 if (!trackIn) continue;
548 if (!trackOut) continue;
549 if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // filter laser
552 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
553 TObject *calibObject;
554 AliTPCseed *seed = 0;
555 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
556 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
558 if (seed) tpcSeeds.AddAt(seed,i);
560 Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
561 if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
562 fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159));
564 if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159));
566 if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159)>300) {
567 TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
568 if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
569 if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
576 if (ntracks<2) return;
580 for (Int_t i=0;i<ntracks;++i) {
581 AliESDtrack *track0 = event->GetTrack(i);
582 // track0 - choosen upper part
583 if (!track0) continue;
584 if (!track0->GetOuterParam()) continue;
585 if (track0->GetOuterParam()->GetAlpha()<0) continue;
587 track0->GetDirection(dir0);
588 for (Int_t j=0;j<ntracks;++j) {
590 AliESDtrack *track1 = event->GetTrack(j);
592 if (!track1) continue;
593 if (!track1->GetOuterParam()) continue;
594 if (track1->GetOuterParam()->GetAlpha()>0) continue;
597 track1->GetDirection(dir1);
599 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
600 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
601 if (! seed0) continue;
602 if (! seed1) continue;
603 Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159);
604 Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159);
606 Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63);
607 Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63);
609 Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159);
610 Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159);
612 Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
613 Float_t d0 = track0->GetLinearD(0,0);
614 Float_t d1 = track1->GetLinearD(0,0);
616 // conservative cuts - convergence to be guarantied
617 // applying before track propagation
618 if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0
619 if (dir>fCutMinDir) continue; // direction vector product
620 Float_t bz = AliTracker::GetBz();
621 Float_t dvertex0[2]; //distance to 0,0
622 Float_t dvertex1[2]; //distance to 0,0
623 track0->GetDZ(0,0,0,bz,dvertex0);
624 track1->GetDZ(0,0,0,bz,dvertex1);
625 if (TMath::Abs(dvertex0[1])>250) continue;
626 if (TMath::Abs(dvertex1[1])>250) continue;
630 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
631 AliExternalTrackParam param0(*track0);
632 AliExternalTrackParam param1(*track1);
634 // Propagate using Magnetic field and correct fo material budget
638 Double_t maxsnp=0.90;
639 AliTracker::PropagateTrackToBxByBz(¶m0,dmax+1,0.0005,3,kTRUE,maxsnp,sign0);
640 AliTracker::PropagateTrackToBxByBz(¶m1,dmax+1,0.0005,3,kTRUE,maxsnp,sign1);
642 // Propagate rest to the 0,0 DCA - z should be ignored
644 Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
645 Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
647 param0.GetDZ(0,0,0,bz,dvertex0);
648 param1.GetDZ(0,0,0,bz,dvertex1);
649 if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
651 Double_t xyz0[3];//,pxyz0[3];
652 Double_t xyz1[3];//,pxyz1[3];
655 Bool_t isPair = IsPair(¶m0,¶m1);
657 if (isPair) FillAcordeHist(track0);
659 // combined track params
661 AliExternalTrackParam *par0U=MakeCombinedTrack(¶m0,¶m1);
662 AliExternalTrackParam *par1U=MakeCombinedTrack(¶m1,¶m0);
667 TTreeSRedirector * cstream = GetDebugStreamer();
668 //printf("My stream=%p\n",(void*)cstream);
669 AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
670 AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
671 AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
672 AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
673 Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
674 Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
675 Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
676 Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
680 FillHistoPerformance(¶m0, ¶m1, ip0, ip1, seed0, seed1,par0U);
681 MaterialBudgetDump(¶m0, ¶m1, ip0, ip1, seed0, seed1,par0U,par1U);
683 (*cstream) << "Track0" <<
684 "run="<<fRun<< // run number
685 "event="<<fEvent<< // event number
686 "time="<<fTime<< // time stamp of event
687 "trigger="<<fTrigger<< // trigger
688 "triggerClass="<<&fTriggerClass<< // trigger
689 "mag="<<fMagF<< // magnetic field
690 "dir="<<dir<< // direction
691 "OK="<<isPair<< // will be accepted
692 "b0="<<b0<< // propagate status
693 "b1="<<b1<< // propagate status
694 "crossI="<<isCrossI<< // cross inner
695 "crossO="<<isCrossO<< // cross outer
697 "Orig0.=" << track0 << // original track 0
698 "Orig1.=" << track1 << // original track 1
699 "Tr0.="<<¶m0<< // track propagated to the DCA 0,0
700 "Tr1.="<<¶m1<< // track propagated to the DCA 0,0
701 "Ip0.="<<ip0<< // inner param - upper
702 "Ip1.="<<ip1<< // inner param - lower
703 "Op0.="<<op0<< // outer param - upper
704 "Op1.="<<op1<< // outer param - lower
705 "Up0.="<<par0U<< // combined track 0
706 "Up1.="<<par1U<< // combined track 1
708 "v00="<<dvertex0[0]<< // distance using kalman
709 "v01="<<dvertex0[1]<< //
710 "v10="<<dvertex1[0]<< //
711 "v11="<<dvertex1[1]<< //
712 "d0="<<d0<< // linear distance to 0,0
713 "d1="<<d1<< // linear distance to 0,0
717 "x00="<<xyz0[0]<< // global position close to vertex
721 "x10="<<xyz1[0]<< // global position close to vertex
727 "dir00="<<dir0[0]<< // direction upper
731 "dir10="<<dir1[0]<< // direction lower
736 "Seed0.=" << seed0 << // original seed 0
737 "Seed1.=" << seed1 << // original seed 1
739 "dedx0="<<dedx0<< // dedx0 - all
740 "dedx1="<<dedx1<< // dedx1 - all
742 "dedx0I="<<dedx0I<< // dedx0 - inner ROC
743 "dedx1I="<<dedx1I<< // dedx1 - inner ROC
745 "dedx0O="<<dedx0O<< // dedx0 - outer ROC
746 "dedx1O="<<dedx1O<< // dedx1 - outer ROC
759 void AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
761 // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
762 if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
764 const Double_t acordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
765 const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
768 upperTrack->GetXYZ(r);
770 upperTrack->GetDirection(d);
772 z = r[2] + (d[2]/d[1])*(acordePlane - r[1]);
773 x = r[0] + (d[0]/d[1])*(acordePlane - r[1]);
776 x = r[0] + (d[0]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
777 z = r[2] + (d[2]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
780 x = r[0] + (d[0]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
781 z = r[2] + (d[2]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
784 fModules->Fill(z, x);
790 Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
792 TIterator* iter = li->MakeIterator();
793 AliTPCcalibCosmic* cal = 0;
795 while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
796 if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
797 //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
801 fHistNTracks->Add(cal->GetHistNTracks());
802 fClusters->Add(cal-> GetHistClusters());
803 fModules->Add(cal->GetHistAcorde());
804 fHistPt->Add(cal->GetHistPt());
805 fDeDx->Add(cal->GetHistDeDx());
806 fDeDxMIP->Add(cal->GetHistMIP());
814 Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
818 // 0. Same direction - OPOSITE - cutDir +cutT
819 TCut cutDir("cutDir","dir<-0.99")
821 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
824 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
828 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
831 const Double_t *p0 = tr0->GetParameter();
832 const Double_t *p1 = tr1->GetParameter();
833 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
834 if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
835 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
837 Double_t d0[3], d1[3];
838 tr0->GetDirection(d0);
839 tr1->GetDirection(d1);
840 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
847 Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
849 TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
850 funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
851 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
852 hist->Fit(funcDoubleGaus);
853 Double_t MIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
855 delete funcDoubleGaus;
864 void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) {
866 // Not implemented yet
873 void AliTPCcalibCosmic::BinLogX(THnSparse *h, Int_t axisDim) {
875 // Method for the correct logarithmic binning of histograms
877 TAxis *axis = h->GetAxis(axisDim);
878 int bins = axis->GetNbins();
880 Double_t from = axis->GetXmin();
881 Double_t to = axis->GetXmax();
882 Double_t *new_bins = new Double_t[bins + 1];
885 Double_t factor = pow(to/from, 1./bins);
887 for (int i = 1; i <= bins; i++) {
888 new_bins[i] = factor * new_bins[i-1];
890 axis->Set(bins, new_bins);
896 void AliTPCcalibCosmic::BinLogX(TH1 *h) {
898 // Method for the correct logarithmic binning of histograms
900 TAxis *axis = h->GetXaxis();
901 int bins = axis->GetNbins();
903 Double_t from = axis->GetXmin();
904 Double_t to = axis->GetXmax();
905 Double_t *new_bins = new Double_t[bins + 1];
908 Double_t factor = pow(to/from, 1./bins);
910 for (int i = 1; i <= bins; i++) {
911 new_bins[i] = factor * new_bins[i-1];
913 axis->Set(bins, new_bins);
919 AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
923 AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
924 par1R->Rotate(track0->GetAlpha());
925 par1R->PropagateTo(track0->GetX(),AliTracker::GetBz());
928 Double_t * param = (Double_t*)par1R->GetParameter();
929 Double_t * covar = (Double_t*)par1R->GetCovariance();
937 covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
938 //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
943 AliExternalTrackParam *AliTPCcalibCosmic::MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
945 // Make combined track
948 AliExternalTrackParam * par1T = MakeTrack(track0,track1);
949 AliExternalTrackParam * par0U = new AliExternalTrackParam(*track0);
951 UpdateTrack(*par0U,*par1T);
957 void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
959 // Update track 1 with track 2
963 TMatrixD vecXk(5,1); // X vector
964 TMatrixD covXk(5,5); // X covariance
965 TMatrixD matHk(5,5); // vector to mesurement
966 TMatrixD measR(5,5); // measurement error
967 TMatrixD vecZk(5,1); // measurement
969 TMatrixD vecYk(5,1); // Innovation or measurement residual
970 TMatrixD matHkT(5,5);
971 TMatrixD matSk(5,5); // Innovation (or residual) covariance
972 TMatrixD matKk(5,5); // Optimal Kalman gain
973 TMatrixD mat1(5,5); // update covariance matrix
974 TMatrixD covXk2(5,5); //
975 TMatrixD covOut(5,5);
977 Double_t *param1=(Double_t*) track1.GetParameter();
978 Double_t *covar1=(Double_t*) track1.GetCovariance();
979 Double_t *param2=(Double_t*) track2.GetParameter();
980 Double_t *covar2=(Double_t*) track2.GetCovariance();
982 // copy data to the matrix
983 for (Int_t ipar=0; ipar<5; ipar++){
984 for (Int_t jpar=0; jpar<5; jpar++){
985 covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
986 measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
990 vecXk(ipar,0) = param1[ipar];
991 vecZk(ipar,0) = param2[ipar];
1000 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
1001 matHkT=matHk.T(); matHk.T();
1002 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
1004 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
1005 vecXk += matKk*vecYk; // updated vector
1006 covXk2 = (mat1-(matKk*matHk));
1007 covOut = covXk2*covXk;
1011 // copy from matrix to parameters
1024 for (Int_t ipar=0; ipar<5; ipar++){
1025 param1[ipar]= vecXk(ipar,0) ;
1026 for (Int_t jpar=0; jpar<5; jpar++){
1027 covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);