]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/Calib/AliTPCcalibV0.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changes by Christian Holm
[u/mrichter/AliRoot.git] / TPC / Calib / AliTPCcalibV0.cxx
CommitLineData
10757ee9 1
2#include <TROOT.h>
3#include <TChain.h>
4#include <TFile.h>
5#include <TH3F.h>
5b00528f 6#include <TH2F.h>
10757ee9 7//
8#include "TParticle.h"
9#include "TDatabasePDG.h"
10#include "AliRunLoader.h"
11#include "AliStack.h"
12
13
14
15#include <TPDGCode.h>
16#include <TStyle.h>
17#include "TLinearFitter.h"
18#include "TMatrixD.h"
19#include "TTreeStream.h"
20#include "TF1.h"
21
22
23
24#include "AliMagF.h"
25#include "AliTracker.h"
57e4988a 26#include "AliESDEvent.h"
10757ee9 27#include "AliESDtrack.h"
28#include "AliESDfriend.h"
29#include "AliESDfriendTrack.h"
b40afa5e 30#include "AliMathBase.h"
10757ee9 31#include "AliTPCseed.h"
32#include "AliTPCclusterMI.h"
33
34#include "AliKFParticle.h"
35#include "AliKFVertex.h"
36
37#include "AliTrackPointArray.h"
10757ee9 38#include "AliTPCcalibV0.h"
39#include "AliV0.h"
46e89793 40#include "TRandom.h"
41#include "TTreeStream.h"
42#include "AliTPCcalibDB.h"
43#include "AliTPCCorrection.h"
44#include "AliGRPObject.h"
45#include "AliTPCTransform.h"
30e59eac 46#include "AliAnalysisManager.h"
10757ee9 47
48
49
50ClassImp(AliTPCcalibV0)
51
52
53AliTPCcalibV0::AliTPCcalibV0() :
57e4988a 54 AliTPCcalibBase(),
46e89793 55 fV0Tree(0),
56 fHPTTree(0),
10757ee9 57 fStack(0),
10757ee9 58 fESD(0),
59 fPdg(0),
60 fParticles(0),
61 fV0s(0),
46e89793 62 fGammas(0)
63{
64
65}
66AliTPCcalibV0::AliTPCcalibV0(const Text_t *name, const Text_t *title):
67 AliTPCcalibBase(),
68 fV0Tree(0),
69 fHPTTree(0),
70 fStack(0),
71 fESD(0),
72 fPdg(0),
73 fParticles(0),
74 fV0s(0),
75 fGammas(0)
10757ee9 76{
57e4988a 77 fPdg = new TDatabasePDG;
10757ee9 78 // create output histograms
46e89793 79 SetName(name);
80 SetTitle(title);
10757ee9 81}
82
83AliTPCcalibV0::~AliTPCcalibV0(){
84 //
85 //
86 //
46e89793 87 delete fV0Tree;
88 delete fHPTTree;
10757ee9 89}
90
91
92
10757ee9 93
94
30e59eac 95void AliTPCcalibV0::ProcessESD(AliESDEvent *esd){
10757ee9 96 //
97 //
98 //
99 fESD = esd;
100 AliKFParticle::SetField(esd->GetMagneticField());
46e89793 101 if (TMath::Abs(AliTracker::GetBz())<1) return;
102 DumpToTree(esd);
103 DumpToTreeHPT(esd);
10757ee9 104}
105
46e89793 106void AliTPCcalibV0::DumpToTreeHPT(AliESDEvent *esd){
107 //
108 // Dump V0s fith full firend information to the
109 //
110 if (TMath::Abs(AliTracker::GetBz())<1) return;
111 const Int_t kMinCluster=110;
30e59eac 112 const Float_t kMinPt =4.;
46e89793 113 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
30e59eac 114// if (!esdFriend) {
115// Printf("ERROR: esdFriend not available");
116// return;
117// }
46e89793 118 //
119 Int_t ntracks=esd->GetNumberOfTracks();
120 for (Int_t i=0;i<ntracks;++i) {
121 Bool_t isOK=kFALSE;
122 AliESDtrack *track = esd->GetTrack(i);
123 if (track->GetTPCncls()<kMinCluster) continue;
46e89793 124 if (TMath::Abs(AliTracker::GetBz())>1){ // cut on momenta if measured
125 if (track->Pt()>kMinPt) isOK=kTRUE;
126 }
127 if (TMath::Abs(AliTracker::GetBz())<1){ // require primary track for the B field OFF data
128 Bool_t isAccepted=kTRUE;
129 if (!track->IsOn(AliESDtrack::kITSrefit)) isAccepted=kFALSE;
130 if (!track->IsOn(AliESDtrack::kTPCrefit)) isAccepted=kFALSE;
131 if (!track->IsOn(AliESDtrack::kTOFout)) isAccepted=kFALSE;
132 Float_t dvertex[2],cvertex[3];
133 track->GetImpactParametersTPC(dvertex,cvertex);
134 if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>20) isAccepted=kFALSE;
135 if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>20) isAccepted=kFALSE;
136 track->GetImpactParameters(dvertex,cvertex);
137 if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>10) isAccepted=kFALSE;
138 if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>10) isAccepted=kFALSE;
139 if (!isAccepted) isOK=kFALSE;
30e59eac 140 }
141 if ( track->GetTPCsignal()>100 && track->GetInnerParam()->Pt()>1 ){
142 if (track->IsOn(AliESDtrack::kITSin)||track->IsOn(AliESDtrack::kTRDout)||track->IsOn(AliESDtrack::kTOFin))
143 isOK=kTRUE;
144 if (isOK){
145 TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
146 Int_t eventNumber = esd->GetEventNumberInFile();
147 Bool_t hasFriend=(esdFriend) ? (esdFriend->GetTrack(i)!=0):0;
148 Bool_t hasITS=(track->GetNcls(0)>2);
149 printf("DUMPIONTrack:%s|%f|%d|%d|%d\n",filename.Data(),track->GetInnerParam()->Pt()*track->GetTPCsignal()/50., eventNumber,hasFriend,hasITS);
150 }
46e89793 151 }
30e59eac 152 if (!isOK) continue;
153 TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
154 Int_t eventNumber = esd->GetEventNumberInFile();
155 Bool_t hasFriend=(esdFriend) ? (esdFriend->GetTrack(i)!=0):0;
156 Bool_t hasITS=(track->GetNcls(0)>2);
157 printf("DUMPHPTTrack:%s|%f|%d|%d|%d\n",filename.Data(),track->Pt(), eventNumber,hasFriend,hasITS);
158 //
159 if (!esdFriend) continue;
160 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
161 if (!friendTrack) continue;
162
46e89793 163 if (!isOK) continue;
164 //
165
166 TObject *calibObject;
167 AliTPCseed *seed = 0;
168 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
169 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
170 }
171 if (!seed) continue;
172 if (!fHPTTree) {
173 fHPTTree = new TTree("HPT","HPT");
174 fHPTTree->SetDirectory(0);
175 }
176 if (fHPTTree->GetEntries()==0){
177 //
178 fHPTTree->SetDirectory(0);
179 fHPTTree->Branch("t.",&track);
180 fHPTTree->Branch("ft.",&friendTrack);
181 fHPTTree->Branch("s.",&seed);
182 }else{
183 fHPTTree->SetBranchAddress("t.",&track);
184 fHPTTree->SetBranchAddress("ft.",&friendTrack);
185 fHPTTree->SetBranchAddress("s.",&seed);
186 }
187 fHPTTree->Fill();
188 //
189 }
190}
191
192
193
194void AliTPCcalibV0::DumpToTree(AliESDEvent *esd){
195 //
196 // Dump V0s fith full firend information to the
197 //
198 Int_t nV0s = fESD->GetNumberOfV0s();
199 const Int_t kMinCluster=110;
200 const Double_t kDownscale=0.01;
46e89793 201 const Float_t kMinPt =1.0;
202 const Float_t kMinMinPt =0.7;
203 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
46e89793 204 //
205
206 for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
207 Bool_t isOK=kFALSE;
208 AliESDv0 * v0 = (AliESDv0*) esd->GetV0(ivertex);
209 AliESDtrack * track0 = fESD->GetTrack(v0->GetIndex(0)); // negative track
210 AliESDtrack * track1 = fESD->GetTrack(v0->GetIndex(1)); // positive track
211 if (track0->GetTPCNcls()<kMinCluster) continue;
212 if (track0->GetKinkIndex(0)>0) continue;
213 if (track1->GetTPCNcls()<kMinCluster) continue;
214 if (track1->GetKinkIndex(0)>0) continue;
215 if (v0->GetOnFlyStatus()==kFALSE) continue;
216 //
217 if (TMath::Min(track0->Pt(),track1->Pt())<kMinMinPt) continue;
218 //
219 //
220 if (TMath::Max(track0->Pt(),track1->Pt())>kMinPt) isOK=kTRUE;
221 if (gRandom->Rndm()<kDownscale) isOK=kTRUE;
222 if (!isOK) continue;
30e59eac 223 //
224 TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
225 Int_t eventNumber = esd->GetEventNumberInFile();
226 Bool_t hasITS=(track0->GetNcls(0)+ track1->GetNcls(0)>4);
227 printf("DUMPHPTV0:%s|%f|%d|%d|%d\n",filename.Data(), (TMath::Min(track0->Pt(),track1->Pt())), eventNumber,(esdFriend!=0), hasITS);
228 //
229 if (!esdFriend) continue;
230 //
231
46e89793 232 //
233 AliESDfriendTrack *ftrack0 = esdFriend->GetTrack(v0->GetIndex(0));
234 if (!ftrack0) continue;
235 AliESDfriendTrack *ftrack1 = esdFriend->GetTrack(v0->GetIndex(1));
236 if (!ftrack1) continue;
237 //
238 TObject *calibObject;
239 AliTPCseed *seed0 = 0;
240 AliTPCseed *seed1 = 0;
241 for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
242 if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
243 }
244 for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
245 if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
246 }
247 if (!seed0) continue;
248 if (!seed1) continue;
249 AliExternalTrackParam * paramIn0 = (AliExternalTrackParam *)track0->GetInnerParam();
250 AliExternalTrackParam * paramIn1 = (AliExternalTrackParam *)track1->GetInnerParam();
251 if (!paramIn0) continue;
252 if (!paramIn1) continue;
253 //
254 //
255 if (!fV0Tree) {
256 fV0Tree = new TTree("V0s","V0s");
257 fV0Tree->SetDirectory(0);
258 }
259 if (fV0Tree->GetEntries()==0){
260 //
261 fV0Tree->SetDirectory(0);
262 fV0Tree->Branch("v0.",&v0);
263 fV0Tree->Branch("t0.",&track0);
264 fV0Tree->Branch("t1.",&track1);
265 fV0Tree->Branch("ft0.",&ftrack0);
266 fV0Tree->Branch("ft1.",&ftrack1);
267 fV0Tree->Branch("s0.",&seed0);
268 fV0Tree->Branch("s1.",&seed1);
269 }else{
270 fV0Tree->SetBranchAddress("v0.",&v0);
271 fV0Tree->SetBranchAddress("t0.",&track0);
272 fV0Tree->SetBranchAddress("t1.",&track1);
273 fV0Tree->SetBranchAddress("ft0.",&ftrack0);
274 fV0Tree->SetBranchAddress("ft1.",&ftrack1);
275 fV0Tree->SetBranchAddress("s0.",&seed0);
276 fV0Tree->SetBranchAddress("s1.",&seed1);
277 }
278 fV0Tree->Fill();
279 }
280}
281
282
283Long64_t AliTPCcalibV0::Merge(TCollection *const li) {
284
285 TIterator* iter = li->MakeIterator();
286 AliTPCcalibV0* cal = 0;
287
288 while ((cal = (AliTPCcalibV0*)iter->Next())) {
289 if (cal->fV0Tree){
290 if (!fV0Tree) {
291 fV0Tree = new TTree("V0s","V0s");
292 fV0Tree->SetDirectory(0);
293 }
294 if (cal->fV0Tree->GetEntries()>0) AliTPCcalibV0::AddTree(cal->fV0Tree);
295 if (cal->fHPTTree->GetEntries()>0) AliTPCcalibV0::AddTreeHPT(cal->fHPTTree);
296 }
297 }
298 return 0;
299}
300
301
302void AliTPCcalibV0::AddTree(TTree * treeInput){
303 //
304 // Add the content of tree:
305 // Notice automatic copy of tree in ROOT does not work for such complicated tree
306 //
30e59eac 307 return ;
46e89793 308 AliESDv0 * v0 = new AliESDv0;
309 Double_t kMinPt=0.8;
310 AliESDtrack * track0 = 0; // negative track
311 AliESDtrack * track1 = 0; // positive track
312 AliESDfriendTrack *ftrack0 = 0;
313 AliESDfriendTrack *ftrack1 = 0;
314 AliTPCseed *seed0 = 0;
315 AliTPCseed *seed1 = 0;
316 treeInput->SetBranchStatus("ft0.",kFALSE);
317 treeInput->SetBranchStatus("ft1.",kFALSE);
318 TDatabasePDG pdg;
319 Double_t massK0= pdg.GetParticle("K0")->Mass();
320 Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
321
322 Int_t entries= treeInput->GetEntries();
323 for (Int_t i=0; i<entries; i++){
324 treeInput->SetBranchAddress("v0.",&v0);
325 treeInput->SetBranchAddress("t0.",&track0);
326 treeInput->SetBranchAddress("t1.",&track1);
327 treeInput->SetBranchAddress("ft0.",&ftrack0);
328 treeInput->SetBranchAddress("ft1.",&ftrack1);
329 treeInput->SetBranchAddress("s0.",&seed0);
330 treeInput->SetBranchAddress("s1.",&seed1);
331 if (fV0Tree->GetEntries()==0){
332 fV0Tree->SetDirectory(0);
333 fV0Tree->Branch("v0.",&v0);
334 fV0Tree->Branch("t0.",&track0);
335 fV0Tree->Branch("t1.",&track1);
336 fV0Tree->Branch("ft0.",&ftrack0);
337 fV0Tree->Branch("ft1.",&ftrack1);
338 fV0Tree->Branch("s0.",&seed0);
339 fV0Tree->Branch("s1.",&seed1);
340 }else{
341 fV0Tree->SetBranchAddress("v0.",&v0);
342 fV0Tree->SetBranchAddress("t0.",&track0);
343 fV0Tree->SetBranchAddress("t1.",&track1);
344 fV0Tree->SetBranchAddress("ft0.",&ftrack0);
345 fV0Tree->SetBranchAddress("ft1.",&ftrack1);
346 fV0Tree->SetBranchAddress("s0.",&seed0);
347 fV0Tree->SetBranchAddress("s1.",&seed1);
348 }
349 //
350 treeInput->GetEntry(i);
f11122f7 351 //ftrack0->GetCalibContainer()->SetOwner(kTRUE);
352 //ftrack1->GetCalibContainer()->SetOwner(kTRUE);
46e89793 353 Bool_t isOK=kTRUE;
354 if (v0->GetOnFlyStatus()==kFALSE) isOK=kFALSE;
355 if (track0->GetTPCncls()<100) isOK=kFALSE;
356 if (track1->GetTPCncls()<100) isOK=kFALSE;
357 if (TMath::Min(seed0->Pt(),seed1->Pt())<kMinPt) isOK=kFALSE;
358 if (TMath::Min(track0->Pt(),track1->Pt())<kMinPt) isOK=kFALSE;
359 Bool_t isV0=kFALSE;
360 if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.05) isV0=kTRUE;
361 if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.05) isV0=kTRUE;
362 if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.05) isV0=kTRUE;
363 if (TMath::Abs(v0->GetEffMass(0,0))<0.02) isV0=kFALSE; //reject electrons
364 if (!isV0) isOK=kFALSE;
365 if (isOK) fV0Tree->Fill();
366 delete v0;
367 delete track0;
368 delete track1;
369 delete ftrack0;
370 delete ftrack1;
371 delete seed0;
372 delete seed1;
373 v0=0;
374 track0=0;
375 track1=0;
376 ftrack0=0;
377 ftrack1=0;
378 seed0=0;
379 seed1=0;
380 }
381}
382
383void AliTPCcalibV0::AddTreeHPT(TTree * treeInput){
384 //
385 // Add the content of tree:
386 // Notice automatic copy of tree in ROOT does not work for such complicated tree
387 //
30e59eac 388 return ;
46e89793 389 AliESDtrack *track = 0;
390 AliESDfriendTrack *friendTrack = 0;
391 AliTPCseed *seed = 0;
392 if (!treeInput) return;
393 if (treeInput->GetEntries()==0) return;
394 //
395 Int_t entries= treeInput->GetEntries();
396 //
397 for (Int_t i=0; i<entries; i++){
398 track=0;
399 friendTrack=0;
400 seed=0;
401 //
402 treeInput->SetBranchAddress("t.",&track);
403 treeInput->SetBranchAddress("ft.",&friendTrack);
404 treeInput->SetBranchAddress("s.",&seed);
405 treeInput->GetEntry(i);
406 //
407 if (fHPTTree->GetEntries()==0){
408 fHPTTree->SetDirectory(0);
409 fHPTTree->Branch("t.",&track);
410 fHPTTree->Branch("ft.",&friendTrack);
411 fHPTTree->Branch("s.",&seed);
412 }else{
413 fHPTTree->SetBranchAddress("t.",&track);
414 fHPTTree->SetBranchAddress("ft.",&friendTrack);
415 fHPTTree->SetBranchAddress("s.",&seed);
416 }
417 Bool_t isOK=kTRUE;
418 if (!track->IsOn(AliESDtrack::kITSrefit)) isOK=kFALSE;
419 if (!track->IsOn(AliESDtrack::kTOFout)) isOK=kFALSE;
420 if (isOK) fHPTTree->Fill();
421 //
422 delete track;
423 delete friendTrack;
424 delete seed;
425 }
426}
427
428
a8ef8a9c 429void AliTPCcalibV0::MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t /*run*/){
46e89793 430 //
431 // Make a fit tree
432 //
433 // 0. Loop over selected tracks
434 // 1. Loop over all transformation - refit the track with and without the
435 // transformtation
436 // 2. Dump the matching paramaeters to the debugStremer
437 //
438
439 //Connect input
440 const Int_t kMinNcl=120;
441 TFile f("TPCV0Objects.root");
442 AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
443 TTree * treeInput = v0TPC->GetHPTTree();
444 TTreeSRedirector *pcstream = new TTreeSRedirector("fitHPT.root");
445 AliESDtrack *track = 0;
446 AliESDfriendTrack *friendTrack = 0;
447 AliTPCseed *seed = 0;
448 if (!treeInput) return;
449 if (treeInput->GetEntries()==0) return;
450 //
451 treeInput->SetBranchAddress("t.",&track);
452 treeInput->SetBranchAddress("ft.",&friendTrack);
453 treeInput->SetBranchAddress("s.",&seed);
454 //
455 Int_t ncorr=0;
456 if (corrArray) ncorr = corrArray->GetEntries();
457 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
f2ebe917 458 // AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
459// AliGRPObject* grp = AliTPCcalibDB::Instance()->GetGRP(run);
460// Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
46e89793 461 //
462 //
463 //
464 Int_t ntracks= treeInput->GetEntries();
465 for (Int_t itrack=0; itrack<ntracks; itrack++){
466 treeInput->GetEntry(itrack);
467 if (!track) continue;
468 if (seed->Pt()<ptCut) continue;
469 if (track->Pt()<ptCut) continue;
470 if (track->GetTPCncls()<kMinNcl) continue;
471 //
472 // Reapply transformation
473 //
474 for (Int_t irow=0; irow<159; irow++){
475 AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
476 if (cluster &&cluster->GetX()>10){
477 Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
478 Int_t index0[1]={cluster->GetDetector()};
479 transform->Transform(x0,index0,0,1);
480 cluster->SetX(x0[0]);
481 cluster->SetY(x0[1]);
482 cluster->SetZ(x0[2]);
483 //
484 }
485 }
486 //
46e89793 487 AliExternalTrackParam* paramInner=0;
488 AliExternalTrackParam* paramOuter=0;
489 AliExternalTrackParam* paramIO=0;
490 Bool_t isOK=kTRUE;
491 for (Int_t icorr=-1; icorr<ncorr; icorr++){
492 //
493 AliTPCCorrection *corr = 0;
494 if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
495 AliExternalTrackParam * trackInner = RefitTrack(seed, corr,85,134,0.1);
496 AliExternalTrackParam * trackIO = RefitTrack(seed, corr,245,85,0.1);
497 AliExternalTrackParam * trackOuter = RefitTrack(seed, corr,245,134,0.1 );
498 if (trackInner&&trackOuter&&trackIO){
499 trackOuter->Rotate(trackInner->GetAlpha());
500 trackOuter->PropagateTo(trackInner->GetX(),AliTracker::GetBz());
501 if (icorr<0) {
502 paramInner=trackInner;
503 paramOuter=trackOuter;
504 paramIO=trackIO;
505 paramIO->Rotate(seed->GetAlpha());
506 paramIO->PropagateTo(seed->GetX(),AliTracker::GetBz());
507 }
508 }else{
509 isOK=kFALSE;
510 }
511
512 }
513 if (paramOuter&& paramInner) {
514 // Bool_t isOK=kTRUE;
515 if (paramInner->GetSigmaY2()>0.01) isOK&=kFALSE;
516 if (paramOuter->GetSigmaY2()>0.01) isOK&=kFALSE;
517 if (paramInner->GetSigmaZ2()>0.01) isOK&=kFALSE;
518 if (paramOuter->GetSigmaZ2()>0.01) isOK&=kFALSE;
519 (*pcstream)<<"fit"<<
520 "s.="<<seed<<
521 "io.="<<paramIO<<
522 "pIn.="<<paramInner<<
523 "pOut.="<<paramOuter;
524 }
525 //
526 }
527 delete pcstream;
528 /*
529 .x ~/rootlogon.C
530 Int_t run=117112;
531 .x ../ConfigCalibTrain.C(run)
532 .L ../AddTaskTPCCalib.C
533 ConfigOCDB(run)
534 TFile fexb("../../RegisterCorrectionExB.root");
535 AliTPCComposedCorrection *cc= (AliTPCComposedCorrection*) fexb.Get("ComposedExB");
536 cc->Init();
537 cc->Print("DA"); // Print used correction classes
538 TObjArray *array = cc->GetCorrections()
539 AliTPCcalibV0::MakeFitTreeTrack(array,2,run);
540
541 */
542}
10757ee9 543
46e89793 544void AliTPCcalibV0::MakeFitTreeV0(const TObjArray * corrArray, Double_t ptCut, Int_t run){
545 //
546 // Make a fit tree
547 //
548 // 0. Loop over selected tracks
549 // 1. Loop over all transformation - refit the track with and without the
550 // transformtation
551 // 2. Dump the matching paramaeters to the debugStremer
552 //
553
554 //Connect input
46e89793 555 TFile f("TPCV0Objects.root");
556 AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
557 TTree * treeInput = v0TPC->GetV0Tree();
558 TTreeSRedirector *pcstream = new TTreeSRedirector("fitV0.root");
559 AliESDv0 *v0 = 0;
560 AliESDtrack *track0 = 0;
561 AliESDfriendTrack *friendTrack0 = 0;
562 AliTPCseed *seed0 = 0;
563 AliTPCseed *s0 = 0;
564 AliESDtrack *track1 = 0;
565 AliESDfriendTrack *friendTrack1 = 0;
566 AliTPCseed *s1 = 0;
567 AliTPCseed *seed1 = 0;
568 if (!treeInput) return;
569 if (treeInput->GetEntries()==0) return;
570 //
571 treeInput->SetBranchAddress("v0.",&v0);
572 treeInput->SetBranchAddress("t0.",&track0);
573 treeInput->SetBranchAddress("ft0.",&friendTrack0);
574 treeInput->SetBranchAddress("s0.",&s0);
575 treeInput->SetBranchAddress("t1.",&track1);
576 treeInput->SetBranchAddress("ft1.",&friendTrack1);
577 treeInput->SetBranchAddress("s1.",&s1);
578 //
579 TDatabasePDG pdg;
580 Int_t ncorr=0;
581 if (corrArray) ncorr = corrArray->GetEntries();
582 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
46e89793 583 Double_t massK0= pdg.GetParticle("K0")->Mass();
584 Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
585 Double_t massPion=pdg.GetParticle("pi+")->Mass();
586 Double_t massProton=pdg.GetParticle("proton")->Mass();
f2ebe917 587 Int_t pdgPion=pdg.GetParticle("pi+")->PdgCode();
588 Int_t pdgProton=pdg.GetParticle("proton")->PdgCode();
589 Double_t rmass0=0;
590 Double_t rmass1=0;
46e89793 591 Double_t massV0=0;
592 Int_t pdg0=0;
593 Int_t pdg1=0;
594 //
595 //
596 //
597 Int_t nv0s= treeInput->GetEntries();
598 for (Int_t iv0=0; iv0<nv0s; iv0++){
599 Int_t v0Type=0;
600 Int_t isK0=0;
601 Int_t isLambda=0;
602 Int_t isAntiLambda=0;
603 treeInput->GetEntry(iv0);
604 if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.03) {isK0=1; v0Type=1;} //select K0s
605 if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.01) {isLambda=1; v0Type=2;} //select Lambda
606 if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.01) {isAntiLambda=1;v0Type=3;} //select Anti Lambda
607 if (isK0+isLambda+isAntiLambda!=1) continue;
f2ebe917 608 rmass0=massPion;
609 rmass1=massPion;
46e89793 610 pdg0=pdgPion;
611 pdg1=pdgPion;
f2ebe917 612 if (isLambda) {rmass0=massProton; pdg0=pdgProton;}
613 if (isAntiLambda) {rmass1=massProton; pdg1=pdgProton;}
46e89793 614 massV0=massK0;
615 if (isK0==0) massV0=massLambda;
616 //
617 if (!s0) continue;
618 seed0=(s0->GetSign()>0)?s0:s1;
619 seed1=(s0->GetSign()>0)?s1:s0;
620 if (seed0->GetZ()*seed1->GetZ()<0) continue; //remove membrane crossed tracks
621 if (seed0->Pt()<ptCut) continue;
622 if (seed1->Pt()<ptCut) continue;
623 //
624 // Reapply transformation
625 //
626 for (Int_t itype=0; itype<2; itype++){
627 AliTPCseed * seed = (itype==0) ? seed0: seed1;
628 for (Int_t irow=0; irow<159; irow++){
629 AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
630 if (cluster &&cluster->GetX()>10){
631 Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
632 Int_t index0[1]={cluster->GetDetector()};
633 transform->Transform(x0,index0,0,1);
634 cluster->SetX(x0[0]);
635 cluster->SetY(x0[1]);
636 cluster->SetZ(x0[2]);
637 //
638 }
639 }
640 }
641 Bool_t isOK=kTRUE;
642 Double_t radius = v0->GetRr();
643 Double_t xyz[3];
644 v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
645 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
646 TObjArray arrayV0in(ncorr+1);
647 TObjArray arrayV0io(ncorr+1);
648 TObjArray arrayT0(ncorr+1);
649 TObjArray arrayT1(ncorr+1);
650 arrayV0in.SetOwner(kTRUE);
651 arrayV0io.SetOwner(kTRUE);
652 //
653 for (Int_t icorr=-1; icorr<ncorr; icorr++){
654 AliTPCCorrection *corr =0;
655 if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
656 //
f2ebe917 657 AliExternalTrackParam * trackInner0 = RefitTrack(seed0, corr,160,85,rmass0);
658 AliExternalTrackParam * trackIO0 = RefitTrack(seed0, corr,245,85,rmass0);
659 AliExternalTrackParam * trackInner1 = RefitTrack(seed1, corr,160,85,rmass1);
660 AliExternalTrackParam * trackIO1 = RefitTrack(seed1, corr,245,85,rmass1);
46e89793 661 if (!trackInner0) isOK=kFALSE;
662 if (!trackInner1) isOK=kFALSE;
663 if (!trackIO0) isOK=kFALSE;
664 if (!trackIO1) isOK=kFALSE;
665 if (isOK){
666 if (!trackInner0->Rotate(alpha)) isOK=kFALSE;
667 if (!trackInner1->Rotate(alpha)) isOK=kFALSE;
668 if (!trackIO0->Rotate(alpha)) isOK=kFALSE;
669 if (!trackIO1->Rotate(alpha)) isOK=kFALSE;
670 //
f2ebe917 671 if (!AliTracker::PropagateTrackToBxByBz(trackInner0, radius, rmass0, 1, kFALSE)) isOK=kFALSE;
672 if (!AliTracker::PropagateTrackToBxByBz(trackInner1, radius, rmass1, 1, kFALSE)) isOK=kFALSE;
673 if (!AliTracker::PropagateTrackToBxByBz(trackIO0, radius, rmass0, 1, kFALSE)) isOK=kFALSE;
674 if (!AliTracker::PropagateTrackToBxByBz(trackIO1, radius, rmass1, 1, kFALSE)) isOK=kFALSE;
46e89793 675 if (!isOK) continue;
676 arrayT0.AddAt(trackIO0->Clone(),icorr+1);
677 arrayT1.AddAt(trackIO1->Clone(),icorr+1);
f2ebe917 678 Int_t charge=TMath::Nint(trackIO0->GetSign());
46e89793 679 AliKFParticle pin0( *trackInner0, pdg0*charge);
680 AliKFParticle pin1( *trackInner1, -pdg1*charge);
681 AliKFParticle pio0( *trackIO0, pdg0*charge);
682 AliKFParticle pio1( *trackIO1, -pdg1*charge);
683 AliKFParticle v0in;
684 AliKFParticle v0io;
685 v0in+=pin0;
686 v0in+=pin1;
687 v0io+=pio0;
688 v0io+=pio1;
689 arrayV0in.AddAt(v0in.Clone(),icorr+1);
690 arrayV0io.AddAt(v0io.Clone(),icorr+1);
691 }
692 }
693 if (!isOK) continue;
694 //
695 //AliKFParticle* pin0= (AliKFParticle*)arrayV0in.At(0);
696 AliKFParticle* pio0= (AliKFParticle*)arrayV0io.At(0);
697 AliExternalTrackParam *param0=(AliExternalTrackParam *)arrayT0.At(0);
698 AliExternalTrackParam *param1=(AliExternalTrackParam *)arrayT1.At(0);
699 Double_t mass0=0, mass0E=0;
700 pio0->GetMass( mass0,mass0E);
701 //
702 Double_t mean=mass0-massV0;
703 if (TMath::Abs(mean)>0.05) continue;
704 Double_t mass[10000];
705 //
706 Int_t dtype=30; // id for V0
707 Int_t ptype=5; // id for invariant mass
708 // Int_t id=TMath::Nint(100.*(param0->Pt()-param1->Pt())/(param0->Pt()+param1->Pt())); // K0s V0 asymetry
f2ebe917 709 Int_t id=Int_t(1000.*(param0->Pt()-param1->Pt())); // K0s V0 asymetry
46e89793 710 Double_t gx,gy,gz, px,py,pz;
711 Double_t pt = v0->Pt();
712 v0->GetXYZ(gx,gy,gz);
713 v0->GetPxPyPz(px,py,pz);
714 Double_t theta=pz/TMath::Sqrt(px*px+py*py);
715 Double_t phi=TMath::ATan2(py,px);
716 Double_t snp=0.5*(seed0->GetSnp()+seed1->GetSnp());
717 Double_t sector=9*phi/TMath::Pi();
718 Double_t dsec=sector-TMath::Nint(sector);
719 Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
720 //Int_t nentries=v0Type;
721 Double_t bz=AliTracker::GetBz();
722 Double_t dRrec=0;
723 (*pcstream)<<"fitDebug"<<
724 "id="<<id<<
725 "v0.="<<v0<<
726 "mean="<<mean<<
727 "rms="<<mass0E<<
728 "pio0.="<<pio0<<
729 "p0.="<<param0<<
730 "p1.="<<param1;
731 (*pcstream)<<"fit"<< // dump valus for fit
732 "run="<<run<< //run number
733 "bz="<<bz<< // magnetic filed used
734 "dtype="<<dtype<< // detector match type 30
735 "ptype="<<ptype<< // parameter type
736 "theta="<<theta<< // theta
737 "phi="<<phi<< // phi
738 "snp="<<snp<< // snp
739 "mean="<<mean<< // mean dist value
740 "rms="<<mass0E<< // rms
741 "sector="<<sector<<
742 "dsec="<<dsec<<
743 //
744 "refX="<<refX<< // reference radius
745 "gx="<<gx<< // global position
746 "gy="<<gy<< // global position
747 "gz="<<gz<< // global position
748 "dRrec="<<dRrec<< // delta Radius in reconstruction
749 "pt="<<pt<< // pt of the particle
750 "id="<<id<< //delta of the momenta
751 "entries="<<v0Type;// type of the V0
752 for (Int_t icorr=0; icorr<ncorr; icorr++){
753 AliTPCCorrection *corr =0;
754 if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
f2ebe917 755 // AliKFParticle* pin= (AliKFParticle*)arrayV0in.At(icorr+1);
46e89793 756 AliKFParticle* pio= (AliKFParticle*)arrayV0io.At(icorr+1);
757 AliExternalTrackParam *par0=(AliExternalTrackParam *)arrayT0.At(icorr+1);
758 AliExternalTrackParam *par1=(AliExternalTrackParam *)arrayT1.At(icorr+1);
759 Double_t massE=0;
760 pio->GetMass( mass[icorr],massE);
761 mass[icorr]-=mass0;
762 (*pcstream)<<"fit"<<
763 Form("%s=",corr->GetName())<<mass[icorr];
764 (*pcstream)<<"fitDebug"<<
765 Form("%s=",corr->GetName())<<mass[icorr]<<
766 Form("%sp0.=",corr->GetName())<<par0<<
767 Form("%sp1=",corr->GetName())<<par1;
768 }
769 (*pcstream)<<"fit"<< "isOK="<<isOK<<"\n";
770 (*pcstream)<<"fitDebug"<< "isOK="<<isOK<<"\n";
771 }
772 delete pcstream;
773}
10757ee9 774
46e89793 775AliExternalTrackParam * AliTPCcalibV0::RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass){
10757ee9 776 //
46e89793 777 // Refit the track:
778 // seed - tpc track with cluster
779 // corr - distrotion/correction class - apllied to the points
780 // xstart - radius to start propagate/update
781 // xstop - radius to stop propagate/update
782 //
783 const Double_t kResetCov=20.;
784 const Double_t kSigma=5.;
785 Double_t covar[15];
786 for (Int_t i=0;i<15;i++) covar[i]=0;
787 covar[0]=kSigma*kSigma;
788 covar[2]=kSigma*kSigma;
789 covar[5]=kSigma*kSigma/Float_t(150*150);
790 covar[9]=kSigma*kSigma/Float_t(150*150);
791 covar[14]=1*1;
792 //
793 AliExternalTrackParam *refit = new AliExternalTrackParam(*seed);
794 refit->PropagateTo(xstart, AliTracker::GetBz());
795 refit->AddCovariance(covar);
796 refit->ResetCovariance(kResetCov);
797 Double_t xmin = TMath::Min(xstart,xstop);
798 Double_t xmax = TMath::Max(xstart,xstop);
799 Int_t ncl=0;
10757ee9 800 //
46e89793 801 Bool_t isOK=kTRUE;
802 for (Int_t index0=0; index0<159; index0++){
803 Int_t irow= (xstart<xstop)? index0:159-index0;
804 AliTPCclusterMI *cluster=seed->GetClusterPointer(irow); //cluster in local system
805 if (!cluster) continue;
806 if (cluster->GetX()<xmin) continue;
807 if (cluster->GetX()>xmax) continue;
808 Double_t alpha = TMath::Pi()*(cluster->GetDetector()+0.5)/9.;
809 if (!refit->Rotate(alpha)) isOK=kFALSE;
810 Double_t x = cluster->GetX();
811 Double_t y = cluster->GetY();
812 Double_t z = cluster->GetZ();
813 if (corr){
814 Float_t xyz[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()}; // original position
815 corr->DistortPointLocal(xyz,cluster->GetDetector());
816 x=xyz[0];
817 y=xyz[1];
818 z=xyz[2];
819 }
820 if (!AliTracker::PropagateTrackToBxByBz(refit, x,mass,1.,kFALSE)) isOK=kFALSE;
821 if (!isOK) continue;
822 Double_t cov[3]={0.01,0.,0.01};
823 Double_t yz[2]={y,z};
824 if (!refit->Update(yz,cov)) isOK=kFALSE;
825 ncl++;
826 }
827 if (!AliTracker::PropagateTrackToBxByBz(refit, xstop, mass,1.,kTRUE)) isOK=kFALSE;
828 //
829 if (!isOK) {
830 delete refit;
831 return 0;
832 }
833 return refit;
834}
835
836
837
838
839
840//
841// Obsolete part
842//
843
844
10757ee9 845
846
57dc06f2 847AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){
10757ee9 848 //
849 // Make KF Particle
850 //
851 AliKFParticle p1( *(v0->GetParamN()), PDG1 );
852 AliKFParticle p2( *(v0->GetParamP()), PDG2 );
853 AliKFParticle *V0 = new AliKFParticle;
854 Double_t x, y, z;
855 v0->GetXYZ(x,y,z );
856 V0->SetVtxGuess(x,y,z);
857 *(V0)+=p1;
858 *(V0)+=p2;
859 return V0;
860}
861
10757ee9 862
863
864
865void AliTPCcalibV0::BinLogX(TH2F *h) {
866 //
867 //
868 //
869 TAxis *axis = h->GetXaxis();
870 int bins = axis->GetNbins();
871
872 Double_t from = axis->GetXmin();
873 Double_t to = axis->GetXmax();
874 Double_t *new_bins = new Double_t[bins + 1];
875
876 new_bins[0] = from;
877 Double_t factor = pow(to/from, 1./bins);
878
879 for (int i = 1; i <= bins; i++) {
880 new_bins[i] = factor * new_bins[i-1];
881 }
882 axis->Set(bins, new_bins);
30e59eac 883 delete [] new_bins;
10757ee9 884}
885
886
887
30e59eac 888void AliTPCcalibV0::FilterV0s(AliESDEvent* event){
889 //
890 //
891 TDatabasePDG pdg;
892 const Double_t kChi2Cut=20;
893 const Double_t kMinR=2;
894 const Double_t ptCut=0.2;
895 const Int_t kMinNcl=110;
896 //
897 Int_t nv0 = event->GetNumberOfV0s();
898 AliESDVertex *vertex= (AliESDVertex *)event->GetPrimaryVertex();
899 AliKFVertex kfvertex=*vertex;
900 //
901 for (Int_t iv0=0;iv0<nv0;iv0++){
902 AliESDv0 *v0 = event->GetV0(iv0);
903 if (!v0) continue;
904 if (v0->GetPindex()<0) continue;
905 if (v0->GetNindex()<0) continue;
906 if (TMath::Max(v0->GetPindex(), v0->GetNindex())>event->GetNumberOfTracks()) continue;
907 //
908 //
909 AliExternalTrackParam pp=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamP())):(*(v0->GetParamN()));
910 AliExternalTrackParam pn=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamN())):(*(v0->GetParamP()));
911 AliKFParticle kfp1( pp, 211 );
912 AliKFParticle kfp2( pn, -211 );
913 AliKFParticle *v0KFK0 = new AliKFParticle(kfp1,kfp2);
914 AliKFParticle *v0KFK0CV = new AliKFParticle(*v0KFK0);
915 v0KFK0CV->SetProductionVertex(kfvertex);
916 v0KFK0CV->TransportToProductionVertex();
917 AliKFParticle *v0KFK0CVM = new AliKFParticle(*v0KFK0CV);
918 v0KFK0CVM->SetMassConstraint(pdg.GetParticle("K_S0")->Mass());
919 Double_t chi2K0 = v0KFK0CV->GetChi2();
920 Double_t chi2K0M= v0KFK0CVM->GetChi2();
921 Double_t massK0=v0KFK0CV->GetMass();
922 if (chi2K0>kChi2Cut) continue;
923 if (v0->GetRr()>kMinR) continue;
924 //
925 Double_t maxPt = TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());
926 Double_t effMass22=v0->GetEffMass(2,2);
927 Double_t effMass42=v0->GetEffMass(4,2);
928 Double_t effMass24=v0->GetEffMass(2,4);
929 Bool_t isV0= kFALSE;
930 isV0|=TMath::Abs(effMass22-pdg.GetParticle("K_S0")->Mass())<0.1;
931 isV0|=TMath::Abs(effMass42-pdg.GetParticle("Lambda0")->Mass())<0.1;
932 isV0|=TMath::Abs(effMass24-pdg.GetParticle("Lambda0")->Mass())<0.1;
933
934 Double_t sign= v0->GetParamP()->GetSign()* v0->GetParamN()->GetSign();
935 if (sign<0&&v0->GetOnFlyStatus()>0.5&&maxPt>ptCut&&isV0){
936 AliESDtrack * trackP = event->GetTrack(v0->GetPindex());
937 AliESDtrack * trackN = event->GetTrack(v0->GetNindex());
938 if (!trackN) continue;
939 if (!trackP) continue;
940 Int_t nclP= (Int_t)trackP->GetTPCClusterInfo(2,1);
941 Int_t nclN= (Int_t)trackN->GetTPCClusterInfo(2,1);
942 if (TMath::Min(nclP,nclN)<kMinNcl) continue;
943 Double_t eta = TMath::Max(TMath::Abs(trackP->Eta()), TMath::Abs(trackN->Eta()));
944 Double_t ncls = TMath::Min(TMath::Abs(trackP->GetNcls(0)), TMath::Abs(trackN->GetNcls(0)));
945 if (eta<0.8&&ncls>2){
946 // printf("%d\t%f\t%f\t%d\t%d\t%f\t%f\n",i, v0->Pt(), maxPt, v0->GetNindex(),v0->GetPindex(),v0->GetRr(),effMass22);
947 (*fDebugStreamer)<<"v0tree"<<
948 "v0.="<<v0<<
949 "tp.="<<trackP<<
950 "tm.="<<trackN<<
951 //
952 "v.="<<vertex<<
953 "ncls="<<ncls<<
954 "maxPt="<<maxPt<<
955 "\n";
956 }
957 }
958 }
959}