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
  48
  49
  50
  51
  52 ClassImp(AliTPCcalibV0)
  53
  54
  55 AliTPCcalibV0::AliTPCcalibV0() :
  56    AliTPCcalibBase(),
  57    fV0Tree(0),
  58    fHPTTree(0),
  59    fStack(0),
  60    fESD(0),
  61    fPdg(0),
  62    fParticles(0),
  63    fV0s(0),
  64    fGammas(0)
  65 {
  66
  67 }
  68 AliTPCcalibV0::AliTPCcalibV0(const Text_t *name, const Text_t *title):
  69    AliTPCcalibBase(),
  70    fV0Tree(0),
  71    fHPTTree(0),
  72    fStack(0),
  73    fESD(0),
  74    fPdg(0),
  75    fParticles(0),
  76    fV0s(0),
  77    fGammas(0)
  78 {
  79   fPdg = new TDatabasePDG;
  80   // create output histograms
  81   SetName(name);
  82   SetTitle(title);
  83 }
  84
  85 AliTPCcalibV0::~AliTPCcalibV0(){
  86   //
  87   //
  88   //
  89   delete fV0Tree;
  90   delete fHPTTree;
  91 }
  92
  93
  94
  95
  96
  97 void  AliTPCcalibV0::ProcessESD(AliESDEvent *esd, AliStack *stack){
  98   //
  99   //
 100   //
 101   fESD = esd;
 102   AliKFParticle::SetField(esd->GetMagneticField());
 103   if (TMath::Abs(AliTracker::GetBz())<1) return;
 104   DumpToTree(esd);
 105   DumpToTreeHPT(esd);
 106   return;
 107   //
 108   MakeV0s();
 109   if (stack) {
 110     fStack = stack;
 111     MakeMC();
 112   }else{
 113     fStack =0;
 114   }
 115 }
 116
 117 void  AliTPCcalibV0::DumpToTreeHPT(AliESDEvent *esd){
 118   //
 119   // Dump V0s fith full firend information to the
 120   //
 121   if (TMath::Abs(AliTracker::GetBz())<1) return;
 122   const Int_t kMinCluster=110;
 123   const Float_t kMinPt   =3.;
 124   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
 125   if (!esdFriend) {
 126     Printf("ERROR: esdFriend not available");
 127     return;
 128   }
 129   //
 130   Int_t ntracks=esd->GetNumberOfTracks();
 131   for (Int_t i=0;i<ntracks;++i) {
 132     Bool_t isOK=kFALSE;
 133     AliESDtrack *track = esd->GetTrack(i);
 134     if (track->GetTPCncls()<kMinCluster) continue;
 135     AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
 136     if (!friendTrack) continue;
 137     if (TMath::Abs(AliTracker::GetBz())>1){ // cut on momenta if measured
 138       if (track->Pt()>kMinPt) isOK=kTRUE;
 139     }
 140     if (TMath::Abs(AliTracker::GetBz())<1){  // require primary track for the B field OFF data
 141       Bool_t isAccepted=kTRUE;
 142       if (!track->IsOn(AliESDtrack::kITSrefit)) isAccepted=kFALSE;
 143       if (!track->IsOn(AliESDtrack::kTPCrefit)) isAccepted=kFALSE;
 144       if (!track->IsOn(AliESDtrack::kTOFout))   isAccepted=kFALSE;
 145       Float_t dvertex[2],cvertex[3];
 146       track->GetImpactParametersTPC(dvertex,cvertex);
 147       if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>20) isAccepted=kFALSE;
 148       if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>20) isAccepted=kFALSE;
 149       track->GetImpactParameters(dvertex,cvertex);
 150       if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>10) isAccepted=kFALSE;
 151       if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>10) isAccepted=kFALSE;
 152       if (!isAccepted) isOK=kFALSE;
 153     }
 154
 155     if (!isOK) continue;
 156     //
 157
 158     TObject *calibObject;
 159     AliTPCseed *seed = 0;
 160     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 161       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 162     }
 163     if (!seed) continue;
 164       if (!fHPTTree) {
 165       fHPTTree = new TTree("HPT","HPT");
 166       fHPTTree->SetDirectory(0);
 167     }
 168     if (fHPTTree->GetEntries()==0){
 169       //
 170       fHPTTree->SetDirectory(0);
 171       fHPTTree->Branch("t.",&track);
 172       fHPTTree->Branch("ft.",&friendTrack);
 173       fHPTTree->Branch("s.",&seed);
 174     }else{
 175       fHPTTree->SetBranchAddress("t.",&track);
 176       fHPTTree->SetBranchAddress("ft.",&friendTrack);
 177       fHPTTree->SetBranchAddress("s.",&seed);
 178     }
 179     fHPTTree->Fill();
 180     //
 181   }
 182 }
 183
 184
 185
 186 void  AliTPCcalibV0::DumpToTree(AliESDEvent *esd){
 187   //
 188   // Dump V0s fith full firend information to the
 189   //
 190   Int_t nV0s  = fESD->GetNumberOfV0s();
 191   const Int_t kMinCluster=110;
 192   const Double_t kDownscale=0.01;
 193   const Float_t kMinPt   =1.0;
 194   const Float_t kMinMinPt   =0.7;
 195   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
 196   if (!esdFriend) {
 197     Printf("ERROR: esdFriend not available");
 198     return;
 199   }
 200   //
 201
 202   for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
 203     Bool_t isOK=kFALSE;
 204     AliESDv0 * v0 = (AliESDv0*) esd->GetV0(ivertex);
 205     AliESDtrack * track0 = fESD->GetTrack(v0->GetIndex(0)); // negative track
 206     AliESDtrack * track1 = fESD->GetTrack(v0->GetIndex(1)); // positive track
 207     if (track0->GetTPCNcls()<kMinCluster) continue;
 208     if (track0->GetKinkIndex(0)>0) continue;
 209     if (track1->GetTPCNcls()<kMinCluster) continue;
 210     if (track1->GetKinkIndex(0)>0) continue;
 211     if (v0->GetOnFlyStatus()==kFALSE) continue;
 212     //
 213     if (TMath::Min(track0->Pt(),track1->Pt())<kMinMinPt) continue;
 214     //
 215     //
 216     if (TMath::Max(track0->Pt(),track1->Pt())>kMinPt) isOK=kTRUE;
 217     if (gRandom->Rndm()<kDownscale) isOK=kTRUE;
 218     if (!isOK) continue;
 219     //
 220     AliESDfriendTrack *ftrack0 = esdFriend->GetTrack(v0->GetIndex(0));
 221     if (!ftrack0) continue;
 222     AliESDfriendTrack *ftrack1 = esdFriend->GetTrack(v0->GetIndex(1));
 223     if (!ftrack1) continue;
 224     //
 225     TObject *calibObject;
 226     AliTPCseed *seed0 = 0;
 227     AliTPCseed *seed1 = 0;
 228     for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
 229       if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
 230     }
 231     for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
 232       if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
 233     }
 234     if (!seed0) continue;
 235     if (!seed1) continue;
 236     AliExternalTrackParam * paramIn0 = (AliExternalTrackParam *)track0->GetInnerParam();
 237     AliExternalTrackParam * paramIn1 = (AliExternalTrackParam *)track1->GetInnerParam();
 238     if (!paramIn0) continue;
 239     if (!paramIn1) continue;
 240     //
 241     //
 242     if (!fV0Tree) {
 243       fV0Tree = new TTree("V0s","V0s");
 244       fV0Tree->SetDirectory(0);
 245     }
 246     if (fV0Tree->GetEntries()==0){
 247       //
 248       fV0Tree->SetDirectory(0);
 249       fV0Tree->Branch("v0.",&v0);
 250       fV0Tree->Branch("t0.",&track0);
 251       fV0Tree->Branch("t1.",&track1);
 252       fV0Tree->Branch("ft0.",&ftrack0);
 253       fV0Tree->Branch("ft1.",&ftrack1);
 254       fV0Tree->Branch("s0.",&seed0);
 255       fV0Tree->Branch("s1.",&seed1);
 256     }else{
 257       fV0Tree->SetBranchAddress("v0.",&v0);
 258       fV0Tree->SetBranchAddress("t0.",&track0);
 259       fV0Tree->SetBranchAddress("t1.",&track1);
 260       fV0Tree->SetBranchAddress("ft0.",&ftrack0);
 261       fV0Tree->SetBranchAddress("ft1.",&ftrack1);
 262       fV0Tree->SetBranchAddress("s0.",&seed0);
 263       fV0Tree->SetBranchAddress("s1.",&seed1);
 264     }
 265     fV0Tree->Fill();
 266   }
 267 }
 268
 269
 270 Long64_t AliTPCcalibV0::Merge(TCollection *const li) {
 271
 272   TIterator* iter = li->MakeIterator();
 273   AliTPCcalibV0* cal = 0;
 274
 275   while ((cal = (AliTPCcalibV0*)iter->Next())) {
 276     if (cal->fV0Tree){
 277       if (!fV0Tree) {
 278         fV0Tree = new TTree("V0s","V0s");
 279         fV0Tree->SetDirectory(0);
 280       }
 281       if (cal->fV0Tree->GetEntries()>0) AliTPCcalibV0::AddTree(cal->fV0Tree);
 282       if (cal->fHPTTree->GetEntries()>0) AliTPCcalibV0::AddTreeHPT(cal->fHPTTree);
 283     }
 284   }
 285   return 0;
 286 }
 287
 288
 289 void AliTPCcalibV0::AddTree(TTree * treeInput){
 290   //
 291   // Add the content of tree:
 292   // Notice automatic copy of tree in ROOT does not work for such complicated tree
 293   //
 294   AliESDv0 * v0 = new AliESDv0;
 295   Double_t kMinPt=0.8;
 296   AliESDtrack * track0 = 0; // negative track
 297   AliESDtrack * track1 = 0; // positive track
 298   AliESDfriendTrack *ftrack0 = 0;
 299   AliESDfriendTrack *ftrack1 = 0;
 300   AliTPCseed *seed0 = 0;
 301   AliTPCseed *seed1 = 0;
 302   treeInput->SetBranchStatus("ft0.",kFALSE);
 303   treeInput->SetBranchStatus("ft1.",kFALSE);
 304   TDatabasePDG pdg;
 305   Double_t massK0= pdg.GetParticle("K0")->Mass();
 306   Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
 307
 308   Int_t entries= treeInput->GetEntries();
 309   for (Int_t i=0; i<entries; i++){
 310     treeInput->SetBranchAddress("v0.",&v0);
 311     treeInput->SetBranchAddress("t0.",&track0);
 312     treeInput->SetBranchAddress("t1.",&track1);
 313     treeInput->SetBranchAddress("ft0.",&ftrack0);
 314     treeInput->SetBranchAddress("ft1.",&ftrack1);
 315     treeInput->SetBranchAddress("s0.",&seed0);
 316     treeInput->SetBranchAddress("s1.",&seed1);
 317     if (fV0Tree->GetEntries()==0){
 318       fV0Tree->SetDirectory(0);
 319       fV0Tree->Branch("v0.",&v0);
 320       fV0Tree->Branch("t0.",&track0);
 321       fV0Tree->Branch("t1.",&track1);
 322       fV0Tree->Branch("ft0.",&ftrack0);
 323       fV0Tree->Branch("ft1.",&ftrack1);
 324       fV0Tree->Branch("s0.",&seed0);
 325       fV0Tree->Branch("s1.",&seed1);
 326     }else{
 327       fV0Tree->SetBranchAddress("v0.",&v0);
 328       fV0Tree->SetBranchAddress("t0.",&track0);
 329       fV0Tree->SetBranchAddress("t1.",&track1);
 330       fV0Tree->SetBranchAddress("ft0.",&ftrack0);
 331       fV0Tree->SetBranchAddress("ft1.",&ftrack1);
 332       fV0Tree->SetBranchAddress("s0.",&seed0);
 333       fV0Tree->SetBranchAddress("s1.",&seed1);
 334     }
 335     //
 336     treeInput->GetEntry(i);
 337     //ftrack0->GetCalibContainer()->SetOwner(kTRUE);
 338     //ftrack1->GetCalibContainer()->SetOwner(kTRUE);
 339     Bool_t isOK=kTRUE;
 340     if (v0->GetOnFlyStatus()==kFALSE) isOK=kFALSE;
 341     if (track0->GetTPCncls()<100) isOK=kFALSE;
 342     if (track1->GetTPCncls()<100) isOK=kFALSE;
 343     if (TMath::Min(seed0->Pt(),seed1->Pt())<kMinPt) isOK=kFALSE;
 344     if (TMath::Min(track0->Pt(),track1->Pt())<kMinPt) isOK=kFALSE;
 345     Bool_t isV0=kFALSE;
 346     if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.05)     isV0=kTRUE;
 347     if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.05) isV0=kTRUE;
 348     if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.05) isV0=kTRUE;
 349     if (TMath::Abs(v0->GetEffMass(0,0))<0.02) isV0=kFALSE; //reject electrons
 350     if (!isV0) isOK=kFALSE;
 351     if (isOK) fV0Tree->Fill();
 352     delete v0;
 353     delete track0;
 354     delete track1;
 355     delete ftrack0;
 356     delete ftrack1;
 357     delete seed0;
 358     delete seed1;
 359     v0=0;
 360     track0=0;
 361     track1=0;
 362     ftrack0=0;
 363     ftrack1=0;
 364     seed0=0;
 365     seed1=0;
 366   }
 367 }
 368
 369 void AliTPCcalibV0::AddTreeHPT(TTree * treeInput){
 370   //
 371   // Add the content of tree:
 372   // Notice automatic copy of tree in ROOT does not work for such complicated tree
 373   //
 374   AliESDtrack *track = 0;
 375   AliESDfriendTrack *friendTrack = 0;
 376   AliTPCseed *seed = 0;
 377   if (!treeInput) return;
 378   if (treeInput->GetEntries()==0) return;
 379   //
 380   Int_t entries= treeInput->GetEntries();
 381   //
 382   for (Int_t i=0; i<entries; i++){
 383     track=0;
 384     friendTrack=0;
 385     seed=0;
 386     //
 387     treeInput->SetBranchAddress("t.",&track);
 388     treeInput->SetBranchAddress("ft.",&friendTrack);
 389     treeInput->SetBranchAddress("s.",&seed);
 390     treeInput->GetEntry(i);
 391     //
 392     if (fHPTTree->GetEntries()==0){
 393       fHPTTree->SetDirectory(0);
 394       fHPTTree->Branch("t.",&track);
 395       fHPTTree->Branch("ft.",&friendTrack);
 396       fHPTTree->Branch("s.",&seed);
 397     }else{
 398       fHPTTree->SetBranchAddress("t.",&track);
 399       fHPTTree->SetBranchAddress("ft.",&friendTrack);
 400       fHPTTree->SetBranchAddress("s.",&seed);
 401     }
 402     Bool_t isOK=kTRUE;
 403     if (!track->IsOn(AliESDtrack::kITSrefit)) isOK=kFALSE;
 404     if (!track->IsOn(AliESDtrack::kTOFout)) isOK=kFALSE;
 405     if (isOK) fHPTTree->Fill();
 406     //
 407     delete track;
 408     delete friendTrack;
 409     delete seed;
 410   }
 411 }
 412
 413
 414 void AliTPCcalibV0::MakeMC(){
 415   //
 416   // MC comparison
 417   //    1. Select interesting particles
 418   //    2. Assign the recosntructed particles
 419   //
 420   //1. Select interesting particles
 421   const Float_t kMinP   = 0.2;
 422   const Float_t kMinPt  = 0.1;
 423   const Float_t kMaxR   = 0.5;
 424   const Float_t kMaxTan = 1.2;
 425   const Float_t kMaxRad = 150;
 426   //
 427   if (!fParticles) fParticles = new TObjArray;
 428   TParticle *part=0;
 429   //
 430   Int_t entries = fStack->GetNtrack();
 431   for (Int_t ipart=0; ipart<entries; ipart++){
 432     part = fStack->Particle(ipart);
 433     if (!part) continue;
 434     if (part->P()<kMinP) continue;
 435     if (part->R()>kMaxR) continue;
 436     if (TMath::Abs(TMath::Tan(part->Theta()-TMath::Pi()*0.5))>kMaxTan) continue;
 437     Bool_t isInteresting = kFALSE;
 438     if (part->GetPdgCode()==22) isInteresting =kTRUE;
 439     if (part->GetPdgCode()==310) isInteresting =kTRUE;
 440     if (part->GetPdgCode()==111) isInteresting =kTRUE;
 441     if (TMath::Abs(part->GetPdgCode()==3122)) isInteresting =kTRUE;
 442
 443     //
 444     if (!isInteresting) continue;
 445     fParticles->AddLast(new TParticle(*part));
 446   }
 447   if (fParticles->GetEntries()<1) {
 448     return;
 449   }
 450   //
 451   //
 452   //
 453   Int_t sentries=fParticles->GetEntries();;
 454   for (Int_t ipart=0; ipart<sentries; ipart++){
 455     part = (TParticle*)fParticles->At(ipart);
 456     TParticle *p0 = 0;
 457     TParticle *p1 = 0;
 458
 459     Int_t nold =0;
 460     Int_t nnew =0;
 461     Int_t id0  = part->GetDaughter(0);
 462     Int_t id1  = part->GetDaughter(1);
 463     if (id0>=fStack->GetNtrack() ) id0*=-1;
 464     if (id1>=fStack->GetNtrack() ) id1*=-1;
 465     Bool_t findable = kTRUE;
 466     if (id0<0 || id1<0) findable = kFALSE;
 467     Int_t charge =0;
 468     if (findable){
 469       p0 = fStack->Particle(id0);
 470       if (p0->R()>kMaxRad) findable = kFALSE;
 471       if (p0->Pt()<kMinPt) findable = kFALSE;
 472       if (p0->Vz()>250) findable= kFALSE;
 473       if (TMath::Abs(TMath::Tan(p0->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE;
 474       if (fPdg->GetParticle(p0->GetPdgCode())==0) findable =kFALSE;
 475       else
 476         if (fPdg->GetParticle(p0->GetPdgCode())->Charge()==0) charge++;
 477
 478       p1 = fStack->Particle(id1);
 479       if (p1->R()>kMaxRad) findable = kFALSE;
 480       if (p1->Pt()<kMinPt) findable = kFALSE;
 481       if (TMath::Abs(p1->Vz())>250) findable= kFALSE;
 482       if (TMath::Abs(TMath::Tan(p1->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE;
 483       if (fPdg->GetParticle(p1->GetPdgCode())==0) findable = kFALSE;
 484       else
 485         if (fPdg->GetParticle(p1->GetPdgCode())->Charge()==0) charge++;
 486
 487     }
 488     //   (*fDebugStream)<<"MC0"<<
 489     //       "P.="<<part<<
 490     //       "findable="<<findable<<
 491     //       "id0="<<id0<<
 492     //       "id1="<<id1<<
 493     //       "\n";
 494     if (!findable) continue;
 495     Float_t minpt = TMath::Min(p0->Pt(), p1->Pt());
 496     Int_t type=-1;
 497
 498     //
 499     //
 500     AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
 501     for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
 502       AliESDv0 * v0 = fESD->GetV0(ivertex);
 503       // select coresponding track
 504       AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0));
 505       if (TMath::Abs(trackN->GetLabel())!=id0 && TMath::Abs(trackN->GetLabel())!=id1) continue;
 506       AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1));
 507       if (TMath::Abs(trackP->GetLabel())!=id0 && TMath::Abs(trackP->GetLabel())!=id1) continue;
 508       TParticle *pn = fStack->Particle(TMath::Abs(trackN->GetLabel()));
 509       TParticle *pp = fStack->Particle(TMath::Abs(trackP->GetLabel()));
 510       //
 511       //
 512       if ( v0->GetOnFlyStatus()) nnew++;
 513       if (!v0->GetOnFlyStatus()) nold++;
 514       if (part->GetPdgCode()==22 && TMath::Abs(pn->GetPdgCode())==11 &&  TMath::Abs(pp->GetPdgCode())==11)
 515         type =1;
 516       if (part->GetPdgCode()==310 && TMath::Abs(pn->GetPdgCode())==211 &&  TMath::Abs(pp->GetPdgCode())==211)
 517         type =0;
 518       if (part->GetPdgCode()==3122){
 519         if (TMath::Abs(pn->GetPdgCode())==210 ) type=2;
 520         else type=3;
 521       }
 522       AliKFParticle *v0kf       = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode());
 523       v0kf->SetProductionVertex( primVtx );
 524       AliKFParticle *v0kfc       = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode());
 525       v0kfc->SetProductionVertex( primVtx );
 526       v0kfc->SetMassConstraint(fPdg->GetParticle(part->GetPdgCode())->Mass());
 527       Float_t chi2 = v0kf->GetChi2();
 528       Float_t chi2C = v0kf->GetChi2();
 529       //
 530       //
 531       TTreeSRedirector *cstream = GetDebugStreamer();
 532       if (cstream){
 533         (*cstream)<<"MCRC"<<
 534           "P.="<<part<<
 535           "type="<<type<<
 536           "chi2="<<chi2<<
 537           "chi2C="<<chi2C<<
 538           "minpt="<<minpt<<
 539           "id0="<<id0<<
 540           "id1="<<id1<<
 541           "Pn.="<<pn<<
 542           "Pp.="<<pp<<
 543           "tn.="<<trackN<<
 544           "tp.="<<trackP<<
 545           "nold.="<<nold<<
 546           "nnew.="<<nnew<<
 547           "v0.="<<v0<<
 548           "v0kf.="<<v0kf<<
 549           "v0kfc.="<<v0kfc<<
 550           "\n";
 551         delete v0kf;
 552         delete v0kfc;
 553         //
 554       }
 555
 556       if (cstream){
 557         (*cstream)<<"MC"<<
 558           "P.="<<part<<
 559           "charge="<<charge<<
 560           "type="<<type<<
 561           "minpt="<<minpt<<
 562           "id0="<<id0<<
 563           "id1="<<id1<<
 564           "P0.="<<p0<<
 565           "P1.="<<p1<<
 566           "nold="<<nold<<
 567           "nnew="<<nnew<<
 568           "\n";
 569       }
 570     }
 571     fParticles->Delete();
 572   }
 573 }
 574
 575 void AliTPCcalibV0::MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t /*run*/){
 576   //
 577   // Make a fit tree
 578   //
 579   // 0. Loop over selected tracks
 580   // 1. Loop over all transformation - refit the track with and without the
 581   //    transformtation
 582   // 2. Dump the matching paramaeters to the debugStremer
 583   //
 584
 585   //Connect input
 586   const Int_t kMinNcl=120;
 587   TFile f("TPCV0Objects.root");
 588   AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
 589   TTree * treeInput = v0TPC->GetHPTTree();
 590   TTreeSRedirector *pcstream = new TTreeSRedirector("fitHPT.root");
 591   AliESDtrack *track = 0;
 592   AliESDfriendTrack *friendTrack = 0;
 593   AliTPCseed *seed = 0;
 594   if (!treeInput) return;
 595   if (treeInput->GetEntries()==0) return;
 596   //
 597   treeInput->SetBranchAddress("t.",&track);
 598   treeInput->SetBranchAddress("ft.",&friendTrack);
 599   treeInput->SetBranchAddress("s.",&seed);
 600   //
 601   Int_t ncorr=0;
 602   if (corrArray) ncorr = corrArray->GetEntries();
 603   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
 604  //  AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
 605 //   AliGRPObject*  grp = AliTPCcalibDB::Instance()->GetGRP(run);
 606 //   Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
 607   //
 608   //
 609   //
 610   Int_t ntracks= treeInput->GetEntries();
 611   for (Int_t itrack=0; itrack<ntracks; itrack++){
 612     treeInput->GetEntry(itrack);
 613     if (!track) continue;
 614     if (seed->Pt()<ptCut) continue;
 615     if (track->Pt()<ptCut) continue;
 616     if (track->GetTPCncls()<kMinNcl) continue;
 617     //
 618     // Reapply transformation
 619     //
 620     for (Int_t irow=0; irow<159; irow++){
 621       AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
 622       if (cluster &&cluster->GetX()>10){
 623         Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
 624         Int_t index0[1]={cluster->GetDetector()};
 625         transform->Transform(x0,index0,0,1);
 626         cluster->SetX(x0[0]);
 627         cluster->SetY(x0[1]);
 628         cluster->SetZ(x0[2]);
 629         //
 630       }
 631     }
 632     //
 633     AliExternalTrackParam* paramInner=0;
 634     AliExternalTrackParam* paramOuter=0;
 635     AliExternalTrackParam* paramIO=0;
 636     Bool_t isOK=kTRUE;
 637     for (Int_t icorr=-1; icorr<ncorr; icorr++){
 638       //
 639       AliTPCCorrection *corr = 0;
 640       if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
 641       AliExternalTrackParam * trackInner = RefitTrack(seed, corr,85,134,0.1);
 642       AliExternalTrackParam * trackIO = RefitTrack(seed, corr,245,85,0.1);
 643       AliExternalTrackParam * trackOuter = RefitTrack(seed, corr,245,134,0.1 );
 644       if (trackInner&&trackOuter&&trackIO){
 645         trackOuter->Rotate(trackInner->GetAlpha());
 646         trackOuter->PropagateTo(trackInner->GetX(),AliTracker::GetBz());
 647         if (icorr<0) {
 648           paramInner=trackInner;
 649           paramOuter=trackOuter;
 650           paramIO=trackIO;
 651           paramIO->Rotate(seed->GetAlpha());
 652           paramIO->PropagateTo(seed->GetX(),AliTracker::GetBz());
 653         }
 654       }else{
 655         isOK=kFALSE;
 656       }
 657
 658     }
 659     if (paramOuter&& paramInner) {
 660       //      Bool_t isOK=kTRUE;
 661       if (paramInner->GetSigmaY2()>0.01) isOK&=kFALSE;
 662       if (paramOuter->GetSigmaY2()>0.01) isOK&=kFALSE;
 663       if (paramInner->GetSigmaZ2()>0.01) isOK&=kFALSE;
 664       if (paramOuter->GetSigmaZ2()>0.01) isOK&=kFALSE;
 665       (*pcstream)<<"fit"<<
 666         "s.="<<seed<<
 667         "io.="<<paramIO<<
 668         "pIn.="<<paramInner<<
 669         "pOut.="<<paramOuter;
 670     }
 671     //
 672   }
 673   delete pcstream;
 674   /*
 675     .x ~/rootlogon.C
 676     Int_t run=117112;
 677     .x ../ConfigCalibTrain.C(run)
 678     .L ../AddTaskTPCCalib.C
 679     ConfigOCDB(run)
 680     TFile fexb("../../RegisterCorrectionExB.root");
 681     AliTPCComposedCorrection *cc=  (AliTPCComposedCorrection*) fexb.Get("ComposedExB");
 682     cc->Init();
 683     cc->Print("DA"); // Print used correction classes
 684     TObjArray *array = cc->GetCorrections()
 685     AliTPCcalibV0::MakeFitTreeTrack(array,2,run);
 686
 687    */
 688 }
 689
 690 void AliTPCcalibV0::MakeFitTreeV0(const TObjArray * corrArray, Double_t ptCut, Int_t run){
 691   //
 692   // Make a fit tree
 693   //
 694   // 0. Loop over selected tracks
 695   // 1. Loop over all transformation - refit the track with and without the
 696   //    transformtation
 697   // 2. Dump the matching paramaeters to the debugStremer
 698   //
 699
 700   //Connect input
 701   TFile f("TPCV0Objects.root");
 702   AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
 703   TTree * treeInput = v0TPC->GetV0Tree();
 704   TTreeSRedirector *pcstream = new TTreeSRedirector("fitV0.root");
 705   AliESDv0 *v0 = 0;
 706   AliESDtrack *track0 = 0;
 707   AliESDfriendTrack *friendTrack0 = 0;
 708   AliTPCseed *seed0 = 0;
 709   AliTPCseed *s0 = 0;
 710   AliESDtrack *track1 = 0;
 711   AliESDfriendTrack *friendTrack1 = 0;
 712   AliTPCseed *s1 = 0;
 713   AliTPCseed *seed1 = 0;
 714   if (!treeInput) return;
 715   if (treeInput->GetEntries()==0) return;
 716   //
 717   treeInput->SetBranchAddress("v0.",&v0);
 718   treeInput->SetBranchAddress("t0.",&track0);
 719   treeInput->SetBranchAddress("ft0.",&friendTrack0);
 720   treeInput->SetBranchAddress("s0.",&s0);
 721   treeInput->SetBranchAddress("t1.",&track1);
 722   treeInput->SetBranchAddress("ft1.",&friendTrack1);
 723   treeInput->SetBranchAddress("s1.",&s1);
 724   //
 725   TDatabasePDG pdg;
 726   Int_t ncorr=0;
 727   if (corrArray) ncorr = corrArray->GetEntries();
 728   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
 729   Double_t massK0= pdg.GetParticle("K0")->Mass();
 730   Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
 731   Double_t massPion=pdg.GetParticle("pi+")->Mass();
 732   Double_t massProton=pdg.GetParticle("proton")->Mass();
 733   Int_t pdgPion=pdg.GetParticle("pi+")->PdgCode();
 734   Int_t pdgProton=pdg.GetParticle("proton")->PdgCode();
 735   Double_t rmass0=0;
 736   Double_t rmass1=0;
 737   Double_t massV0=0;
 738   Int_t    pdg0=0;
 739   Int_t    pdg1=0;
 740   //
 741   //
 742   //
 743   Int_t nv0s= treeInput->GetEntries();
 744   for (Int_t iv0=0; iv0<nv0s; iv0++){
 745     Int_t  v0Type=0;
 746     Int_t isK0=0;
 747     Int_t isLambda=0;
 748     Int_t isAntiLambda=0;
 749     treeInput->GetEntry(iv0);
 750     if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.03) {isK0=1; v0Type=1;} //select K0s
 751     if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.01) {isLambda=1; v0Type=2;} //select Lambda
 752     if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.01) {isAntiLambda=1;v0Type=3;} //select Anti Lambda
 753     if (isK0+isLambda+isAntiLambda!=1) continue;
 754     rmass0=massPion;
 755     rmass1=massPion;
 756     pdg0=pdgPion;
 757     pdg1=pdgPion;
 758     if (isLambda) {rmass0=massProton; pdg0=pdgProton;}
 759     if (isAntiLambda) {rmass1=massProton; pdg1=pdgProton;}
 760     massV0=massK0;
 761     if (isK0==0) massV0=massLambda;
 762     //
 763     if (!s0) continue;
 764     seed0=(s0->GetSign()>0)?s0:s1;
 765     seed1=(s0->GetSign()>0)?s1:s0;
 766     if (seed0->GetZ()*seed1->GetZ()<0) continue; //remove membrane crossed tracks
 767     if (seed0->Pt()<ptCut) continue;
 768     if (seed1->Pt()<ptCut) continue;
 769     //
 770     // Reapply transformation
 771     //
 772     for  (Int_t itype=0; itype<2; itype++){
 773       AliTPCseed * seed = (itype==0) ? seed0: seed1;
 774       for (Int_t irow=0; irow<159; irow++){
 775         AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
 776         if (cluster &&cluster->GetX()>10){
 777           Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
 778           Int_t index0[1]={cluster->GetDetector()};
 779           transform->Transform(x0,index0,0,1);
 780           cluster->SetX(x0[0]);
 781           cluster->SetY(x0[1]);
 782           cluster->SetZ(x0[2]);
 783           //
 784         }
 785       }
 786     }
 787     Bool_t isOK=kTRUE;
 788     Double_t radius = v0->GetRr();
 789     Double_t xyz[3];
 790     v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
 791     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
 792     TObjArray arrayV0in(ncorr+1);
 793     TObjArray arrayV0io(ncorr+1);
 794     TObjArray arrayT0(ncorr+1);
 795     TObjArray arrayT1(ncorr+1);
 796     arrayV0in.SetOwner(kTRUE);
 797     arrayV0io.SetOwner(kTRUE);
 798     //
 799     for (Int_t icorr=-1; icorr<ncorr; icorr++){
 800       AliTPCCorrection *corr =0;
 801       if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
 802       //
 803       AliExternalTrackParam * trackInner0 = RefitTrack(seed0, corr,160,85,rmass0);
 804       AliExternalTrackParam * trackIO0    = RefitTrack(seed0, corr,245,85,rmass0);
 805       AliExternalTrackParam * trackInner1 = RefitTrack(seed1, corr,160,85,rmass1);
 806       AliExternalTrackParam * trackIO1    = RefitTrack(seed1, corr,245,85,rmass1);
 807       if (!trackInner0) isOK=kFALSE;
 808       if (!trackInner1) isOK=kFALSE;
 809       if (!trackIO0)    isOK=kFALSE;
 810       if (!trackIO1)    isOK=kFALSE;
 811       if (isOK){
 812         if (!trackInner0->Rotate(alpha)) isOK=kFALSE;
 813         if (!trackInner1->Rotate(alpha)) isOK=kFALSE;
 814         if (!trackIO0->Rotate(alpha)) isOK=kFALSE;
 815         if (!trackIO1->Rotate(alpha)) isOK=kFALSE;
 816         //
 817         if (!AliTracker::PropagateTrackToBxByBz(trackInner0, radius, rmass0, 1, kFALSE)) isOK=kFALSE;
 818         if (!AliTracker::PropagateTrackToBxByBz(trackInner1, radius, rmass1, 1, kFALSE)) isOK=kFALSE;
 819         if (!AliTracker::PropagateTrackToBxByBz(trackIO0, radius, rmass0, 1, kFALSE)) isOK=kFALSE;
 820         if (!AliTracker::PropagateTrackToBxByBz(trackIO1, radius, rmass1, 1, kFALSE)) isOK=kFALSE;
 821         if (!isOK) continue;
 822         arrayT0.AddAt(trackIO0->Clone(),icorr+1);
 823         arrayT1.AddAt(trackIO1->Clone(),icorr+1);
 824         Int_t charge=TMath::Nint(trackIO0->GetSign());
 825         AliKFParticle pin0( *trackInner0,  pdg0*charge);
 826         AliKFParticle pin1( *trackInner1, -pdg1*charge);
 827         AliKFParticle pio0( *trackIO0,  pdg0*charge);
 828         AliKFParticle pio1( *trackIO1, -pdg1*charge);
 829         AliKFParticle v0in;
 830         AliKFParticle v0io;
 831         v0in+=pin0;
 832         v0in+=pin1;
 833         v0io+=pio0;
 834         v0io+=pio1;
 835         arrayV0in.AddAt(v0in.Clone(),icorr+1);
 836         arrayV0io.AddAt(v0io.Clone(),icorr+1);
 837       }
 838     }
 839     if (!isOK) continue;
 840     //
 841     //AliKFParticle* pin0= (AliKFParticle*)arrayV0in.At(0);
 842     AliKFParticle* pio0= (AliKFParticle*)arrayV0io.At(0);
 843     AliExternalTrackParam *param0=(AliExternalTrackParam *)arrayT0.At(0);
 844     AliExternalTrackParam *param1=(AliExternalTrackParam *)arrayT1.At(0);
 845     Double_t mass0=0, mass0E=0;
 846     pio0->GetMass( mass0,mass0E);
 847     //
 848     Double_t mean=mass0-massV0;
 849     if (TMath::Abs(mean)>0.05) continue;
 850     Double_t mass[10000];
 851     //
 852     Int_t dtype=30;  // id for V0
 853     Int_t ptype=5;   // id for invariant mass
 854     //    Int_t id=TMath::Nint(100.*(param0->Pt()-param1->Pt())/(param0->Pt()+param1->Pt()));      // K0s V0 asymetry
 855     Int_t id=Int_t(1000.*(param0->Pt()-param1->Pt()));      // K0s V0 asymetry
 856     Double_t gx,gy,gz, px,py,pz;
 857     Double_t pt = v0->Pt();
 858     v0->GetXYZ(gx,gy,gz);
 859     v0->GetPxPyPz(px,py,pz);
 860     Double_t theta=pz/TMath::Sqrt(px*px+py*py);
 861     Double_t phi=TMath::ATan2(py,px);
 862     Double_t snp=0.5*(seed0->GetSnp()+seed1->GetSnp());
 863     Double_t sector=9*phi/TMath::Pi();
 864     Double_t dsec=sector-TMath::Nint(sector);
 865     Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
 866     //Int_t nentries=v0Type;
 867     Double_t bz=AliTracker::GetBz();
 868     Double_t dRrec=0;
 869     (*pcstream)<<"fitDebug"<<
 870       "id="<<id<<
 871       "v0.="<<v0<<
 872       "mean="<<mean<<
 873       "rms="<<mass0E<<
 874       "pio0.="<<pio0<<
 875       "p0.="<<param0<<
 876       "p1.="<<param1;
 877     (*pcstream)<<"fit"<<  // dump valus for fit
 878       "run="<<run<<       //run number
 879       "bz="<<bz<<         // magnetic filed used
 880       "dtype="<<dtype<<   // detector match type 30
 881       "ptype="<<ptype<<   // parameter type
 882       "theta="<<theta<<   // theta
 883       "phi="<<phi<<       // phi
 884       "snp="<<snp<<       // snp
 885       "mean="<<mean<<     // mean dist value
 886       "rms="<<mass0E<<       // rms
 887       "sector="<<sector<<
 888       "dsec="<<dsec<<
 889       //
 890       "refX="<<refX<<      // reference radius
 891       "gx="<<gx<<         // global position
 892       "gy="<<gy<<         // global position
 893       "gz="<<gz<<         // global position
 894       "dRrec="<<dRrec<<      // delta Radius in reconstruction
 895       "pt="<<pt<<         // pt of the particle
 896       "id="<<id<<     //delta of the momenta
 897       "entries="<<v0Type;//  type of the V0
 898     for (Int_t icorr=0; icorr<ncorr; icorr++){
 899       AliTPCCorrection *corr =0;
 900       if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
 901       //      AliKFParticle* pin= (AliKFParticle*)arrayV0in.At(icorr+1);
 902       AliKFParticle* pio= (AliKFParticle*)arrayV0io.At(icorr+1);
 903       AliExternalTrackParam *par0=(AliExternalTrackParam *)arrayT0.At(icorr+1);
 904       AliExternalTrackParam *par1=(AliExternalTrackParam *)arrayT1.At(icorr+1);
 905       Double_t massE=0;
 906       pio->GetMass( mass[icorr],massE);
 907       mass[icorr]-=mass0;
 908       (*pcstream)<<"fit"<<
 909         Form("%s=",corr->GetName())<<mass[icorr];
 910       (*pcstream)<<"fitDebug"<<
 911         Form("%s=",corr->GetName())<<mass[icorr]<<
 912         Form("%sp0.=",corr->GetName())<<par0<<
 913         Form("%sp1=",corr->GetName())<<par1;
 914     }
 915     (*pcstream)<<"fit"<< "isOK="<<isOK<<"\n";
 916     (*pcstream)<<"fitDebug"<< "isOK="<<isOK<<"\n";
 917   }
 918   delete pcstream;
 919 }
 920
 921 AliExternalTrackParam * AliTPCcalibV0::RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass){
 922   //
 923   // Refit the track:
 924   //    seed   - tpc track with cluster
 925   //    corr   - distrotion/correction class  - apllied to the points
 926   //    xstart - radius to start propagate/update
 927   //    xstop  - radius to stop propagate/update
 928   //
 929   const Double_t kResetCov=20.;
 930   const Double_t kSigma=5.;
 931   Double_t covar[15];
 932   for (Int_t i=0;i<15;i++) covar[i]=0;
 933   covar[0]=kSigma*kSigma;
 934   covar[2]=kSigma*kSigma;
 935   covar[5]=kSigma*kSigma/Float_t(150*150);
 936   covar[9]=kSigma*kSigma/Float_t(150*150);
 937   covar[14]=1*1;
 938   //
 939   AliExternalTrackParam *refit  = new AliExternalTrackParam(*seed);
 940   refit->PropagateTo(xstart, AliTracker::GetBz());
 941   refit->AddCovariance(covar);
 942   refit->ResetCovariance(kResetCov);
 943   Double_t xmin = TMath::Min(xstart,xstop);
 944   Double_t xmax = TMath::Max(xstart,xstop);
 945   Int_t ncl=0;
 946   //
 947   Bool_t isOK=kTRUE;
 948   for (Int_t index0=0; index0<159; index0++){
 949     Int_t irow= (xstart<xstop)? index0:159-index0;
 950     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);  //cluster in local system
 951     if (!cluster) continue;
 952     if (cluster->GetX()<xmin) continue;
 953     if (cluster->GetX()>xmax) continue;
 954     Double_t alpha = TMath::Pi()*(cluster->GetDetector()+0.5)/9.;
 955     if (!refit->Rotate(alpha)) isOK=kFALSE;
 956     Double_t x     = cluster->GetX();
 957     Double_t y     = cluster->GetY();
 958     Double_t z     = cluster->GetZ();
 959     if (corr){
 960       Float_t xyz[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};  // original position
 961       corr->DistortPointLocal(xyz,cluster->GetDetector());
 962       x=xyz[0];
 963       y=xyz[1];
 964       z=xyz[2];
 965     }
 966     if (!AliTracker::PropagateTrackToBxByBz(refit, x,mass,1.,kFALSE)) isOK=kFALSE;
 967     if (!isOK) continue;
 968     Double_t cov[3]={0.01,0.,0.01};
 969     Double_t yz[2]={y,z};
 970     if (!refit->Update(yz,cov)) isOK=kFALSE;
 971     ncl++;
 972   }
 973   if (!AliTracker::PropagateTrackToBxByBz(refit, xstop, mass,1.,kTRUE)) isOK=kFALSE;
 974   //
 975   if (!isOK) {
 976     delete refit;
 977     return 0;
 978   }
 979   return refit;
 980 }
 981
 982
 983
 984
 985
 986 //
 987 // Obsolete part
 988 //
 989
 990
 991 void AliTPCcalibV0::MakeV0s(){
 992   //
 993   //
 994   //
 995   const Int_t kMinCluster=40;
 996   const Float_t kMinR    =0;
 997   if (! fV0s) fV0s = new TObjArray(10);
 998   fV0s->Clear();
 999   //
1000   // Old V0 finder
1001   //
1002   for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
1003     AliESDv0 * v0 = fESD->GetV0(ivertex);
1004     if (v0->GetOnFlyStatus()) continue;
1005     fV0s->AddLast(v0);
1006   }
1007   ProcessV0(0);
1008   fV0s->Clear(0);
1009   //
1010   // MI V0 finder
1011   //
1012   for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
1013     AliESDv0 * v0 = fESD->GetV0(ivertex);
1014     if (!v0->GetOnFlyStatus()) continue;
1015     fV0s->AddLast(v0);
1016   }
1017   ProcessV0(1);
1018   fV0s->Clear();
1019   return;
1020   //
1021   // combinatorial
1022   //
1023   Int_t ntracks = fESD->GetNumberOfTracks();
1024   for (Int_t itrack0=0; itrack0<ntracks; itrack0++){
1025     AliESDtrack * track0 = fESD->GetTrack(itrack0);
1026     if (track0->GetSign()>0) continue;
1027     if ( track0->GetTPCNcls()<kMinCluster) continue;
1028     if (track0->GetKinkIndex(0)>0) continue;
1029     //
1030     for (Int_t itrack1=0; itrack1<ntracks; itrack1++){
1031       AliESDtrack * track1 = fESD->GetTrack(itrack1);
1032       if (track1->GetSign()<0) continue;
1033       if ( track1->GetTPCNcls()<kMinCluster) continue;
1034       if (track1->GetKinkIndex(0)>0) continue;
1035       //
1036       //      AliExternalTrackParam param0(*track0);
1037       // AliExternalTrackParam param1(*track1);
1038       AliV0 vertex;
1039       vertex.SetParamN(*track0);
1040       vertex.SetParamP(*track1);
1041       Float_t xyz[3];
1042       xyz[0] = fESD->GetPrimaryVertex()->GetXv();
1043       xyz[1] = fESD->GetPrimaryVertex()->GetYv();
1044       xyz[2] = fESD->GetPrimaryVertex()->GetZv();
1045       vertex.Update(xyz);
1046       if (vertex.GetRr()<kMinR) continue;
1047       if (vertex.GetDcaV0Daughters()>1.) continue;
1048       if (vertex.GetDcaV0Daughters()>0.3*vertex.GetRr()) continue;
1049       // if (vertex.GetPointAngle()<0.9) continue;
1050       vertex.SetIndex(0,itrack0);
1051       vertex.SetIndex(1,itrack1);
1052       fV0s->AddLast(new AliV0(vertex));
1053     }
1054   }
1055   ProcessV0(2);
1056   for (Int_t i=0;i<fV0s->GetEntries(); i++) delete fV0s->At(i);
1057   fV0s->Clear();
1058 }
1059
1060
1061
1062
1063
1064
1065
1066
1067 void AliTPCcalibV0::ProcessV0(Int_t ftype){
1068   //
1069   // Obsolete
1070   //
1071   if (! fGammas) fGammas = new TObjArray(10);
1072   fGammas->Clear();
1073   Int_t nV0s  = fV0s->GetEntries();
1074   if (nV0s==0) return;
1075   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
1076   //
1077   for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
1078     AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex);
1079     AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); // negative track
1080     AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); // positive track
1081
1082     const AliExternalTrackParam * paramInNeg = trackN->GetInnerParam();
1083     const AliExternalTrackParam * paramInPos = trackP->GetInnerParam();
1084
1085     if (!paramInPos) continue; // in case the inner paramters do not exist
1086     if (!paramInNeg) continue;
1087     //
1088     //
1089     //
1090     AliKFParticle *v0K0       = Fit(primVtx,v0,-211,211);
1091     AliKFParticle *v0Gamma    = Fit(primVtx,v0,11,-11);
1092     AliKFParticle *v0Lambda42 = Fit(primVtx,v0,-2212,211);
1093     AliKFParticle *v0Lambda24 = Fit(primVtx,v0,-211,2212);
1094     //Set production vertex
1095     v0K0->SetProductionVertex( primVtx );
1096     v0Gamma->SetProductionVertex( primVtx );
1097     v0Lambda42->SetProductionVertex( primVtx );
1098     v0Lambda24->SetProductionVertex( primVtx );
1099     Double_t massK0, massGamma, massLambda42,massLambda24, massSigma;
1100     v0K0->GetMass( massK0,massSigma);
1101     v0Gamma->GetMass( massGamma,massSigma);
1102     v0Lambda42->GetMass( massLambda42,massSigma);
1103     v0Lambda24->GetMass( massLambda24,massSigma);
1104     Float_t chi2K0       = v0K0->GetChi2()/v0K0->GetNDF();
1105     Float_t chi2Gamma    = v0Gamma->GetChi2()/v0Gamma->GetNDF();
1106     Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF();
1107     Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF();
1108     //
1109     // Mass Contrained params
1110     //
1111     AliKFParticle *v0K0C       = Fit(primVtx,v0,-211,211);
1112     AliKFParticle *v0GammaC    = Fit(primVtx,v0,11,-11);
1113     AliKFParticle *v0Lambda42C = Fit(primVtx,v0,-2212,211); //lambdaBar
1114     AliKFParticle *v0Lambda24C = Fit(primVtx,v0,-211,2212); //lambda
1115     //
1116     v0K0C->SetProductionVertex( primVtx );
1117     v0GammaC->SetProductionVertex( primVtx );
1118     v0Lambda42C->SetProductionVertex( primVtx );
1119     v0Lambda24C->SetProductionVertex( primVtx );
1120
1121     v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass());
1122     v0GammaC->SetMassConstraint(0);
1123     v0Lambda42C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass());
1124     v0Lambda24C->SetMassConstraint(fPdg->GetParticle(3122)->Mass());
1125     //
1126     Double_t timeK0, sigmaTimeK0;
1127     Double_t timeLambda42, sigmaTimeLambda42;
1128     Double_t timeLambda24, sigmaTimeLambda24;
1129     v0K0C->GetLifeTime(timeK0, sigmaTimeK0);
1130     //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0);
1131     v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42);
1132     v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24);
1133
1134
1135     //
1136     Float_t chi2K0C       = v0K0C->GetChi2()/v0K0C->GetNDF();
1137     if (chi2K0C<0) chi2K0C=100;
1138     Float_t chi2GammaC    = v0GammaC->GetChi2()/v0GammaC->GetNDF();
1139     if (chi2GammaC<0) chi2GammaC=100;
1140     Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF();
1141     if (chi2Lambda42C<0) chi2Lambda42C=100;
1142     Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF();
1143     if (chi2Lambda24C<0) chi2Lambda24C=100;
1144     //
1145     Float_t  minChi2C=99;
1146     Int_t   type   =-1;
1147     if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;}
1148     if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;}
1149     if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;}
1150     if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;}
1151     Float_t  minChi2=99;
1152     Int_t   type0   =-1;
1153     if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;}
1154     if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;}
1155     if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;}
1156     if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;}
1157
1158      // 0 is  negative particle; 1 is positive particle
1159     Float_t betaGamma0 = 0;
1160     Float_t betaGamma1 = 0;
1161
1162     switch (type) {
1163      case 0:
1164       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass();
1165       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass();
1166       break;
1167      case 1:
1168       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(11)->Mass();
1169       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(-11)->Mass();
1170       break;
1171      case 2:
1172       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-2212)->Mass();
1173       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass();
1174       break;
1175      case 3:
1176       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass();
1177       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(2212)->Mass();
1178       break;
1179     }
1180
1181     // cuts and histogram filling
1182     Int_t numCand = 0; // number of particle types which have a chi2 < 10*minChi2
1183
1184     if (minChi2C < 2 && ftype == 1) {
1185      //
1186      if (chi2K0C < 10*minChi2C) numCand++;
1187      if (chi2GammaC < 10*minChi2C) numCand++;
1188      if (chi2Lambda42C < 10*minChi2C) numCand++;
1189      if (chi2Lambda24C < 10*minChi2C) numCand++;
1190      //
1191     }
1192
1193     //
1194     //
1195     // write output tree
1196     if (minChi2>50) continue;
1197     TTreeSRedirector *cstream = GetDebugStreamer();
1198     if (cstream){
1199       (*cstream)<<"V0"<<
1200         "ftype="<<ftype<<
1201         "v0.="<<v0<<
1202         "trackN.="<<trackN<<
1203         "trackP.="<<trackP<<
1204         //
1205         "betaGamma0="<<betaGamma0<<
1206         "betaGamma1="<<betaGamma1<<
1207         //
1208         "type="<<type<<
1209         "chi2C="<<minChi2C<<
1210         "v0K0.="<<v0K0<<
1211         "v0Gamma.="<<v0Gamma<<
1212         "v0Lambda42.="<<v0Lambda42<<
1213         "v0Lambda24.="<<v0Lambda24<<
1214         //
1215         "chi20K0.="<<chi2K0<<
1216         "chi2Gamma.="<<chi2Gamma<<
1217         "chi2Lambda42.="<<chi2Lambda42<<
1218         "chi2Lambda24.="<<chi2Lambda24<<
1219         //
1220         "chi20K0c.="<<chi2K0C<<
1221         "chi2Gammac.="<<chi2GammaC<<
1222         "chi2Lambda42c.="<<chi2Lambda42C<<
1223         "chi2Lambda24c.="<<chi2Lambda24C<<
1224         //
1225         "v0K0C.="<<v0K0C<<
1226         "v0GammaC.="<<v0GammaC<<
1227         "v0Lambda42C.="<<v0Lambda42C<<
1228         "v0Lambda24C.="<<v0Lambda24C<<
1229         //
1230         "massK0="<<massK0<<
1231         "massGamma="<<massGamma<<
1232         "massLambda42="<<massLambda42<<
1233         "massLambda24="<<massLambda24<<
1234         //
1235         "timeK0="<<timeK0<<
1236         "timeLambda42="<<timeLambda42<<
1237         "timeLambda24="<<timeLambda24<<
1238         "\n";
1239     }
1240     if (type==1) fGammas->AddLast(v0);
1241     //
1242     //
1243     //
1244     delete v0K0;
1245     delete v0Gamma;
1246     delete v0Lambda42;
1247     delete v0Lambda24;
1248     delete v0K0C;
1249     delete v0GammaC;
1250     delete v0Lambda42C;
1251     delete v0Lambda24C;
1252   }
1253   ProcessPI0();
1254 }
1255
1256
1257
1258 void AliTPCcalibV0::ProcessPI0(){
1259   //
1260   //
1261   //
1262   Int_t nentries = fGammas->GetEntries();
1263   if (nentries<2) return;
1264   //
1265   Double_t m0[3], m1[3];
1266   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
1267   for (Int_t i0=0; i0<nentries; i0++){
1268     AliESDv0 *v00 = (AliESDv0*)fGammas->At(i0);
1269     v00->GetPxPyPz (m0[0], m0[1], m0[2]);
1270     AliKFParticle *p00 = Fit(primVtx, v00, 11,-11);
1271     p00->SetProductionVertex( primVtx );
1272     p00->SetMassConstraint(0);
1273     //
1274     for (Int_t i1=i0; i1<nentries; i1++){
1275       AliESDv0 *v01 = (AliESDv0*)fGammas->At(i1);
1276       v01->GetPxPyPz (m1[0], m1[1], m1[2]);
1277       AliKFParticle *p01 = Fit(primVtx, v01, 11,-11);
1278       p01->SetProductionVertex( primVtx );
1279       p01->SetMassConstraint(0);
1280       if (v00->GetIndex(0) != v01->GetIndex(0) &&
1281           v00->GetIndex(1) != v01->GetIndex(1)){
1282         AliKFParticle pi0( *p00,*p01);
1283         pi0.SetProductionVertex(primVtx);
1284         Double_t n1 = TMath::Sqrt (m0[0]*m0[0] + m0[1]*m0[1] + m0[2]*m0[2]);
1285         Double_t n2 = TMath::Sqrt (m1[0]*m1[0] + m1[1]*m1[1] + m1[2]*m1[2]);
1286         Double_t mass = TMath::Sqrt(2.*(n1*n2 - (m0[0]*m1[0] + m0[1]*m1[1] + m0[2]*m1[2])));
1287         TTreeSRedirector *cstream = GetDebugStreamer();
1288         if (cstream){
1289           (*cstream)<<"PI0"<<
1290             "v00.="<<v00<<
1291             "v01.="<<v01<<
1292             "mass="<<mass<<
1293             "p00.="<<p00<<
1294             "p01.="<<p01<<
1295             "pi0.="<<&pi0<<
1296             "\n";
1297         }
1298       }
1299       delete p01;
1300     }
1301     delete p00;
1302   }
1303 }
1304
1305
1306
1307
1308
1309 AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){
1310   //
1311   // Make KF Particle
1312   //
1313   AliKFParticle p1( *(v0->GetParamN()), PDG1 );
1314   AliKFParticle p2( *(v0->GetParamP()), PDG2 );
1315   AliKFParticle *V0 = new AliKFParticle;
1316   Double_t x, y, z;
1317   v0->GetXYZ(x,y,z );
1318   V0->SetVtxGuess(x,y,z);
1319   *(V0)+=p1;
1320   *(V0)+=p2;
1321   return V0;
1322 }
1323
1324
1325
1326
1327 void AliTPCcalibV0::BinLogX(TH2F *h) {
1328   //
1329   //
1330   //
1331    TAxis *axis = h->GetXaxis();
1332    int bins = axis->GetNbins();
1333
1334    Double_t from = axis->GetXmin();
1335    Double_t to = axis->GetXmax();
1336    Double_t *new_bins = new Double_t[bins + 1];
1337
1338    new_bins[0] = from;
1339    Double_t factor = pow(to/from, 1./bins);
1340
1341    for (int i = 1; i <= bins; i++) {
1342      new_bins[i] = factor * new_bins[i-1];
1343    }
1344    axis->Set(bins, new_bins);
1345    delete [] new_bins;
1346
1347 }
1348
1349
1350
1351