TPC/Calib/AliTPCcalibV0.cxx

   1
   2 #include <TROOT.h>
   3 #include <TChain.h>
   4 #include <TFile.h>
   5 #include <TH3F.h>
   6 #include <TH2F.h>
   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"
  26 #include "AliESDEvent.h"
  27 #include "AliESDtrack.h"
  28 #include "AliESDfriend.h"
  29 #include "AliESDfriendTrack.h"
  30 #include "AliMathBase.h"
  31 #include "AliTPCseed.h"
  32 #include "AliTPCclusterMI.h"
  33
  34 #include "AliKFParticle.h"
  35 #include "AliKFVertex.h"
  36
  37 #include "AliTrackPointArray.h"
  38 #include "AliTPCcalibV0.h"
  39 #include "AliV0.h"
  40 #include "TRandom.h"
  41 #include "TTreeStream.h"
  42 #include "AliTPCcalibDB.h"
  43 #include "AliTPCCorrection.h"
  44 #include "AliGRPObject.h"
  45 #include "AliTPCTransform.h"
  46 #include "AliAnalysisManager.h"
  47
  48
  49
  50 ClassImp(AliTPCcalibV0)
  51
  52
  53 AliTPCcalibV0::AliTPCcalibV0() :
  54    AliTPCcalibBase(),
  55    fV0Tree(0),
  56    fHPTTree(0),
  57    fStack(0),
  58    fESD(0),
  59    fPdg(0),
  60    fParticles(0),
  61    fV0s(0),
  62    fGammas(0)
  63 {
  64
  65 }
  66 AliTPCcalibV0::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)
  76 {
  77   fPdg = new TDatabasePDG;
  78   // create output histograms
  79   SetName(name);
  80   SetTitle(title);
  81 }
  82
  83 AliTPCcalibV0::~AliTPCcalibV0(){
  84   //
  85   //
  86   //
  87   delete fV0Tree;
  88   delete fHPTTree;
  89 }
  90
  91
  92
  93
  94
  95 void  AliTPCcalibV0::ProcessESD(AliESDEvent *esd){
  96   //
  97   //
  98   //
  99   fESD = esd;
 100   AliKFParticle::SetField(esd->GetMagneticField());
 101   if (TMath::Abs(AliTracker::GetBz())<1) return;
 102   DumpToTree(esd);
 103   DumpToTreeHPT(esd);
 104 }
 105
 106 void  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;
 112   const Float_t kMinPt   =4.;
 113   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
 114 //   if (!esdFriend) {
 115 //     Printf("ERROR: esdFriend not available");
 116 //     return;
 117 //   }
 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;
 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;
 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       }
 151     }
 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
 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
 194 void  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;
 201   const Float_t kMinPt   =1.0;
 202   const Float_t kMinMinPt   =0.7;
 203   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
 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;
 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
 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
 283 Long64_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
 302 void 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   //
 307   return ;
 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);
 351     //ftrack0->GetCalibContainer()->SetOwner(kTRUE);
 352     //ftrack1->GetCalibContainer()->SetOwner(kTRUE);
 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
 383 void 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   //
 388   return ;
 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
 429 void AliTPCcalibV0::MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t /*run*/){
 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() ;
 458  //  AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
 459 //   AliGRPObject*  grp = AliTPCcalibDB::Instance()->GetGRP(run);
 460 //   Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
 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     //
 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 }
 543
 544 void 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
 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() ;
 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();
 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;
 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;
 608     rmass0=massPion;
 609     rmass1=massPion;
 610     pdg0=pdgPion;
 611     pdg1=pdgPion;
 612     if (isLambda) {rmass0=massProton; pdg0=pdgProton;}
 613     if (isAntiLambda) {rmass1=massProton; pdg1=pdgProton;}
 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       //
 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);
 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         //
 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;
 675         if (!isOK) continue;
 676         arrayT0.AddAt(trackIO0->Clone(),icorr+1);
 677         arrayT1.AddAt(trackIO1->Clone(),icorr+1);
 678         Int_t charge=TMath::Nint(trackIO0->GetSign());
 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
 709     Int_t id=Int_t(1000.*(param0->Pt()-param1->Pt()));      // K0s V0 asymetry
 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);
 755       //      AliKFParticle* pin= (AliKFParticle*)arrayV0in.At(icorr+1);
 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 }
 774
 775 AliExternalTrackParam * AliTPCcalibV0::RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass){
 776   //
 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;
 800   //
 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
 845
 846
 847 AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){
 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
 862
 863
 864
 865 void 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);
 883    delete [] new_bins;
 884 }
 885
 886
 887
 888 void 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 }