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