TPC/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 "TCint.h"
  39 #include "AliTPCcalibV0.h"
  40 #include "AliV0.h"
  41 #include "TRandom.h"
  42 #include "TTreeStream.h"
  43 #include "AliTPCcalibDB.h"
  44 #include "AliTPCCorrection.h"
  45 #include "AliGRPObject.h"
  46 #include "AliTPCTransform.h"
  47 #include "AliAnalysisManager.h"
  48
  49
  50
  51 ClassImp(AliTPCcalibV0)
  52
  53
  54 AliTPCcalibV0::AliTPCcalibV0() :
  55    AliTPCcalibBase(),
  56    fV0Tree(0),
  57    fHPTTree(0),
  58    fStack(0),
  59    fESD(0),
  60    fPdg(0),
  61    fParticles(0),
  62    fV0s(0),
  63    fGammas(0)
  64 {
  65
  66 }
  67 AliTPCcalibV0::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)
  77 {
  78   fPdg = new TDatabasePDG;
  79   // create output histograms
  80   SetName(name);
  81   SetTitle(title);
  82 }
  83
  84 AliTPCcalibV0::~AliTPCcalibV0(){
  85   //
  86   //
  87   //
  88   delete fV0Tree;
  89   delete fHPTTree;
  90 }
  91
  92
  93
  94
  95
  96 void  AliTPCcalibV0::ProcessESD(AliESDEvent *esd){
  97   //
  98   //
  99   //
 100   fESD = esd;
 101   AliKFParticle::SetField(esd->GetMagneticField());
 102   if (TMath::Abs(AliTracker::GetBz())<1) return;
 103   DumpToTree(esd);
 104   DumpToTreeHPT(esd);
 105 }
 106
 107 void  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;
 113   const Float_t kMinPt   =4.;
 114   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
 115 //   if (!esdFriend) {
 116 //     Printf("ERROR: esdFriend not available");
 117 //     return;
 118 //   }
 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;
 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;
 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       }
 152     }
 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
 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
 195 void  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;
 202   const Float_t kMinPt   =1.0;
 203   const Float_t kMinMinPt   =0.7;
 204   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
 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;
 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
 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
 284 Long64_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
 303 void 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   //
 308   return ;
 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);
 352     //ftrack0->GetCalibContainer()->SetOwner(kTRUE);
 353     //ftrack1->GetCalibContainer()->SetOwner(kTRUE);
 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
 384 void 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   //
 389   return ;
 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
 430 void AliTPCcalibV0::MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t /*run*/){
 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() ;
 459  //  AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
 460 //   AliGRPObject*  grp = AliTPCcalibDB::Instance()->GetGRP(run);
 461 //   Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
 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     //
 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 }
 544
 545 void 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
 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() ;
 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();
 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;
 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;
 609     rmass0=massPion;
 610     rmass1=massPion;
 611     pdg0=pdgPion;
 612     pdg1=pdgPion;
 613     if (isLambda) {rmass0=massProton; pdg0=pdgProton;}
 614     if (isAntiLambda) {rmass1=massProton; pdg1=pdgProton;}
 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       //
 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);
 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         //
 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;
 676         if (!isOK) continue;
 677         arrayT0.AddAt(trackIO0->Clone(),icorr+1);
 678         arrayT1.AddAt(trackIO1->Clone(),icorr+1);
 679         Int_t charge=TMath::Nint(trackIO0->GetSign());
 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
 710     Int_t id=Int_t(1000.*(param0->Pt()-param1->Pt()));      // K0s V0 asymetry
 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);
 756       //      AliKFParticle* pin= (AliKFParticle*)arrayV0in.At(icorr+1);
 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 }
 775
 776 AliExternalTrackParam * AliTPCcalibV0::RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass){
 777   //
 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;
 801   //
 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
 846
 847
 848 AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){
 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
 863
 864
 865
 866 void 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);
 884    delete [] new_bins;
 885 }
 886
 887
 888
 889 void 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 }