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 "AliESDVertex.h"
  31
  32 #include "AliVEvent.h"
  33 #include "AliVTrack.h"
  34 #include "AliVfriendTrack.h"
  35 #include "AliVfriendEvent.h"
  36
  37 #include "AliMathBase.h"
  38 #include "AliTPCseed.h"
  39 #include "AliTPCclusterMI.h"
  40
  41 #include "AliKFParticle.h"
  42 #include "AliKFVertex.h"
  43
  44 #include "AliTrackPointArray.h"
  45 #include "AliTPCcalibV0.h"
  46 #include "AliV0.h"
  47 #include "TRandom.h"
  48 #include "TTreeStream.h"
  49 #include "AliTPCcalibDB.h"
  50 #include "AliTPCCorrection.h"
  51 #include "AliGRPObject.h"
  52 #include "AliTPCTransform.h"
  53 #include "AliAnalysisManager.h"
  54
  55
  56
  57 ClassImp(AliTPCcalibV0)
  58
  59
  60 AliTPCcalibV0::AliTPCcalibV0() :
  61    AliTPCcalibBase(),
  62    fV0Tree(0),
  63    fHPTTree(0),
  64    fStack(0),
  65    fEvent(0),
  66    fPdg(0),
  67    fParticles(0),
  68    fV0s(0),
  69    fGammas(0)
  70 {
  71
  72 }
  73 AliTPCcalibV0::AliTPCcalibV0(const Text_t *name, const Text_t *title):
  74    AliTPCcalibBase(),
  75    fV0Tree(0),
  76    fHPTTree(0),
  77    fStack(0),
  78    fEvent(0),
  79    fPdg(0),
  80    fParticles(0),
  81    fV0s(0),
  82    fGammas(0)
  83 {
  84   fPdg = new TDatabasePDG;
  85   // create output histograms
  86   SetName(name);
  87   SetTitle(title);
  88 }
  89
  90 AliTPCcalibV0::~AliTPCcalibV0(){
  91   //
  92   //
  93   //
  94   delete fV0Tree;
  95   delete fHPTTree;
  96 }
  97
  98
  99
 100
 101
 102 void  AliTPCcalibV0::ProcessESD(AliVEvent *event){
 103   //
 104   //
 105   //
 106   fEvent = event;
 107   AliKFParticle::SetField(event->GetMagneticField());
 108   if (TMath::Abs(AliTracker::GetBz())<1) return;
 109   DumpToTree(event);
 110   DumpToTreeHPT(event);
 111 }
 112
 113 void  AliTPCcalibV0::DumpToTreeHPT(AliVEvent *event){
 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;
 119   const Float_t kMinPt   =4.;
 120   AliVfriendEvent *friendEvent=event->FindFriend();
 121 //   if (!esdFriend) {
 122 //     Printf("ERROR: esdFriend not available");
 123 //     return;
 124 //   }
 125   //
 126   Int_t ntracks=event->GetNumberOfTracks();
 127   for (Int_t i=0;i<ntracks;++i) {
 128     Bool_t isOK=kFALSE;
 129     AliVTrack *track = event->GetVTrack(i);
 130     if (track->GetTPCncls()<kMinCluster) continue;
 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;
 136       if (!track->IsOn(AliVTrack::kITSrefit)) isAccepted=kFALSE;
 137       if (!track->IsOn(AliVTrack::kTPCrefit)) isAccepted=kFALSE;
 138       if (!track->IsOn(AliVTrack::kTOFout))   isAccepted=kFALSE;
 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;
 147     }
 148     if ( track->GetTPCsignal()>100 && track->GetInnerParam()->Pt()>1 ){
 149       if (track->IsOn(AliVTrack::kITSin)||track->IsOn(AliVTrack::kTRDout)||track->IsOn(AliVTrack::kTOFin))
 150         isOK=kTRUE;
 151       if (isOK){
 152         TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
 153     Int_t eventNumber = event->GetEventNumberInFile();
 154     Bool_t hasFriend=(friendEvent) ? (friendEvent->GetTrack(i)!=0):0;
 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       }
 158     }
 159     if (!isOK) continue;
 160     TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
 161     Int_t eventNumber = event->GetEventNumberInFile();
 162     Bool_t hasFriend=(friendEvent) ? (friendEvent->GetTrack(i)!=0):0;
 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     //
 166     if (!friendEvent) continue;
 167     const AliVfriendTrack *friendTrack = friendEvent->GetTrack(i);
 168     if (!friendTrack) continue;
 169
 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     }
 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
 192     if (fHPTTree->GetEntries()==0){
 193       //
 194       fHPTTree->SetDirectory(0);
 195       fHPTTree->Branch("t.",&dummytrack);
 196       fHPTTree->Branch("ft.",&dummyfriendTrack);
 197       fHPTTree->Branch("s.",&seed);
 198     }else{
 199       fHPTTree->SetBranchAddress("t.",&dummytrack);
 200       fHPTTree->SetBranchAddress("ft.",&dummyfriendTrack);
 201       fHPTTree->SetBranchAddress("s.",&seed);
 202     }
 203     fHPTTree->Fill();
 204     //
 205   }
 206 }
 207
 208
 209
 210 void  AliTPCcalibV0::DumpToTree(AliVEvent *event){
 211   //
 212   // Dump V0s fith full firend information to the
 213   //
 214   Int_t nV0s  = fEvent->GetNumberOfV0s();
 215   const Int_t kMinCluster=110;
 216   const Double_t kDownscale=0.01;
 217   const Float_t kMinPt   =1.0;
 218   const Float_t kMinMinPt   =0.7;
 219   AliVfriendEvent *friendEvent=event->FindFriend();
 220   //
 221
 222   for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
 223     Bool_t isOK=kFALSE;
 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
 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;
 242     //
 243     TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
 244     Int_t eventNumber = event->GetEventNumberInFile();
 245     Bool_t hasITS=(track0->GetNcls(0)+ track1->GetNcls(0)>4);
 246     printf("DUMPHPTV0:%s|%f|%d|%d|%d\n",filename.Data(), (TMath::Min(track0->Pt(),track1->Pt())), eventNumber,(friendEvent!=0), hasITS);
 247     //
 248     if (!friendEvent) continue;
 249     //
 250
 251     //
 252     const AliVfriendTrack *ftrack0 = friendEvent->GetTrack(v0->GetIndex(0));
 253     if (!ftrack0) continue;
 254     const AliVfriendTrack *ftrack1 = friendEvent->GetTrack(v0->GetIndex(1));
 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     }
 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
 288     if (fV0Tree->GetEntries()==0){
 289       //
 290       fV0Tree->SetDirectory(0);
 291       fV0Tree->Branch("v0.",&v0);
 292       fV0Tree->Branch("t0.",&dummytrack0);
 293       fV0Tree->Branch("t1.",&dummytrack1);
 294       fV0Tree->Branch("ft0.",&dummyftrack0);
 295       fV0Tree->Branch("ft1.",&dummyftrack1);
 296       fV0Tree->Branch("s0.",&seed0);
 297       fV0Tree->Branch("s1.",&seed1);
 298     }else{
 299       fV0Tree->SetBranchAddress("v0.",&v0);
 300       fV0Tree->SetBranchAddress("t0.",&dummytrack0);
 301       fV0Tree->SetBranchAddress("t1.",&dummytrack1);
 302       fV0Tree->SetBranchAddress("ft0.",&dummyftrack0);
 303       fV0Tree->SetBranchAddress("ft1.",&dummyftrack1);
 304       fV0Tree->SetBranchAddress("s0.",&seed0);
 305       fV0Tree->SetBranchAddress("s1.",&seed1);
 306     }
 307     fV0Tree->Fill();
 308   }
 309 }
 310
 311
 312 Long64_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
 331 void 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   //
 336   return ;
 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);
 380     //ftrack0->GetCalibContainer()->SetOwner(kTRUE);
 381     //ftrack1->GetCalibContainer()->SetOwner(kTRUE);
 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
 412 void 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   //
 417   return ;
 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
 458 void AliTPCcalibV0::MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t /*run*/){
 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() ;
 487  //  AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
 488 //   AliGRPObject*  grp = AliTPCcalibDB::Instance()->GetGRP(run);
 489 //   Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
 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){
 506         Double_t x0[3]={ static_cast<Double_t>(cluster->GetRow()),cluster->GetPad(),cluster->GetTimeBin()};
 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     //
 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 }
 572
 573 void 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
 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() ;
 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();
 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;
 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;
 637     rmass0=massPion;
 638     rmass1=massPion;
 639     pdg0=pdgPion;
 640     pdg1=pdgPion;
 641     if (isLambda) {rmass0=massProton; pdg0=pdgProton;}
 642     if (isAntiLambda) {rmass1=massProton; pdg1=pdgProton;}
 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){
 660           Double_t x0[3]={ static_cast<Double_t>(cluster->GetRow()),cluster->GetPad(),cluster->GetTimeBin()};
 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       //
 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);
 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         //
 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;
 704         if (!isOK) continue;
 705         arrayT0.AddAt(trackIO0->Clone(),icorr+1);
 706         arrayT1.AddAt(trackIO1->Clone(),icorr+1);
 707         Int_t charge=TMath::Nint(trackIO0->GetSign());
 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
 738     Int_t id=Int_t(1000.*(param0->Pt()-param1->Pt()));      // K0s V0 asymetry
 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);
 784       //      AliKFParticle* pin= (AliKFParticle*)arrayV0in.At(icorr+1);
 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 }
 803
 804 AliExternalTrackParam * AliTPCcalibV0::RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass){
 805   //
 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;
 829   //
 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
 874
 875
 876 AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){
 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
 891
 892
 893
 894 void 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);
 912    delete [] new_bins;
 913 }
 914
 915
 916
 917 void AliTPCcalibV0::FilterV0s(AliVEvent *event){
 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   //
 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;
 938   AliKFVertex kfvertex=*vertex;
 939
 940   //
 941   for (Int_t iv0=0;iv0<nv0;iv0++){
 942     AliESDv0 dummyv0;
 943     event->GetV0(dummyv0,iv0);
 944     AliESDv0 *v0=&dummyv0;
 945
 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){
 979       AliVTrack * trackP = event->GetVTrack(v0->GetPindex());
 980       AliVTrack * trackN = event->GetVTrack(v0->GetNindex());
 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           //
 995       //"v.="<<vertex<<
 996           "ncls="<<ncls<<
 997           "maxPt="<<maxPt<<
 998           "\n";
 999       }
1000     }
1001   }
1002 }