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