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