PWG1/AliRecInfoMaker.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //  Time Projection Chamber                                                  //
  20 //  Comparison macro for ESD                                                 //
  21 //  responsible:
  22 //  marian.ivanov@cern.ch                                                    //
  23 //
  24 //
  25
  26 /*
  27 marian.ivanov@cern.ch
  28 Usage:
  29
  30
  31
  32 gSystem->Load("libPWG1.so");
  33 //
  34 AliRecInfoMaker *t2 = new AliRecInfoMaker("genTracks.root","cmpESDTracks.root","galice.root",-1,0,0);
  35 t2->Exec();
  36
  37
  38 TFile f("cmpESDTracks.root");
  39 AliTreeDraw comp;
  40 TTree * tree = (TTree*)f.Get("ESDcmpTracks");
  41 comp.SetTree(tree)
  42
  43
  44
  45 //
  46 //some cuts definition
  47 TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01")
  48 //TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.5&&abs(MC.fVDist[2])<0.5")
  49 //TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<3.9");
  50 TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5");
  51 TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
  52
  53
  54 TCut crec("crec","fReconstructed==1");
  55 TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25");
  56 TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1");
  57
  58 TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
  59 TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
  60 TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
  61 TCut cchi2("cchi2","fESDtrack.fITSchi2MIP[0]<7.&&fESDtrack.fITSchi2MIP[1]<5.&&fESDtrack.fITSchi2MIP[2]<7.&&fESDtrack.fITSchi2MIP[3]<7.5&&fESDtrack.fITSchi2MIP[4]<6.")
  62
  63
  64 //
  65 //example
  66 comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
  67 comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
  68 comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
  69 comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
  70 comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
  71 comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
  72 comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
  73
  74
  75 TH1F his("his","his",100,0,20);
  76 TH1F hpools("hpools","hpools",100,-7,7);
  77 TH1F hfake("hfake","hfake",1000,0,150);
  78 TProfile profp0("profp0","profp0",20,-0.4,0.9)
  79
  80 comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1"+cprim,"1",4,0.2,1.5,-0.06,0.06)
  81 comp.fRes->Draw();
  82 comp.fMean->Draw();
  83
  84 comp.DrawXY("fITSinP0[3]","fITSDelta[4]/fITSinP1[3]","fReconstructed==1&&fITSOn"+cprim,"1",4,0.2,1.5,-0.06,0.06)
  85 comp.fRes->Draw();
  86
  87 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
  88 comp.fRes->Draw();
  89
  90 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDtrack.fITSFakeRatio<0.1",10,0.2,1.5)
  91 comp.fRes->Draw();
  92 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDtrack.fITSFakeRatio>0.1",10,0.2,1.5)
  93 comp.fRes->Draw();
  94
  95 comp.fTree->Draw("fESDtrack.fITSsignal/fESDtrack.fTPCsignal","fITSOn&&fTPCOn&&fESDtrack.fITSFakeRatio==0")
  96
  97 TH1F his("his","his",100,0,20);
  98 TH1F hpools("hpools","hpools",100,-7,7);
  99
 100 TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,550); hdedx0->SetMarkerColor(1);
 101 TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,550); hdedx1->SetMarkerColor(4);
 102 TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,550); hdedx2->SetMarkerColor(3);
 103 TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,550); hdedx3->SetMarkerColor(2);
 104
 105 comp.fTree->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&abs(fPdg)==211&&fITStrack.fN==6"+cprim)
 106 comp.fTree->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim)
 107 comp.fTree->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&abs(fPdg)==321&&fITStrack.fN==6"+cprim)
 108 comp.fTree->Draw("fESDtrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&abs(fPdg)==11&&fITStrack.fN==6"+cprim)
 109
 110
 111 comp.fTree->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&fStatus[2]>1")
 112 comp.fTree->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&fStatus[2]>1")
 113 comp.fTree->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&fStatus[2]>1")
 114 comp.fTree->Draw("fESDtrack.fTRDsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&abs(fPdg)==11&&fTRDtrack.fN>40&&fStatus[2]>1")
 115
 116 comp.fTree->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx0","fTPCOn&&abs(fPdg)==211&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1);
 117 comp.fTree->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx1","fTPCOn&&abs(fPdg)==2212&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1);
 118 comp.fTree->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx2","fTPCOn&&abs(fPdg)==321&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1);
 119 comp.fTree->Draw("fESDtrack.fTPCsignal:fTPCinP0[4]>>dEdx3","fTPCOn&&abs(fPdg)==11&&fESDtrack.fTPCncls>180&&fESDtrack.fTPCsignal>10"+cteta1);
 120
 121 hdedx3->SetXTitle("P(GeV/c)");
 122 hdedx3->SetYTitle("dEdx(unit)");
 123 hdedx3->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx0->Draw("same");
 124
 125 comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
 126
 127 TProfile prof("prof","prof",10,0.5,5);
 128
 129
 130
 131
 132 */
 133
 134
 135 #include <stdio.h>
 136 #include <string.h>
 137 //ROOT includes
 138 #include "Rtypes.h"
 139 #include "TFile.h"
 140 #include "TTree.h"
 141 #include "TChain.h"
 142 #include "TCut.h"
 143 #include "TString.h"
 144 #include "TBenchmark.h"
 145 #include "TStopwatch.h"
 146 #include "TParticle.h"
 147 #include "TSystem.h"
 148 #include "TTimer.h"
 149 #include "TVector3.h"
 150 #include "TPad.h"
 151 #include "TCanvas.h"
 152 #include "TH1F.h"
 153 #include "TH2F.h"
 154 #include "TF1.h"
 155 #include "TText.h"
 156 #include "Getline.h"
 157 #include "TStyle.h"
 158 //ALIROOT includes
 159 #include "AliRun.h"
 160 #include "AliStack.h"
 161 #include "AliESDtrack.h"
 162 #include "AliSimDigits.h"
 163 #include "AliTPCParam.h"
 164 #include "AliTPC.h"
 165 #include "AliTPCLoader.h"
 166 #include "AliDetector.h"
 167 #include "AliTrackReference.h"
 168 #include "AliRun.h"
 169 #include "AliTPCParamSR.h"
 170 #include "AliTracker.h"
 171 #include "AliComplexCluster.h"
 172 #include "AliMagF.h"
 173 #include "AliESD.h"
 174 #include "AliESDfriend.h"
 175 #include "AliESDtrack.h"
 176 #include "AliTPCseed.h"
 177 #include "AliITStrackMI.h"
 178 #include "AliTRDtrack.h"
 179 #include "AliHelix.h"
 180 #include "AliESDVertex.h"
 181 #include "AliExternalTrackParam.h"
 182 #include "AliESDkink.h"
 183 #include "AliESDv0.h"
 184 #include "AliV0.h"
 185 //
 186 #include "AliTreeDraw.h"
 187 #include "AliGenInfo.h"
 188 #include "AliRecInfo.h"
 189 #include "AliRecInfoMaker.h"
 190
 191
 192
 193 ClassImp(AliRecInfoMaker)
 194
 195
 196
 197
 198 AliTPCParam * AliRecInfoMaker::GetTPCParam(){
 199   AliTPCParamSR * par = new AliTPCParamSR;
 200   par->Update();
 201   return par;
 202 }
 203
 204
 205
 206 void MakeAliases(AliTreeDraw&comp)
 207 {
 208   //
 209   // aliases definition
 210   //
 211   comp.T()->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
 212   comp.T()->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
 213   comp.T()->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
 214   comp.T()->SetAlias("theta","MC.fTrackRef.Theta()");
 215   comp.T()->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
 216   comp.T()->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
 217   comp.T()->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
 218
 219   comp.T()->SetAlias("trddedx","(RC.fESDtrack.fTRDsignals[0]+RC.fESDtrack.fTRDsignals[1]+RC.fESDtrack.fTRDsignals[2]+RC.fESDtrack.fTRDsignals[3]+RC.fESDtrack.fTRDsignals[4]+RC.fESDtrack.fTRDsignals[5])/6.");
 220
 221   comp.T()->SetAlias("dtofmc2","fESDtrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)");
 222   comp.T()->SetAlias("dtofrc2","(fESDtrack.fTrackTime[2]-fESDtrack.fTOFsignal)");
 223
 224   comp.T()->SetAlias("psum","fESDtrack.fTOFr[4]+fESDtrack.fTOFr[3]+fESDtrack.fTOFr[2]+fESDtrack.fTOFr[1]+fESDtrack.fTOFr[0]");
 225   comp.T()->SetAlias("P0","fESDtrack.fTOFr[0]/psum");
 226   comp.T()->SetAlias("P1","fESDtrack.fTOFr[1]/psum");
 227   comp.T()->SetAlias("P2","fESDtrack.fTOFr[2]/psum");
 228   comp.T()->SetAlias("P3","fESDtrack.fTOFr[3]/psum");
 229   comp.T()->SetAlias("P4","fESDtrack.fTOFr[4]/psum");
 230   comp.T()->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)");
 231 }
 232
 233
 234 ////////////////////////////////////////////////////////////////////////
 235 AliRecInfoMaker::AliRecInfoMaker()
 236 {
 237   Reset();
 238 }
 239
 240 ////////////////////////////////////////////////////////////////////////
 241 AliRecInfoMaker::AliRecInfoMaker(const char* fnGenTracks,
 242                    const char* fnCmp,
 243                    const char* fnGalice, Int_t direction,
 244                    Int_t nEvents, Int_t firstEvent)
 245 {
 246   Reset();
 247   //  fFnGenTracks = fnGenTracks;
 248   //  fFnCmp = fnCmp;
 249   sprintf(fFnGenTracks,"%s",fnGenTracks);
 250   sprintf(fFnCmp,"%s",fnCmp);
 251
 252   fFirstEventNr = firstEvent;
 253   fEventNr = firstEvent;
 254   fNEvents = nEvents;
 255   fDirection = direction;
 256   //
 257   fLoader = AliRunLoader::Open(fnGalice);
 258   if (gAlice){
 259     //delete gAlice->GetRunLoader();
 260     delete gAlice;
 261     gAlice = 0x0;
 262   }
 263   if (fLoader->LoadgAlice()){
 264     cerr<<"Error occured while l"<<endl;
 265   }
 266   Int_t nall = fLoader->GetNumberOfEvents();
 267   if (nEvents==0) {
 268     nEvents =nall;
 269     fNEvents=nall;
 270     fFirstEventNr=0;
 271   }
 272
 273   if (nall<=0){
 274     cerr<<"no events available"<<endl;
 275     fEventNr = 0;
 276     return;
 277   }
 278   if (firstEvent+nEvents>nall) {
 279     fEventNr = nall-firstEvent;
 280     cerr<<"restricted number of events availaible"<<endl;
 281   }
 282   AliMagF * magf = gAlice->Field();
 283   AliTracker::SetFieldMap(magf,0);
 284
 285 }
 286
 287
 288 ////////////////////////////////////////////////////////////////////////
 289 AliRecInfoMaker::~AliRecInfoMaker()
 290 {
 291   if (fLoader) {
 292     delete fLoader;
 293   }
 294 }
 295
 296 //////////////////////////////////////////////////////////////
 297 Int_t AliRecInfoMaker::SetIO()
 298 {
 299   //
 300   //
 301   CreateTreeCmp();
 302   if (!fTreeCmp) return 1;
 303   fParamTPC = GetTPCParam();
 304   //
 305   if (!ConnectGenTree()) {
 306     cerr<<"Cannot connect tree with generated tracks"<<endl;
 307     return 1;
 308   }
 309   return 0;
 310 }
 311
 312 //////////////////////////////////////////////////////////////
 313
 314 Int_t AliRecInfoMaker::SetIO(Int_t eventNr)
 315 {
 316   //
 317   //
 318   // SET INPUT
 319   //
 320   TFile f("AliESDs.root");
 321   //
 322
 323   TTree* tree = (TTree*) f.Get("esdTree");
 324   if (!tree) {
 325     Char_t ename[100];
 326     sprintf(ename,"%d",eventNr);
 327     fEvent = (AliESD*)f.Get(ename);
 328     if (!fEvent){
 329       sprintf(ename,"ESD%d",eventNr);
 330       fEvent = (AliESD*)f.Get(ename);
 331     }
 332   }
 333   else{
 334     tree->SetBranchStatus("*",1);
 335     tree->SetBranchAddress("ESD", &fEvent);
 336     tree->SetBranchAddress("ESDfriend.",&fESDfriend);
 337     tree->GetEntry(eventNr);
 338     fEvent->SetESDfriend(fESDfriend);
 339   }
 340
 341
 342   /*
 343   Char_t ename[100];
 344   sprintf(ename,"%d",eventNr);
 345   fEvent = (AliESD*)f.Get(ename);
 346   if (!fEvent){
 347     sprintf(ename,"ESD%d",eventNr);
 348     fEvent = (AliESD*)f.Get(ename);
 349   }
 350
 351   TTree* tree = (TTree*) f.Get("esdTree");
 352   if (!tree) {
 353     Error("CheckESD", "no ESD tree found");
 354     return kFALSE;
 355   }
 356   tree->SetBranchAddress("ESD", &fEvent);
 357   tree->GetEntry(eventNr);
 358   */
 359
 360   if (!fEvent) return 1;
 361
 362   return 0;
 363 }
 364
 365
 366
 367 ////////////////////////////////////////////////////////////////////////
 368 void AliRecInfoMaker::Reset()
 369 {
 370   fEventNr = 0;
 371   fNEvents = 0;
 372   fTreeCmp = 0;
 373   fTreeCmpKinks =0;
 374   fTreeCmpV0 =0;
 375   //  fFnCmp = "cmpTracks.root";
 376   fFileGenTracks = 0;
 377   fDebug = 0;
 378   //
 379   fParamTPC = 0;
 380   fEvent =0;
 381 }
 382
 383 ////////////////////////////////////////////////////////////////////////
 384 Int_t AliRecInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
 385 {
 386   fNEvents = nEvents;
 387   fFirstEventNr = firstEventNr;
 388   return Exec();
 389 }
 390
 391 ////////////////////////////////////////////////////////////////////////
 392 Int_t AliRecInfoMaker::Exec()
 393 {
 394   TStopwatch timer;
 395   timer.Start();
 396
 397   if (SetIO()==1)
 398     return 1;
 399
 400   fNextTreeGenEntryToRead = 0;
 401   fNextKinkToRead = 0;
 402   fNextV0ToRead   =0;
 403   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
 404   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
 405        eventNr++) {
 406     fEventNr = eventNr;
 407     SetIO(fEventNr);
 408     fNParticles = gAlice->GetEvent(fEventNr);
 409
 410     fIndexRecTracks = new Short_t[fNParticles*20];  //write at maximum 4 tracks corresponding to particle
 411     fIndexRecKinks  = new Short_t[fNParticles*20];  //write at maximum 20 tracks corresponding to particle
 412     fIndexRecV0  = new Short_t[fNParticles*20];  //write at maximum 20 tracks corresponding to particle
 413
 414     fFakeRecTracks = new Short_t[fNParticles];
 415     fMultiRecTracks = new Short_t[fNParticles];
 416     fMultiRecKinks = new Short_t[fNParticles];
 417     fMultiRecV0 = new Short_t[fNParticles];
 418
 419     for (Int_t i = 0; i<fNParticles; i++) {
 420       for (Int_t j=0;j<20;j++){
 421         fIndexRecTracks[20*i+j] = -1;
 422         fIndexRecKinks[20*i+j]  = -1;
 423         fIndexRecV0[20*i+j]  = -1;
 424       }
 425       fFakeRecTracks[i] = 0;
 426       fMultiRecTracks[i] = 0;
 427       fMultiRecKinks[i] = 0;
 428       fMultiRecV0[i] = 0;
 429     }
 430
 431     cout<<"Start to process event "<<fEventNr<<endl;
 432     cout<<"\tfNParticles = "<<fNParticles<<endl;
 433     if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
 434     if (TreeTLoop()>0) return 1;
 435
 436     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
 437     if (TreeGenLoop(eventNr)>0) return 1;
 438     BuildKinkInfo0(eventNr);
 439     //BuildV0Info(eventNr); // no V0 info for a moment
 440     fRecArray->Delete();
 441
 442     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
 443
 444     delete [] fIndexRecTracks;
 445     delete [] fIndexRecKinks;
 446     delete [] fIndexRecV0;
 447     delete [] fFakeRecTracks;
 448     delete [] fMultiRecTracks;
 449     delete [] fMultiRecKinks;
 450     delete [] fMultiRecV0;
 451   }
 452
 453   CloseOutputFile();
 454
 455   cerr<<"Exec finished"<<endl;
 456   timer.Stop();
 457   timer.Print();
 458   return 0;
 459
 460 }
 461 ////////////////////////////////////////////////////////////////////////
 462 Bool_t AliRecInfoMaker::ConnectGenTree()
 463 {
 464 //
 465 // connect all branches from the genTracksTree
 466 // use the same variables as for the new cmp tree, it may work
 467 //
 468   fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
 469   if (!fFileGenTracks) {
 470     cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
 471     return kFALSE;
 472   }
 473   fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
 474   if (!fTreeGenTracks) {
 475     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
 476         <<fFnGenTracks<<endl;
 477     return kFALSE;
 478   }
 479   //
 480   fMCInfo = new AliMCInfo;
 481   fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
 482   //
 483   //
 484   fTreeGenKinks = (TTree*)fFileGenTracks->Get("genKinksTree");
 485   if (!fTreeGenKinks) {
 486     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
 487         <<fFnGenTracks<<endl;
 488     //return kFALSE;
 489   }
 490   else{
 491     fGenKinkInfo = new AliGenKinkInfo;
 492     fTreeGenKinks->SetBranchAddress("MC",&fGenKinkInfo);
 493   }
 494
 495   fTreeGenV0 = (TTree*)fFileGenTracks->Get("genV0Tree");
 496   if (!fTreeGenV0) {
 497     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
 498         <<fFnGenTracks<<endl;
 499     //return kFALSE;
 500   }
 501   else{
 502     fGenV0Info = new AliGenV0Info;
 503     fTreeGenV0->SetBranchAddress("MC",&fGenV0Info);
 504   }
 505   //
 506   if (fDebug > 1) {
 507     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
 508   }
 509   return kTRUE;
 510 }
 511
 512
 513 ////////////////////////////////////////////////////////////////////////
 514 void AliRecInfoMaker::CreateTreeCmp()
 515 {
 516   fFileCmp = TFile::Open(fFnCmp,"RECREATE");
 517   if (!fFileCmp) {
 518     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
 519     return;
 520   }
 521   //
 522   //
 523   fTreeCmp    = new TTree("ESDcmpTracks","ESDcmpTracks");
 524   fMCInfo = new AliMCInfo;
 525   fRecInfo = new AliESDRecInfo;
 526   AliESDtrack * esdTrack = new AliESDtrack;
 527   //  AliITStrackMI * itsTrack = new AliITStrackMI;
 528   fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo,256000);
 529   fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo,256000);
 530   //  fTreeCmp->Branch("ITS","AliITStrackMI",&itsTrack);
 531   delete esdTrack;
 532   //
 533   //
 534   fTreeCmpKinks    = new TTree("ESDcmpKinks","ESDcmpKinks");
 535   fGenKinkInfo     = new AliGenKinkInfo;
 536   fRecKinkInfo     = new AliESDRecKinkInfo;
 537   fTreeCmpKinks->Branch("MC.","AliGenKinkInfo",&fGenKinkInfo,256000);
 538   fTreeCmpKinks->Branch("RC.","AliESDRecKinkInfo",&fRecKinkInfo,256000);
 539   //
 540   //
 541   fTreeCmpV0       = new TTree("ESDcmpV0","ESDcmpV0");
 542   fGenV0Info     = new AliGenV0Info;
 543   fRecV0Info     = new AliESDRecV0Info;
 544   fTreeCmpV0->Branch("MC.","AliGenV0Info",   &fGenV0Info,256000);
 545   fTreeCmpV0->Branch("RC.","AliESDRecV0Info",&fRecV0Info,256000);
 546   //
 547   fTreeCmp->AutoSave();
 548   fTreeCmpKinks->AutoSave();
 549   fTreeCmpV0->AutoSave();
 550 }
 551 ////////////////////////////////////////////////////////////////////////
 552 void AliRecInfoMaker::CloseOutputFile()
 553 {
 554   if (!fFileCmp) {
 555     cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
 556     return;
 557   }
 558   fFileCmp->cd();
 559   fTreeCmp->Write();
 560   delete fTreeCmp;
 561
 562   fFileCmp->Close();
 563   delete fFileCmp;
 564   return;
 565 }
 566 ////////////////////////////////////////////////////////////////////////
 567
 568 TVector3 AliRecInfoMaker::TR2Local(AliTrackReference *trackRef,
 569                             AliTPCParam *paramTPC) {
 570
 571   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
 572   Int_t index[4];
 573   paramTPC->Transform0to1(x,index);
 574   paramTPC->Transform1to2Ideal(x,index);
 575   return TVector3(x);
 576 }
 577 ////////////////////////////////////////////////////////////////////////
 578
 579 Int_t AliRecInfoMaker::TreeTLoop()
 580 {
 581   //
 582   // loop over all ESD reconstructed tracks and store info in memory
 583   //
 584   // + loop over all reconstructed kinks
 585   TStopwatch  timer;
 586   timer.Start();
 587   //
 588   Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks();
 589   Int_t nKinks = (Int_t) fEvent->GetNumberOfKinks();
 590   Int_t nV0MIs = (Int_t) fEvent->GetNumberOfV0s();
 591   fSignedKinks = new Short_t[nKinks];
 592   fSignedV0    = new Short_t[nV0MIs];
 593   //
 594   // load kinks to the memory
 595   for (Int_t i=0; i<nKinks;i++){
 596     AliESDkink * kink =fEvent->GetKink(i);
 597     fSignedKinks[i]=0;
 598   }
 599   //
 600   for (Int_t i=0; i<nV0MIs;i++){
 601     AliV0 * v0MI = (AliV0*)fEvent->GetV0(i);
 602     fSignedV0[i]=0;
 603   }
 604
 605   //
 606   //
 607   AliESDtrack * track=0;
 608   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
 609     //track = (AliESDtrack*)fTracks->UncheckedAt(iEntry);
 610     track = (AliESDtrack*)fEvent->GetTrack(iEntry);
 611     //
 612     Int_t label = track->GetLabel();
 613     Int_t absLabel = abs(label);
 614     if (absLabel < fNParticles) {
 615       //      fIndexRecTracks[absLabel] =  iEntry;
 616       if (label < 0) fFakeRecTracks[absLabel]++;
 617       if (fMultiRecTracks[absLabel]>0){
 618         if (fMultiRecTracks[absLabel]<20)
 619           fIndexRecTracks[absLabel*20+fMultiRecTracks[absLabel]] =  iEntry;
 620       }
 621       else
 622         fIndexRecTracks[absLabel*20] =  iEntry;
 623       fMultiRecTracks[absLabel]++;
 624     }
 625   }
 626   // sort reconstructed kinks
 627   //
 628   AliESDkink * kink=0;
 629   for (Int_t iEntry=0; iEntry<nKinks;iEntry++){
 630     kink = (AliESDkink*)fEvent->GetKink(iEntry);
 631     if (!kink) continue;
 632     //
 633     Int_t label0 = TMath::Abs(kink->GetLabel(0));
 634     Int_t label1 = TMath::Abs(kink->GetLabel(1));
 635     Int_t absLabel = TMath::Min(label0,label1);
 636     if (absLabel < fNParticles) {
 637       if (fMultiRecKinks[absLabel]>0){
 638         if (fMultiRecKinks[absLabel]<20)
 639           fIndexRecKinks[absLabel*20+fMultiRecKinks[absLabel]] =  iEntry;
 640       }
 641       else
 642         fIndexRecKinks[absLabel*20] =  iEntry;
 643       fMultiRecKinks[absLabel]++;
 644     }
 645   }
 646   // --sort reconstructed V0
 647   //
 648 //   AliV0 * v0MI=0;
 649 //   for (Int_t iEntry=0; iEntry<nV0MIs;iEntry++){
 650 //     v0MI = (AliV0*)fEvent->GetV0(iEntry);
 651 //     if (!v0MI) continue;
 652 //     //
 653 //     Int_t label0 = TMath::Abs(v0MI->GetLabel(0));
 654 //     Int_t label1 = TMath::Abs(v0MI->GetLabel(1));
 655 //     //
 656 //     for (Int_t i=0;i<2;i++){
 657 //       Int_t absLabel =  TMath::Abs(v0MI->GetLabel(i));
 658 //       if (absLabel < fNParticles) {
 659 //      if (fMultiRecV0[absLabel]>0){
 660 //        if (fMultiRecV0[absLabel]<20)
 661 //          fIndexRecV0[absLabel*20+fMultiRecV0[absLabel]] =  iEntry;
 662 //      }
 663 //      else
 664 //        fIndexRecV0[absLabel*20] =  iEntry;
 665 //      fMultiRecV0[absLabel]++;
 666 //       }
 667 //     }
 668 //   }
 669
 670
 671   printf("Time spended in TreeTLoop\n");
 672   timer.Print();
 673
 674   if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
 675   return 0;
 676 }
 677
 678 ////////////////////////////////////////////////////////////////////////
 679 Int_t AliRecInfoMaker::TreeGenLoop(Int_t eventNr)
 680 {
 681 //
 682 // loop over all entries for a given event, find corresponding
 683 // rec. track and store in the fTreeCmp
 684 //
 685   TStopwatch timer;
 686   timer.Start();
 687   Int_t entry = fNextTreeGenEntryToRead;
 688   Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
 689   cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
 690       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
 691   TBranch * branch = fTreeCmp->GetBranch("RC");
 692   TBranch * branchF = fTreeCmp->GetBranch("F");
 693
 694   branch->SetAddress(&fRecInfo); // set all pointers
 695   fRecArray = new TObjArray(fNParticles);
 696   AliESDtrack dummytrack;  //
 697   AliESDfriendTrack dummytrackF;  //
 698
 699   while (entry < nParticlesTR) {
 700     fTreeGenTracks->GetEntry(entry);
 701     entry++;
 702     if (eventNr < fMCInfo->fEventNr) continue;
 703     if (eventNr > fMCInfo->fEventNr) continue;;
 704     //
 705     fNextTreeGenEntryToRead = entry-1;
 706     if (fDebug > 2 && fMCInfo->fLabel < 10) {
 707       cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
 708     }
 709     //    if (fMCInfo->fNTPCRef<1) continue; // not TPCref
 710     //
 711     fRecInfo->Reset();
 712     AliESDtrack * track=0;
 713     fRecInfo->fReconstructed =0;
 714     TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
 715     local.GetXYZ(fRecInfo->fTRLocalCoord);
 716     //
 717     if (fIndexRecTracks[fMCInfo->fLabel*20] >= 0) {
 718       //track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
 719       track= (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20]);
 720       //
 721       //
 722       // find nearest track if multifound
 723       //Int_t sign = Int_t(track->GetSign()*fMCInfo->fCharge);
 724       //
 725       Int_t status = 0;
 726       if  ((track->GetStatus()&AliESDtrack::kITSrefit)>0) status++;
 727       if  ((track->GetStatus()&AliESDtrack::kTPCrefit)>0) status++;
 728       if  ((track->GetStatus()&AliESDtrack::kTRDrefit)>0) status++;
 729
 730       //
 731       if (fIndexRecTracks[fMCInfo->fLabel*20+1]>0){
 732         //
 733         Double_t p[3];
 734         track->GetInnerPxPyPz(p);
 735         Float_t maxp = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
 736         //
 737         for (Int_t i=1;i<20;i++){
 738           if (fIndexRecTracks[fMCInfo->fLabel*20+i]>=0){
 739             AliESDtrack * track2 = (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20+i]);
 740             if (!track2) continue;
 741             //Int_t sign2 = track->GetSign()*fMCInfo->fCharge; //
 742             //if (sign2<0) continue;
 743             track2->GetInnerPxPyPz(p);
 744             Float_t mom = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
 745             /*
 746             if (sign<0){
 747               sign = sign2;
 748               track = track2;
 749               maxp = mom;
 750               continue;
 751             }
 752             */
 753             //
 754             Int_t status2 = 0;
 755             if  ((track2->GetStatus()&AliESDtrack::kITSrefit)>0) status2++;
 756             if  ((track2->GetStatus()&AliESDtrack::kTPCrefit)>0) status2++;
 757             if  ((track2->GetStatus()&AliESDtrack::kTRDrefit)>0) status2++;
 758             if (status2<status) continue;
 759             //
 760             if (mom<maxp) continue;
 761             maxp = mom;
 762             track = track2;
 763             //
 764           }
 765         }
 766       }
 767       //
 768       if (track) {
 769         fRecInfo->SetESDtrack(track);
 770       }else{
 771         fRecInfo->SetESDtrack(&dummytrack);
 772       }
 773       //
 774
 775       fRecInfo->fReconstructed = 1;
 776       fRecInfo->fFake     = fFakeRecTracks[fMCInfo->fLabel];
 777       fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
 778       //
 779       fRecInfo->Update(fMCInfo,fParamTPC,kTRUE, fEvent);
 780     }
 781     else{
 782       fRecInfo->SetESDtrack(&dummytrack);
 783       fRecInfo->Update(fMCInfo,fParamTPC,kFALSE, fEvent);
 784     }
 785     fRecArray->AddAt(new AliESDRecInfo(*fRecInfo),fMCInfo->fLabel);
 786     fTreeCmp->Fill();
 787   }
 788   fTreeCmp->AutoSave();
 789   //fTracks->Delete();
 790   printf("Time spended in TreeGenLoop\n");
 791   timer.Print();
 792   if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
 793
 794   return 0;
 795 }
 796
 797
 798
 799 ////////////////////////////////////////////////////////////////////////
 800 ////////////////////////////////////////////////////////////////////////
 801 ////////////////////////////////////////////////////////////////////////
 802 Int_t AliRecInfoMaker::BuildKinkInfo0(Int_t eventNr)
 803 {
 804 //
 805 // loop over all entries for a given event, find corresponding
 806 // rec. track and store in the fTreeCmp
 807 //
 808   TStopwatch timer;
 809   timer.Start();
 810   Int_t entry = fNextKinkToRead;
 811   Double_t nParticlesTR = fTreeGenKinks->GetEntriesFast();
 812   cerr<<"fNParticles, nParticlesTR, fNextKinkToRead: "<<fNParticles<<" "
 813       <<nParticlesTR<<" "<<fNextKinkToRead<<endl;
 814   //
 815   TBranch * branch = fTreeCmpKinks->GetBranch("RC.");
 816   branch->SetAddress(&fRecKinkInfo); // set all pointers
 817
 818   //
 819   while (entry < nParticlesTR) {
 820     fTreeGenKinks->GetEntry(entry);
 821     entry++;
 822     if (eventNr < fGenKinkInfo->GetMinus().fEventNr) continue;
 823     if (eventNr > fGenKinkInfo->GetMinus().fEventNr) continue;;
 824     //
 825     fNextKinkToRead = entry-1;
 826     //
 827     //
 828     AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->GetMinus().fLabel);
 829     AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->GetPlus().fLabel);
 830     fRecKinkInfo->fT1 = (*fRecInfo1);
 831     fRecKinkInfo->fT2 = (*fRecInfo2);
 832     fRecKinkInfo->fStatus =0;
 833     if (fRecInfo1 && fRecInfo1->fTPCOn) fRecKinkInfo->fStatus+=1;
 834     if (fRecInfo2 && fRecInfo2->fTPCOn) fRecKinkInfo->fStatus+=2;
 835     if (fRecKinkInfo->fStatus==3&&fRecInfo1->fSign!=fRecInfo2->fSign) fRecKinkInfo->fStatus*=-1;
 836
 837     if (fRecKinkInfo->fStatus==3){
 838       fRecKinkInfo->Update();
 839     }
 840     Int_t label =  TMath::Min(fGenKinkInfo->GetMinus().fLabel,fGenKinkInfo->GetPlus().fLabel);
 841     Int_t label2 = TMath::Max(fGenKinkInfo->GetMinus().fLabel,fGenKinkInfo->GetPlus().fLabel);
 842
 843     AliESDkink *kink=0;
 844     fRecKinkInfo->fRecStatus   =0;
 845     fRecKinkInfo->fMultiple    = fMultiRecKinks[label];
 846     fRecKinkInfo->fKinkMultiple=0;
 847     //
 848     if (fMultiRecKinks[label]>0){
 849
 850       //      for (Int_t j=0;j<TMath::Min(fMultiRecKinks[label],100);j++){
 851       for (Int_t j=TMath::Min(fMultiRecKinks[label],Short_t(20))-1;j>=0;j--){
 852         Int_t index = fIndexRecKinks[label*20+j];
 853         //AliESDkink *kink2  = (AliESDkink*)fKinks->At(index);
 854         AliESDkink *kink2  = (AliESDkink*)fEvent->GetKink(index);
 855         if (TMath::Abs(kink2->GetLabel(0))==label &&TMath::Abs(kink2->GetLabel(1))==label2) {
 856           fRecKinkInfo->fKinkMultiple++;
 857           fSignedKinks[index]=1;
 858           Int_t c0=0;
 859           if (kink){
 860             //      if (kink->fTRDOn) c0++;
 861             //if (kink->fITSOn) c0++;
 862             if (kink->GetStatus(2)>0) c0++;
 863             if (kink->GetStatus(0)>0) c0++;
 864           }
 865           Int_t c2=0;
 866           //      if (kink2->fTRDOn) c2++;
 867           //if (kink2->fITSOn) c2++;
 868           if (kink2->GetStatus(2)>0) c2++;
 869           if (kink2->GetStatus(0)>0) c2++;
 870
 871           if (c2<c0) continue;
 872           kink =kink2;
 873         }
 874         if (TMath::Abs(kink2->GetLabel(1))==label &&TMath::Abs(kink2->GetLabel(0))==label2) {
 875           fRecKinkInfo->fKinkMultiple++;
 876           fSignedKinks[index]=1;
 877           Int_t c0=0;
 878           if (kink){
 879             //if (kink->fTRDOn) c0++;
 880             //if (kink->fITSOn) c0++;
 881             if (kink->GetStatus(2)>0) c0++;
 882             if (kink->GetStatus(0)>0) c0++;
 883
 884           }
 885           Int_t c2=0;
 886           //      if (kink2->fTRDOn) c2++;
 887           //if (kink2->fITSOn) c2++;
 888           if (kink2->GetStatus(2)>0) c2++;
 889           if (kink2->GetStatus(0)>0) c2++;
 890
 891           if (c2<c0) continue;
 892           kink =kink2;
 893         }
 894       }
 895     }
 896     if (kink){
 897       fRecKinkInfo->fKink = *kink;
 898       fRecKinkInfo->fRecStatus=1;
 899     }
 900     fTreeCmpKinks->Fill();
 901   }
 902   //  Int_t nkinks = fKinks->GetEntriesFast();
 903   Int_t nkinks = fEvent->GetNumberOfKinks();
 904   for (Int_t i=0;i<nkinks;i++){
 905     if (fSignedKinks[i]==0){
 906       //      AliESDkink *kink  = (AliESDkink*)fKinks->At(i);
 907       AliESDkink *kink  = (AliESDkink*)fEvent->GetKink(i);
 908       if (!kink) continue;
 909       //
 910       fRecKinkInfo->fKink = *kink;
 911       fRecKinkInfo->fRecStatus =-2;
 912       //
 913       AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->GetLabel(0)));
 914       AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->GetLabel(1)));
 915       if (fRecInfo1 && fRecInfo2){
 916         fRecKinkInfo->fT1 = (*fRecInfo1);
 917         fRecKinkInfo->fT2 = (*fRecInfo2);
 918         fRecKinkInfo->fRecStatus =-1;
 919       }
 920       fTreeCmpKinks->Fill();
 921     }
 922   }
 923
 924
 925   fTreeCmpKinks->AutoSave();
 926   printf("Time spended in BuilKinkInfo Loop\n");
 927   timer.Print();
 928   if (fDebug > 2) cerr<<"end of BuildKinkInfo Loop"<<endl;
 929   return 0;
 930 }
 931
 932
 933
 934 ////////////////////////////////////////////////////////////////////////
 935 ////////////////////////////////////////////////////////////////////////
 936 ////////////////////////////////////////////////////////////////////////
 937
 938
 939
 940 Int_t AliRecInfoMaker::BuildV0Info(Int_t eventNr)
 941 {
 942 //
 943 // loop over all entries for a given event, find corresponding
 944 // rec. track and store in the fTreeCmp
 945 //
 946   TStopwatch timer;
 947   timer.Start();
 948   Int_t entry = fNextV0ToRead;
 949   Double_t nParticlesTR = fTreeGenV0->GetEntriesFast();
 950   cerr<<"fNParticles, nParticlesTR, fNextV0ToRead: "<<fNParticles<<" "
 951       <<nParticlesTR<<" "<<fNextV0ToRead<<endl;
 952   //
 953   TBranch * branch = fTreeCmpV0->GetBranch("RC.");
 954   branch->SetAddress(&fRecV0Info); // set all pointers
 955   const AliESDVertex * esdvertex = fEvent->GetVertex();
 956   Float_t vertex[3]= {esdvertex->GetXv(), esdvertex->GetYv(),esdvertex->GetZv()};
 957
 958   //
 959   while (entry < nParticlesTR) {
 960     fTreeGenV0->GetEntry(entry);
 961     entry++;
 962     if (eventNr < fGenV0Info->GetMinus().fEventNr) continue;
 963     if (eventNr > fGenV0Info->GetMinus().fEventNr) continue;;
 964     //
 965     fNextV0ToRead = entry-1;
 966     //
 967     //
 968     AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->GetMinus().fLabel);
 969     AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->GetPlus().fLabel);
 970     if (fGenV0Info->GetMinus().fCharge*fGenV0Info->GetPlus().fCharge>0) continue;  // interactions
 971     if (!fRecInfo1 || !fRecInfo2) continue;
 972     fRecV0Info->fT1 = (*fRecInfo1);
 973     fRecV0Info->fT2 = (*fRecInfo2);
 974     fRecV0Info->fV0Status =0;
 975     if (fRecInfo1 && fRecInfo1->fStatus[1]>0) fRecV0Info->fV0Status+=1;
 976     if (fRecInfo2 && fRecInfo2->fStatus[1]>0) fRecV0Info->fV0Status+=2;
 977
 978     if (fRecV0Info->fV0Status==3&&fRecInfo1->fSign==fRecInfo2->fSign) fRecV0Info->fV0Status*=-1;
 979
 980
 981     if (abs(fRecV0Info->fV0Status)==3){
 982       fRecV0Info->Update(vertex);
 983       {
 984         //
 985         // TPC V0 Info
 986         Double_t x,alpha, param[5],cov[15];
 987         if ( fRecV0Info->fT1.GetESDtrack()->GetInnerParam() && fRecV0Info->fT2.GetESDtrack()->GetInnerParam()){
 988           fRecV0Info->fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
 989           fRecV0Info->fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
 990           AliExternalTrackParam paramP(x,alpha,param,cov);
 991           //
 992           fRecV0Info->fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
 993           fRecV0Info->fT2.GetESDtrack()->GetInnerExternalCovariance(cov);
 994           AliExternalTrackParam paramM(x,alpha,param,cov);
 995           //
 996           fRecV0Info->fV0tpc->SetParamN(paramM);
 997           fRecV0Info->fV0tpc->SetParamP(paramP);
 998           Double_t pid1[5],pid2[5];
 999           fRecV0Info->fT1.GetESDtrack()->GetESDpid(pid1);
1000           fRecV0Info->fT1.GetESDtrack()->GetESDpid(pid2);
1001           //
1002           //fRecV0Info->fV0tpc.UpdatePID(pid1,pid2);
1003           fRecV0Info->fV0tpc->Update(vertex);
1004
1005           //
1006           //
1007           fRecV0Info->fT1.GetESDtrack()->GetExternalParameters(x,param);
1008           fRecV0Info->fT1.GetESDtrack()->GetExternalCovariance(cov);
1009           alpha = fRecV0Info->fT1.GetESDtrack()->GetAlpha();
1010           new (&paramP) AliExternalTrackParam(x,alpha,param,cov);
1011           //
1012           fRecV0Info->fT2.GetESDtrack()->GetExternalParameters(x,param);
1013           fRecV0Info->fT2.GetESDtrack()->GetExternalCovariance(cov);
1014           alpha = fRecV0Info->fT2.GetESDtrack()->GetAlpha();
1015           new (&paramM) AliExternalTrackParam(x,alpha,param,cov);
1016           //
1017           fRecV0Info->fV0its->SetParamN(paramM);
1018           fRecV0Info->fV0its->SetParamP(paramP);
1019           //      fRecV0Info->fV0its.UpdatePID(pid1,pid2);
1020           fRecV0Info->fV0its->Update(vertex);
1021         }
1022       }
1023       if (TMath::Abs(fGenV0Info->GetMinus().fPdg)==11 &&TMath::Abs(fGenV0Info->GetPlus().fPdg)==11){
1024         if (fRecV0Info->fDist2>10){
1025           fRecV0Info->Update(vertex);
1026         }
1027         if (fRecV0Info->fDist2>10){
1028           fRecV0Info->Update(vertex);
1029         }
1030       }
1031     }
1032     //
1033     // take the V0 from reconstruction
1034
1035     Int_t label =  TMath::Min(fGenV0Info->GetMinus().fLabel,fGenV0Info->GetPlus().fLabel);
1036     Int_t label2 = TMath::Max(fGenV0Info->GetMinus().fLabel,fGenV0Info->GetPlus().fLabel);
1037     AliV0 *v0MI=0;
1038     fRecV0Info->fRecStatus   =0;
1039     fRecV0Info->fMultiple    = fMultiRecV0[label];
1040     fRecV0Info->fV0Multiple=0;
1041     //
1042     if (fMultiRecV0[label]>0 || fMultiRecV0[label2]>0){
1043
1044       //      for (Int_t j=0;j<TMath::Min(fMultiRecV0s[label],100);j++){
1045      //  for (Int_t j=TMath::Min(fMultiRecV0[label],Short_t(20))-1;j>=0;j--){
1046 //      Int_t index = fIndexRecV0[label*20+j];
1047 //      if (index<0) continue;
1048 //      AliV0 *v0MI2  = (AliV0*)fEvent->GetV0(index);
1049 //      if (TMath::Abs(v0MI2->GetLabel(0))==label &&TMath::Abs(v0MI2->GetLabel(1))==label2) {
1050 //        v0MI =v0MI2;
1051 //        fRecV0Info->fV0Multiple++;
1052 //        fSignedV0[index]=1;
1053 //      }
1054 //      if (TMath::Abs(v0MI2->GetLabel(1))==label &&TMath::Abs(v0MI2->GetLabel(0))==label2) {
1055 //        v0MI =v0MI2;
1056 //        fRecV0Info->fV0Multiple++;
1057 //        fSignedV0[index]=1;
1058 //      }
1059 //       }
1060     }
1061     if (v0MI){
1062       fRecV0Info->fV0rec = v0MI;
1063       fRecV0Info->fRecStatus=1;
1064     }
1065
1066     fTreeCmpV0->Fill();
1067   }
1068   //
1069   // write fake v0s
1070
1071   Int_t nV0MIs = fEvent->GetNumberOfV0s();
1072   for (Int_t i=0;i<nV0MIs;i++){
1073     if (fSignedV0[i]==0){
1074       AliV0 *v0MI  = (AliV0*)fEvent->GetV0(i);
1075       if (!v0MI) continue;
1076       //
1077       fRecV0Info->fV0rec = v0MI;
1078       fRecV0Info->fV0Status  =-10;
1079       fRecV0Info->fRecStatus =-2;
1080       //
1081  //      AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->GetLabel(0)));
1082 //       AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->GetLabel(1)));
1083 //       if (fRecInfo1 && fRecInfo2){
1084 //      fRecV0Info->fT1 = (*fRecInfo1);
1085 //      fRecV0Info->fT2 = (*fRecInfo2);
1086 //      fRecV0Info->fRecStatus =-1;
1087 //       }
1088       fRecV0Info->Update(vertex);
1089       fTreeCmpV0->Fill();
1090     }
1091   }
1092
1093
1094
1095   fTreeCmpV0->AutoSave();
1096   printf("Time spended in BuilV0Info Loop\n");
1097   timer.Print();
1098   if (fDebug > 2) cerr<<"end of BuildV0Info Loop"<<endl;
1099   return 0;
1100 }
1101 ////////////////////////////////////////////////////////////////////////
1102 ////////////////////////////////////////////////////////////////////////
1103
1104