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