PWGPP/TPC/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("libPWGPP.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   fEventNr(0),                 //! current event number
 232   fNEvents(0),                 //! number of events to process
 233   fFirstEventNr(0),            //! first event to process
 234   fFileCmp(0),                //! output file with cmp tracks
 235   fTreeCmp(0),                //! output tree with cmp tracks
 236   fTreeCmpKinks(0),                //! output tree with cmp Kinks
 237   fTreeCmpV0(0),                //! output tree with cmp V0
 238   //
 239   fFileGenTracks(0),                //! input files with generated tracks
 240   fTreeGenTracks(0),           //! tree with generated tracks
 241   fTreeGenKinks(0),            // tree with gen kinks
 242   fTreeGenV0(0),            // tree with gen V0
 243   //
 244   fLoader(0),         //! pointer to the run loader
 245   //
 246   fIndexRecTracks(0),         //! index of particle label in the TreeT_ESD
 247   fFakeRecTracks(0),          //! number of fake tracks
 248   fMultiRecTracks(0),         //! number of multiple reconstructions
 249   //
 250   fIndexRecKinks(0),         //! index of particle label in treeesd
 251   fMultiRecKinks(0),         //! number of multiple reconstructions
 252   fSignedKinks(0),           //! indicator that kink was not fake
 253   //
 254   fIndexRecV0(0),         //! index of particle label in treeesd
 255   fMultiRecV0(0),         //! number of multiple reconstructions
 256   fSignedV0(0),                //! indicator that kink was not fake
 257   //
 258   fRecArray(0),           // container with rec infos
 259   fEvent(0),             //!event
 260   fESDfriend(0),              //!event friend
 261   //
 262   fParamTPC(0),         //! AliTPCParam
 263   fNParticles(0),              //! number of particles in the input tree genTracks
 264   fDebug(0),                   //! debug flag
 265   fNextTreeGenEntryToRead(0),    //! last entry already read from genTracks tree
 266   fNextKinkToRead(0),            //! last entry already read from genKinks tree
 267   fNextV0ToRead(0),            //! last entry already read from genV0 tree
 268   //
 269   fMCInfo(0),           //! MC information writen per particle
 270   fGenKinkInfo(0),      //! MC information writen per Kink
 271   fGenV0Info(0),      //! MC information writen per Kink
 272   fRecInfo(0),          //! Rec. information writen per particle
 273   fFriend(0),          //! friend track
 274   fRecKinkInfo(0),    //! reconstructed kink info
 275   fRecV0Info(0)    //! reconstructed kink info
 276 {
 277   // AliRecInfoMaker - connencts the MC information with reconstructed information
 278   // fnGenTracks  - file with MC to be created before using AliGenInfoMaker
 279   // fnCmp        - file name  to be created
 280   // fnGalice     - file with Loaders - usualy galice.root
 281   //
 282   // nEvent       - number of event s to be analyzed
 283   // AliRecInfoMaker *t2 = new AliRecInfoMaker("genTracks.root","cmpESDTracks.root","galice.root",0,0);
 284   //
 285
 286
 287   Reset();
 288   //  fFnGenTracks = fnGenTracks;
 289   //  fFnCmp = fnCmp;
 290
 291   memset(fFnGenTracks,0,sizeof(fFnGenTracks));
 292   memset(fFnCmp,0,sizeof(fFnCmp));
 293
 294   snprintf(fFnGenTracks,1000,"%s",fnGenTracks);
 295   snprintf(fFnCmp,1000,"%s",fnCmp);
 296
 297   fFirstEventNr = firstEvent;
 298   fEventNr = firstEvent;
 299   fNEvents = nEvents;
 300   //
 301   fLoader = AliRunLoader::Open(fnGalice);
 302   if (gAlice){
 303     //delete AliRunLoader::Instance();
 304     delete gAlice;
 305     gAlice = 0x0;
 306   }
 307   if (fLoader->LoadgAlice()){
 308     cerr<<"Error occured while l"<<endl;
 309   }
 310   Int_t nall = fLoader->GetNumberOfEvents();
 311   if (nEvents==0) {
 312     nEvents =nall;
 313     fNEvents=nall;
 314     fFirstEventNr=0;
 315   }
 316
 317   if (nall<=0){
 318     cerr<<"no events available"<<endl;
 319     fEventNr = 0;
 320     return;
 321   }
 322   if (firstEvent+nEvents>nall) {
 323     fEventNr = nall-firstEvent;
 324     cerr<<"restricted number of events availaible"<<endl;
 325   }
 326   TGeoManager::Import("geometry.root");
 327
 328
 329 }
 330 ////////////////////////////////////////////////////////////////////////
 331 AliRecInfoMaker::AliRecInfoMaker(const AliRecInfoMaker& /*info*/):
 332
 333   fEventNr(0),                 //! current event number
 334   fNEvents(0),                 //! number of events to process
 335   fFirstEventNr(0),            //! first event to process
 336   fFileCmp(0),                //! output file with cmp tracks
 337   fTreeCmp(0),                //! output tree with cmp tracks
 338   fTreeCmpKinks(0),                //! output tree with cmp Kinks
 339   fTreeCmpV0(0),                //! output tree with cmp V0
 340   //
 341   fFileGenTracks(0),                //! input files with generated tracks
 342   fTreeGenTracks(0),           //! tree with generated tracks
 343   fTreeGenKinks(0),            // tree with gen kinks
 344   fTreeGenV0(0),            // tree with gen V0
 345   //
 346   fLoader(0),         //! pointer to the run loader
 347   //
 348   fIndexRecTracks(0),         //! index of particle label in the TreeT_ESD
 349   fFakeRecTracks(0),          //! number of fake tracks
 350   fMultiRecTracks(0),         //! number of multiple reconstructions
 351   //
 352   fIndexRecKinks(0),         //! index of particle label in treeesd
 353   fMultiRecKinks(0),         //! number of multiple reconstructions
 354   fSignedKinks(0),           //! indicator that kink was not fake
 355   //
 356   fIndexRecV0(0),         //! index of particle label in treeesd
 357   fMultiRecV0(0),         //! number of multiple reconstructions
 358   fSignedV0(0),                //! indicator that kink was not fake
 359   //
 360   fRecArray(0),           // container with rec infos
 361   fEvent(0),             //!event
 362   fESDfriend(0),              //!event friend
 363   //
 364   fParamTPC(0),         //! AliTPCParam
 365   fNParticles(0),              //! number of particles in the input tree genTracks
 366   fDebug(0),                   //! debug flag
 367   fNextTreeGenEntryToRead(0),    //! last entry already read from genTracks tree
 368   fNextKinkToRead(0),            //! last entry already read from genKinks tree
 369   fNextV0ToRead(0),            //! last entry already read from genV0 tree
 370   //
 371   fMCInfo(0),           //! MC information writen per particle
 372   fGenKinkInfo(0),      //! MC information writen per Kink
 373   fGenV0Info(0),      //! MC information writen per Kink
 374   fRecInfo(0),          //! Rec. information writen per particle
 375   fFriend(0),          //! friend track
 376   fRecKinkInfo(0),    //! reconstructed kink info
 377   fRecV0Info(0)    //! reconstructed kink info
 378 {
 379   //
 380   // Dummy copu constructor
 381   //
 382   memset(fFnGenTracks,0,sizeof(fFnGenTracks));
 383   memset(fFnCmp,0,sizeof(fFnCmp));
 384 }
 385
 386
 387
 388
 389 ////////////////////////////////////////////////////////////////////////
 390 AliRecInfoMaker::~AliRecInfoMaker()
 391 {
 392   //
 393   // Destructor
 394   //
 395   if (fLoader) {
 396     delete fLoader;
 397   }
 398 }
 399
 400 //////////////////////////////////////////////////////////////
 401 Int_t AliRecInfoMaker::SetIO()
 402 {
 403   //
 404   // SetIO  - Create the input trees
 405   //
 406   CreateTreeCmp();
 407   if (!fTreeCmp) return 1;
 408   fParamTPC = GetTPCParam();
 409   //
 410   if (!ConnectGenTree()) {
 411     cerr<<"Cannot connect tree with generated tracks"<<endl;
 412     return 1;
 413   }
 414   return 0;
 415 }
 416
 417 //////////////////////////////////////////////////////////////
 418
 419 Int_t AliRecInfoMaker::SetIO(Int_t eventNr)
 420 {
 421   //
 422   //
 423   // SET INPUT
 424   //
 425   TFile f("AliESDs.root");
 426   //
 427
 428   TTree* tree = (TTree*) f.Get("esdTree");
 429   tree->SetBranchStatus("*",1);
 430   fEvent = new AliESDEvent;
 431
 432   if (tree->GetBranch("ESD")){
 433     //    tree->SetBranchAddress("ESD", &fEvent);
 434     // tree->SetBranchAddress("ESDfriend.",&fESDfriend);
 435     // tree->GetEntry(eventNr);
 436     // fEvent->SetESDfriend(fESDfriend);
 437   }else{
 438     fEvent->ReadFromTree(tree);
 439     fESDfriend = (AliESDfriend*)fEvent->FindListObject("AliESDfriend");
 440     tree->GetEntry(eventNr);
 441     fEvent->SetESDfriend(fESDfriend);
 442   }
 443
 444
 445
 446   if (!fEvent) return 1;
 447
 448   return 0;
 449 }
 450
 451
 452
 453 ////////////////////////////////////////////////////////////////////////
 454 void AliRecInfoMaker::Reset()
 455 {
 456   //
 457   // Reset the class
 458   //
 459   fEventNr = 0;
 460   fNEvents = 0;
 461   fTreeCmp = 0;
 462   fTreeCmpKinks =0;
 463   fTreeCmpV0 =0;
 464   //  fFnCmp = "cmpTracks.root";
 465   fFileGenTracks = 0;
 466   fDebug = 0;
 467   //
 468   fParamTPC = 0;
 469   fEvent =0;
 470 }
 471
 472 ////////////////////////////////////////////////////////////////////////
 473 Int_t AliRecInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
 474 {
 475   //
 476   // Exec comparison for subrange of events
 477   //
 478   fNEvents = nEvents;
 479   fFirstEventNr = firstEventNr;
 480   return Exec();
 481 }
 482
 483 ////////////////////////////////////////////////////////////////////////
 484 Int_t AliRecInfoMaker::Exec()
 485 {
 486   //
 487   // Exec comparison
 488   //
 489   TStopwatch timer;
 490   timer.Start();
 491
 492   if (SetIO()==1)
 493     return 1;
 494
 495   fNextTreeGenEntryToRead = 0;
 496   fNextKinkToRead = 0;
 497   fNextV0ToRead   =0;
 498   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
 499   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
 500        eventNr++) {
 501     fEventNr = eventNr;
 502     SetIO(fEventNr);
 503     fNParticles = gAlice->GetEvent(fEventNr);
 504
 505     fIndexRecTracks = new Short_t[fNParticles*20];  //write at maximum 4 tracks corresponding to particle
 506     fIndexRecKinks  = new Short_t[fNParticles*20];  //write at maximum 20 tracks corresponding to particle
 507     fIndexRecV0  = new Short_t[fNParticles*20];  //write at maximum 20 tracks corresponding to particle
 508
 509     fFakeRecTracks = new Short_t[fNParticles];
 510     fMultiRecTracks = new Short_t[fNParticles];
 511     fMultiRecKinks = new Short_t[fNParticles];
 512     fMultiRecV0 = new Short_t[fNParticles];
 513
 514     for (Int_t i = 0; i<fNParticles; i++) {
 515       for (Int_t j=0;j<20;j++){
 516         fIndexRecTracks[20*i+j] = -1;
 517         fIndexRecKinks[20*i+j]  = -1;
 518         fIndexRecV0[20*i+j]  = -1;
 519       }
 520       fFakeRecTracks[i] = 0;
 521       fMultiRecTracks[i] = 0;
 522       fMultiRecKinks[i] = 0;
 523       fMultiRecV0[i] = 0;
 524     }
 525
 526     cout<<"Start to process event "<<fEventNr<<endl;
 527     cout<<"\tfNParticles = "<<fNParticles<<endl;
 528     if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
 529     if (TreeTLoop()>0) return 1;
 530
 531     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
 532     if (TreeGenLoop(eventNr)>0) return 1;
 533     //BuildKinkInfo0(eventNr);
 534     BuildV0Info(eventNr); // no V0 info for a moment
 535     fRecArray->Delete();
 536
 537     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
 538
 539     delete [] fIndexRecTracks;
 540     delete [] fIndexRecKinks;
 541     delete [] fIndexRecV0;
 542     delete [] fFakeRecTracks;
 543     delete [] fMultiRecTracks;
 544     delete [] fMultiRecKinks;
 545     delete [] fMultiRecV0;
 546   }
 547
 548   CloseOutputFile();
 549
 550   cerr<<"Exec finished"<<endl;
 551   timer.Stop();
 552   timer.Print();
 553   return 0;
 554
 555 }
 556 ////////////////////////////////////////////////////////////////////////
 557 Bool_t AliRecInfoMaker::ConnectGenTree()
 558 {
 559 //
 560 // connect all branches from the genTracksTree
 561 // use the same variables as for the new cmp tree, it may work
 562 //
 563   fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
 564   if (!fFileGenTracks) {
 565     cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
 566     return kFALSE;
 567   }
 568   fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
 569   if (!fTreeGenTracks) {
 570     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
 571         <<fFnGenTracks<<endl;
 572     return kFALSE;
 573   }
 574   //
 575   fMCInfo = new AliMCInfo;
 576   fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
 577   //
 578   //
 579   fTreeGenKinks = (TTree*)fFileGenTracks->Get("genKinksTree");
 580   if (!fTreeGenKinks) {
 581     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
 582         <<fFnGenTracks<<endl;
 583     //return kFALSE;
 584   }
 585   else{
 586     fGenKinkInfo = new AliGenKinkInfo;
 587     fTreeGenKinks->SetBranchAddress("MC",&fGenKinkInfo);
 588   }
 589
 590   fTreeGenV0 = (TTree*)fFileGenTracks->Get("genV0Tree");
 591   if (!fTreeGenV0) {
 592     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
 593         <<fFnGenTracks<<endl;
 594     //return kFALSE;
 595   }
 596   else{
 597     fGenV0Info = new AliGenV0Info;
 598     fTreeGenV0->SetBranchAddress("MC",&fGenV0Info);
 599   }
 600   //
 601   if (fDebug > 1) {
 602     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
 603   }
 604   return kTRUE;
 605 }
 606
 607
 608 ////////////////////////////////////////////////////////////////////////
 609 void AliRecInfoMaker::CreateTreeCmp()
 610 {
 611   //
 612   // Create file and tree with comparison information
 613   //
 614   fFileCmp = TFile::Open(fFnCmp,"RECREATE");
 615   if (!fFileCmp) {
 616     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
 617     return;
 618   }
 619   //
 620   //
 621   fTreeCmp    = new TTree("ESDcmpTracks","ESDcmpTracks");
 622   fMCInfo = new AliMCInfo;
 623   fRecInfo = new AliESDRecInfo;
 624   AliESDtrack * esdTrack = new AliESDtrack;
 625   //  AliITStrackMI * itsTrack = new AliITStrackMI;
 626   fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo,256000);
 627   fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo,256000);
 628   //  fTreeCmp->Branch("ITS","AliITStrackMI",&itsTrack);
 629   delete esdTrack;
 630   //
 631   //
 632   fTreeCmpKinks    = new TTree("ESDcmpKinks","ESDcmpKinks");
 633   fGenKinkInfo     = new AliGenKinkInfo;
 634   fRecKinkInfo     = new AliESDRecKinkInfo;
 635   fTreeCmpKinks->Branch("MC.","AliGenKinkInfo",&fGenKinkInfo,256000);
 636   fTreeCmpKinks->Branch("RC.","AliESDRecKinkInfo",&fRecKinkInfo,256000);
 637   //
 638   //
 639   fTreeCmpV0       = new TTree("ESDcmpV0","ESDcmpV0");
 640   fGenV0Info     = new AliGenV0Info;
 641   fRecV0Info     = new AliESDRecV0Info;
 642   fTreeCmpV0->Branch("MC.","AliGenV0Info",   &fGenV0Info,256000);
 643   fTreeCmpV0->Branch("RC.","AliESDRecV0Info",&fRecV0Info,256000);
 644   //
 645   fTreeCmp->AutoSave();
 646   fTreeCmpKinks->AutoSave();
 647   fTreeCmpV0->AutoSave();
 648 }
 649
 650 ////////////////////////////////////////////////////////////////////////
 651 void AliRecInfoMaker::CloseOutputFile()
 652 {
 653   //
 654   // Close output file
 655   //
 656
 657   if (!fFileCmp) {
 658     cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
 659     return;
 660   }
 661   fFileCmp->cd();
 662   fTreeCmp->Write();
 663   delete fTreeCmp;
 664
 665   fFileCmp->Close();
 666   delete fFileCmp;
 667   return;
 668 }
 669 ////////////////////////////////////////////////////////////////////////
 670
 671 TVector3 AliRecInfoMaker::TR2Local(AliTrackReference *trackRef,
 672                             AliTPCParam *paramTPC) {
 673
 674   //
 675   // Transform position to the local coord frame
 676   //
 677
 678   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
 679   Int_t index[4];
 680   paramTPC->Transform0to1(x,index);
 681   paramTPC->Transform1to2Ideal(x,index);
 682   return TVector3(x);
 683 }
 684 ////////////////////////////////////////////////////////////////////////
 685
 686 Int_t AliRecInfoMaker::TreeTLoop()
 687 {
 688   //
 689   // loop over all ESD reconstructed tracks and store info in memory
 690   //
 691   // + loop over all reconstructed kinks
 692   TStopwatch  timer;
 693   timer.Start();
 694   //
 695   Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks();
 696   Int_t nKinks = (Int_t) fEvent->GetNumberOfKinks();
 697   Int_t nV0MIs = (Int_t) fEvent->GetNumberOfV0s();
 698   fSignedKinks = new Short_t[nKinks];
 699   fSignedV0    = new Short_t[nV0MIs];
 700   //
 701   // load kinks to the memory
 702   for (Int_t i=0; i<nKinks;i++){
 703     //    AliESDkink * kink =
 704     fEvent->GetKink(i);
 705     fSignedKinks[i]=0;
 706   }
 707   //
 708   for (Int_t i=0; i<nV0MIs;i++){
 709     //AliV0 * v0MI =
 710     (AliV0*)fEvent->GetV0(i);
 711     fSignedV0[i]=0;
 712   }
 713
 714   //
 715   //
 716   AliESDtrack * track=0;
 717   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
 718     track = (AliESDtrack*)fEvent->GetTrack(iEntry);
 719     //
 720     Int_t label = track->GetLabel();
 721     Int_t absLabel = abs(label);
 722     if (absLabel < fNParticles) {
 723       //      fIndexRecTracks[absLabel] =  iEntry;
 724       if (label < 0) fFakeRecTracks[absLabel]++;
 725       if (fMultiRecTracks[absLabel]>0){
 726         if (fMultiRecTracks[absLabel]<20)
 727           fIndexRecTracks[absLabel*20+fMultiRecTracks[absLabel]] =  iEntry;
 728       }
 729       else
 730         fIndexRecTracks[absLabel*20] =  iEntry;
 731       fMultiRecTracks[absLabel]++;
 732     }
 733   }
 734   // sort reconstructed kinks
 735   //
 736   AliESDkink * kink=0;
 737   for (Int_t iEntry=0; iEntry<nKinks;iEntry++){
 738     kink = (AliESDkink*)fEvent->GetKink(iEntry);
 739     if (!kink) continue;
 740     //
 741     Int_t label0 = TMath::Abs(kink->GetLabel(0));
 742     Int_t label1 = TMath::Abs(kink->GetLabel(1));
 743     Int_t absLabel = TMath::Min(label0,label1);
 744     if (absLabel < fNParticles) {
 745       if (fMultiRecKinks[absLabel]>0){
 746         if (fMultiRecKinks[absLabel]<20)
 747           fIndexRecKinks[absLabel*20+fMultiRecKinks[absLabel]] =  iEntry;
 748       }
 749       else
 750         fIndexRecKinks[absLabel*20] =  iEntry;
 751       fMultiRecKinks[absLabel]++;
 752     }
 753   }
 754   // --sort reconstructed V0
 755   //
 756   AliV0 * v0MI=0;
 757   for (Int_t iEntry=0; iEntry<nV0MIs;iEntry++){
 758     v0MI = (AliV0*)fEvent->GetV0(iEntry);
 759     if (!v0MI) continue;
 760     //
 761     //
 762     //
 763     //Int_t label0 = TMath::Abs(v0MI->GetLabel(0));
 764     //Int_t label1 = TMath::Abs(v0MI->GetLabel(1));
 765     AliESDtrack * trackn = fEvent->GetTrack((v0MI->GetNindex()));
 766     AliESDtrack * trackp = fEvent->GetTrack((v0MI->GetPindex()));
 767     Int_t labels[2]={-1,-1};
 768     labels[0] = (trackn==0) ? -1 : TMath::Abs(trackn->GetLabel());
 769     labels[1] = (trackp==0) ? -1 : TMath::Abs(trackp->GetLabel());
 770     //
 771     for (Int_t i=0;i<2;i++){
 772       Int_t absLabel =  TMath::Abs(labels[i]);
 773       if (absLabel < fNParticles) {
 774         if (fMultiRecV0[absLabel]>0){
 775           if (fMultiRecV0[absLabel]<20)
 776             fIndexRecV0[absLabel*20+fMultiRecV0[absLabel]] =  iEntry;
 777         }
 778         else
 779           fIndexRecV0[absLabel*20] =  iEntry;
 780         fMultiRecV0[absLabel]++;
 781       }
 782     }
 783   }
 784
 785
 786   printf("Time spended in TreeTLoop\n");
 787   timer.Print();
 788
 789   if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
 790   return 0;
 791 }
 792
 793 ////////////////////////////////////////////////////////////////////////
 794 Int_t AliRecInfoMaker::TreeGenLoop(Int_t eventNr)
 795 {
 796 //
 797 // loop over all entries for a given event, find corresponding
 798 // rec. track and store in the fTreeCmp
 799 //
 800   TStopwatch timer;
 801   timer.Start();
 802   Int_t entry = fNextTreeGenEntryToRead;
 803   Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
 804   cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
 805       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
 806   TBranch * branch = fTreeCmp->GetBranch("RC");
 807   //  TBranch * branchF = fTreeCmp->GetBranch("F");
 808
 809   branch->SetAddress(&fRecInfo); // set all pointers
 810   fRecArray = new TObjArray(fNParticles);
 811   AliESDtrack dummytrack;  //
 812   AliESDfriendTrack dummytrackF;  //
 813
 814   while (entry < nParticlesTR) {
 815     fTreeGenTracks->GetEntry(entry);
 816     entry++;
 817     if (eventNr < fMCInfo->fEventNr) continue;
 818     if (eventNr > fMCInfo->fEventNr) continue;
 819     if (fMCInfo->GetCharge()==0) continue;
 820     //
 821     fNextTreeGenEntryToRead = entry-1;
 822     if (fDebug > 2 && fMCInfo->fLabel < 10) {
 823       cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
 824     }
 825     //    if (fMCInfo->fNTPCRef<1) continue; // not TPCref
 826     //
 827     fRecInfo->Reset();
 828     AliESDtrack * track=0;
 829     fRecInfo->fReconstructed =0;
 830     TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
 831     local.GetXYZ(fRecInfo->fTRLocalCoord);
 832     //
 833     if (fIndexRecTracks[fMCInfo->fLabel*20] >= 0) {
 834       track= (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20]);
 835       if(!track) continue;
 836       //
 837       //
 838       // find nearest track if multifound
 839       //Int_t sign = Int_t(track->GetSign()*fMCInfo->fCharge);
 840       //
 841       Int_t status = 0;
 842       if  ((track->GetStatus()&AliESDtrack::kITSrefit)>0) status++;
 843       if  ((track->GetStatus()&AliESDtrack::kTPCrefit)>0) status++;
 844       if  ((track->GetStatus()&AliESDtrack::kTRDrefit)>0) status++;
 845
 846       //
 847       if (fIndexRecTracks[fMCInfo->fLabel*20+1]>0){
 848         //
 849         Double_t p[3];
 850         track->GetInnerPxPyPz(p);
 851         Float_t maxp = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
 852         //
 853         for (Int_t i=1;i<20;i++){
 854           if (fIndexRecTracks[fMCInfo->fLabel*20+i]>=0){
 855             AliESDtrack * track2 = (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20+i]);
 856             if (!track2) continue;
 857             //Int_t sign2 = track->GetSign()*fMCInfo->fCharge; //
 858             //if (sign2<0) continue;
 859             track2->GetInnerPxPyPz(p);
 860             Float_t mom = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
 861             /*
 862             if (sign<0){
 863               sign = sign2;
 864               track = track2;
 865               maxp = mom;
 866               continue;
 867             }
 868             */
 869             //
 870             Int_t status2 = 0;
 871             if  ((track2->GetStatus()&AliESDtrack::kITSrefit)>0) status2++;
 872             if  ((track2->GetStatus()&AliESDtrack::kTPCrefit)>0) status2++;
 873             if  ((track2->GetStatus()&AliESDtrack::kTRDrefit)>0) status2++;
 874             if (status2<status) continue;
 875             //
 876             if (mom<maxp) continue;
 877             maxp = mom;
 878             track = track2;
 879             //
 880           }
 881         }
 882       }
 883       //
 884       fRecInfo->SetESDtrack(track);
 885
 886       /*
 887       if (track) {
 888         fRecInfo->SetESDtrack(track);
 889       }else{
 890         fRecInfo->SetESDtrack(&dummytrack);
 891       }
 892       */
 893       //
 894
 895       fRecInfo->fReconstructed = 1;
 896       fRecInfo->fFake     = fFakeRecTracks[fMCInfo->fLabel];
 897       fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
 898       //
 899       fRecInfo->Update(fMCInfo,fParamTPC,kTRUE);
 900     }
 901     else{
 902       fRecInfo->SetESDtrack(&dummytrack);
 903       fRecInfo->Update(fMCInfo,fParamTPC,kFALSE);
 904     }
 905     fRecArray->AddAt(new AliESDRecInfo(*fRecInfo),fMCInfo->fLabel);
 906     fTreeCmp->Fill();
 907   }
 908   fTreeCmp->AutoSave();
 909   printf("Time spended in TreeGenLoop\n");
 910   timer.Print();
 911   if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
 912
 913   return 0;
 914 }
 915
 916
 917
 918 ////////////////////////////////////////////////////////////////////////
 919 ////////////////////////////////////////////////////////////////////////
 920 ////////////////////////////////////////////////////////////////////////
 921 Int_t AliRecInfoMaker::BuildKinkInfo0(Int_t eventNr)
 922 {
 923 //
 924 // loop over all entries for a given event, find corresponding
 925 // rec. track and store in the fTreeCmp
 926 //
 927   TStopwatch timer;
 928   timer.Start();
 929   Int_t entry = fNextKinkToRead;
 930   Double_t nParticlesTR = fTreeGenKinks->GetEntriesFast();
 931   cerr<<"fNParticles, nParticlesTR, fNextKinkToRead: "<<fNParticles<<" "
 932       <<nParticlesTR<<" "<<fNextKinkToRead<<endl;
 933   //
 934   TBranch * branch = fTreeCmpKinks->GetBranch("RC.");
 935   branch->SetAddress(&fRecKinkInfo); // set all pointers
 936
 937   //
 938   while (entry < nParticlesTR) {
 939     fTreeGenKinks->GetEntry(entry);
 940     entry++;
 941     if (eventNr < fGenKinkInfo->GetMinus().fEventNr) continue;
 942     if (eventNr > fGenKinkInfo->GetMinus().fEventNr) continue;;
 943     //
 944     fNextKinkToRead = entry-1;
 945     //
 946     //
 947     AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->GetMinus().fLabel);
 948     if(!fRecInfo1) continue;
 949     AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->GetPlus().fLabel);
 950     if(!fRecInfo2) continue;
 951     fRecKinkInfo->fT1 = (*fRecInfo1);
 952     fRecKinkInfo->fT2 = (*fRecInfo2);
 953     fRecKinkInfo->fStatus =0;
 954     if (fRecInfo1 && fRecInfo1->fTPCOn) fRecKinkInfo->fStatus+=1;
 955     if (fRecInfo2 && fRecInfo2->fTPCOn) fRecKinkInfo->fStatus+=2;
 956     if (fRecKinkInfo->fStatus==3&&fRecInfo1->fSign!=fRecInfo2->fSign) fRecKinkInfo->fStatus*=-1;
 957
 958     if (fRecKinkInfo->fStatus==3){
 959       fRecKinkInfo->Update();
 960     }
 961     Int_t label =  TMath::Min(fGenKinkInfo->GetMinus().fLabel,fGenKinkInfo->GetPlus().fLabel);
 962     Int_t label2 = TMath::Max(fGenKinkInfo->GetMinus().fLabel,fGenKinkInfo->GetPlus().fLabel);
 963
 964     AliESDkink *kink=0;
 965     fRecKinkInfo->fRecStatus   =0;
 966     fRecKinkInfo->fMultiple    = fMultiRecKinks[label];
 967     fRecKinkInfo->fKinkMultiple=0;
 968     //
 969     if (fMultiRecKinks[label]>0){
 970
 971       //      for (Int_t j=0;j<TMath::Min(fMultiRecKinks[label],100);j++){
 972       for (Int_t j=TMath::Min(fMultiRecKinks[label],Short_t(20))-1;j>=0;j--){
 973         Int_t index = fIndexRecKinks[label*20+j];
 974         //AliESDkink *kink2  = (AliESDkink*)fKinks->At(index);
 975         AliESDkink *kink2  = (AliESDkink*)fEvent->GetKink(index);
 976         if (TMath::Abs(kink2->GetLabel(0))==label &&TMath::Abs(kink2->GetLabel(1))==label2) {
 977           fRecKinkInfo->fKinkMultiple++;
 978           fSignedKinks[index]=1;
 979           Int_t c0=0;
 980           if (kink){
 981             //      if (kink->fTRDOn) c0++;
 982             //if (kink->fITSOn) c0++;
 983             if (kink->GetStatus(2)>0) c0++;
 984             if (kink->GetStatus(0)>0) c0++;
 985           }
 986           Int_t c2=0;
 987           //      if (kink2->fTRDOn) c2++;
 988           //if (kink2->fITSOn) c2++;
 989           if (kink2->GetStatus(2)>0) c2++;
 990           if (kink2->GetStatus(0)>0) c2++;
 991
 992           if (c2<c0) continue;
 993           kink =kink2;
 994         }
 995         if (TMath::Abs(kink2->GetLabel(1))==label &&TMath::Abs(kink2->GetLabel(0))==label2) {
 996           fRecKinkInfo->fKinkMultiple++;
 997           fSignedKinks[index]=1;
 998           Int_t c0=0;
 999           if (kink){
1000             //if (kink->fTRDOn) c0++;
1001             //if (kink->fITSOn) c0++;
1002             if (kink->GetStatus(2)>0) c0++;
1003             if (kink->GetStatus(0)>0) c0++;
1004
1005           }
1006           Int_t c2=0;
1007           //      if (kink2->fTRDOn) c2++;
1008           //if (kink2->fITSOn) c2++;
1009           if (kink2->GetStatus(2)>0) c2++;
1010           if (kink2->GetStatus(0)>0) c2++;
1011
1012           if (c2<c0) continue;
1013           kink =kink2;
1014         }
1015       }
1016     }
1017     if (kink){
1018       fRecKinkInfo->fKink = *kink;
1019       fRecKinkInfo->fRecStatus=1;
1020     }
1021     fTreeCmpKinks->Fill();
1022   }
1023   //  Int_t nkinks = fKinks->GetEntriesFast();
1024   Int_t nkinks = fEvent->GetNumberOfKinks();
1025   for (Int_t i=0;i<nkinks;i++){
1026     if (fSignedKinks[i]==0){
1027       //      AliESDkink *kink  = (AliESDkink*)fKinks->At(i);
1028       AliESDkink *kink  = (AliESDkink*)fEvent->GetKink(i);
1029       if (!kink) continue;
1030       //
1031       fRecKinkInfo->fKink = *kink;
1032       fRecKinkInfo->fRecStatus =-2;
1033       //
1034       AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->GetLabel(0)));
1035       AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->GetLabel(1)));
1036       if (fRecInfo1 && fRecInfo2){
1037         fRecKinkInfo->fT1 = (*fRecInfo1);
1038         fRecKinkInfo->fT2 = (*fRecInfo2);
1039         fRecKinkInfo->fRecStatus =-1;
1040       }
1041       fTreeCmpKinks->Fill();
1042     }
1043   }
1044
1045
1046   fTreeCmpKinks->AutoSave();
1047   printf("Time spended in BuilKinkInfo Loop\n");
1048   timer.Print();
1049   if (fDebug > 2) cerr<<"end of BuildKinkInfo Loop"<<endl;
1050   return 0;
1051 }
1052
1053
1054
1055 ////////////////////////////////////////////////////////////////////////
1056 ////////////////////////////////////////////////////////////////////////
1057 ////////////////////////////////////////////////////////////////////////
1058
1059
1060
1061 Int_t AliRecInfoMaker::BuildV0Info(Int_t eventNr)
1062 {
1063 //
1064 // loop over all entries for a given event, find corresponding
1065 // rec. track and store in the fTreeCmp
1066 //
1067   static TDatabasePDG pdgtable;
1068
1069   TStopwatch timer;
1070   timer.Start();
1071   Int_t entry = fNextV0ToRead;
1072   Double_t nParticlesTR = fTreeGenV0->GetEntriesFast();
1073   cerr<<"fNParticles, nParticlesTR, fNextV0ToRead: "<<fNParticles<<" "
1074       <<nParticlesTR<<" "<<fNextV0ToRead<<endl;
1075   //
1076   TBranch * branch = fTreeCmpV0->GetBranch("RC.");
1077   branch->SetAddress(&fRecV0Info); // set all pointers
1078   const AliESDVertex * esdvertex = fEvent->GetVertex();
1079   Float_t vertex[3]= {esdvertex->GetXv(), esdvertex->GetYv(),esdvertex->GetZv()};
1080
1081   //
1082   while (entry < nParticlesTR) {
1083     fTreeGenV0->GetEntry(entry);
1084     entry++;
1085     fRecV0Info->Reset();  //reset all variables
1086     if (eventNr < fGenV0Info->GetMinus().fEventNr) continue;
1087     if (eventNr > fGenV0Info->GetMinus().fEventNr) continue;;
1088     //
1089     fNextV0ToRead = entry-1;
1090     //
1091     //
1092     AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->GetMinus().fLabel);
1093     AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->GetPlus().fLabel);
1094     if (fGenV0Info->GetMinus().fCharge*fGenV0Info->GetPlus().fCharge>0) continue;  // interactions
1095     if (!fRecInfo1 || !fRecInfo2) continue;
1096     fRecV0Info->fT1 = (*fRecInfo1);
1097     fRecV0Info->fT2 = (*fRecInfo2);
1098     fRecV0Info->fV0Status =0;
1099     if (fRecInfo1 && fRecInfo1->fStatus[1]>0) fRecV0Info->fV0Status+=1;
1100     if (fRecInfo2 && fRecInfo2->fStatus[1]>0) fRecV0Info->fV0Status+=2;
1101
1102     if (fRecV0Info->fV0Status==3&&fRecInfo1->fSign==fRecInfo2->fSign) fRecV0Info->fV0Status*=-1;
1103
1104
1105     if (abs(fRecV0Info->fV0Status)==3){
1106       fRecV0Info->Update(vertex);
1107       {
1108         //
1109         // TPC V0 Info
1110         Double_t x,alpha, param[5],cov[15];
1111         if ( fRecV0Info->fT1.GetESDtrack()->GetInnerParam() && fRecV0Info->fT2.GetESDtrack()->GetInnerParam()){
1112           fRecV0Info->fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
1113           fRecV0Info->fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
1114           AliExternalTrackParam paramP(x,alpha,param,cov);
1115           //
1116           fRecV0Info->fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
1117           fRecV0Info->fT2.GetESDtrack()->GetInnerExternalCovariance(cov);
1118           AliExternalTrackParam paramM(x,alpha,param,cov);
1119           //
1120           fRecV0Info->fV0tpc->SetParamN(paramM);
1121           fRecV0Info->fV0tpc->SetParamP(paramP);
1122           Double_t pid1[5],pid2[5];
1123           fRecV0Info->fT1.GetESDtrack()->GetESDpid(pid1);
1124           fRecV0Info->fT1.GetESDtrack()->GetESDpid(pid2);
1125           //
1126           //fRecV0Info->fV0tpc.UpdatePID(pid1,pid2);
1127           fRecV0Info->fV0tpc->Update(vertex);
1128
1129           //
1130           //
1131           fRecV0Info->fT1.GetESDtrack()->GetExternalParameters(x,param);
1132           fRecV0Info->fT1.GetESDtrack()->GetExternalCovariance(cov);
1133           alpha = fRecV0Info->fT1.GetESDtrack()->GetAlpha();
1134           new (&paramP) AliExternalTrackParam(x,alpha,param,cov);
1135           //
1136           fRecV0Info->fT2.GetESDtrack()->GetExternalParameters(x,param);
1137           fRecV0Info->fT2.GetESDtrack()->GetExternalCovariance(cov);
1138           alpha = fRecV0Info->fT2.GetESDtrack()->GetAlpha();
1139           new (&paramM) AliExternalTrackParam(x,alpha,param,cov);
1140           //
1141           fRecV0Info->fV0its->SetParamN(paramM);
1142           fRecV0Info->fV0its->SetParamP(paramP);
1143           //      fRecV0Info->fV0its.UpdatePID(pid1,pid2);
1144           fRecV0Info->fV0its->Update(vertex);
1145         }
1146       }
1147       //
1148       // ????
1149       //
1150       if (TMath::Abs(fGenV0Info->GetMinus().fPdg)==11 &&TMath::Abs(fGenV0Info->GetPlus().fPdg)==11){
1151         if (fRecV0Info->fDist2>10){
1152           fRecV0Info->Update(vertex);
1153         }
1154         if (fRecV0Info->fDist2>10){
1155           fRecV0Info->Update(vertex);
1156         }
1157       }
1158     }
1159     //
1160     // take the V0 from reconstruction
1161
1162     Int_t label =  TMath::Min(fGenV0Info->GetMinus().fLabel,fGenV0Info->GetPlus().fLabel);
1163     Int_t label2 = TMath::Max(fGenV0Info->GetMinus().fLabel,fGenV0Info->GetPlus().fLabel);
1164     AliV0 *v0MI=0;
1165     AliV0 *v0MIOff=0;
1166     fRecV0Info->fRecStatus   =0;
1167     fRecV0Info->fMultiple    = fMultiRecV0[label];
1168     fRecV0Info->fV0MultipleOn=0;
1169     fRecV0Info->fV0MultipleOff=0;
1170     //
1171     if (fMultiRecV0[label]>0 || fMultiRecV0[label2]>0){
1172
1173       //      for (Int_t j=0;j<TMath::Min(fMultiRecV0s[label],100);j++){
1174       for (Int_t j=TMath::Min(fMultiRecV0[label],Short_t(20))-1;j>=0;j--){
1175         Int_t index = fIndexRecV0[label*20+j];
1176         if (index<0) continue;
1177         AliV0 *v0MI2  = (AliV0*)fEvent->GetV0(index);
1178         // get track labels
1179         AliESDtrack * trackn = fEvent->GetTrack((v0MI2->GetNindex()));
1180         AliESDtrack * trackp = fEvent->GetTrack((v0MI2->GetPindex()));
1181         Int_t vlabeln = (trackn==0) ? -1 : trackn->GetLabel();
1182         Int_t vlabelp = (trackp==0) ? -1 : trackp->GetLabel();
1183         fRecV0Info->fLab[0]=TMath::Abs(vlabelp);
1184         fRecV0Info->fLab[1]=TMath::Abs(vlabeln);
1185         //
1186         if (TMath::Abs(vlabeln)==label &&TMath::Abs(vlabelp)==label2) {
1187           if (v0MI2->GetOnFlyStatus()) {
1188             v0MI =v0MI2;
1189             fRecV0Info->fV0MultipleOn++;
1190           }else  {
1191             v0MIOff = v0MI2;
1192             fRecV0Info->fV0MultipleOff++;
1193           }
1194           fSignedV0[index]=1;
1195         }
1196         if (TMath::Abs(vlabelp)==label &&TMath::Abs(vlabeln)==label2) {
1197           if (v0MI2->GetOnFlyStatus()){
1198             v0MI =v0MI2;
1199             fRecV0Info->fV0MultipleOn++;
1200           }else  {
1201             v0MIOff = v0MI2;
1202             fRecV0Info->fV0MultipleOff++;
1203           }
1204           fSignedV0[index]=1;
1205         }
1206       }
1207     }
1208     if (v0MI){
1209       new (fRecV0Info->fV0rec) AliV0(*v0MI);
1210       fRecV0Info->fRecStatus=1;
1211     }
1212     if (v0MIOff){
1213       new (fRecV0Info->fV0recOff) AliV0(*v0MIOff);
1214       fRecV0Info->fRecStatus=1;
1215     }
1216     Int_t mpdg = fGenV0Info->GetMother().GetPdgCode();
1217     Float_t mass = ( pdgtable.GetParticle(mpdg)==0) ? 0 :pdgtable.GetParticle(mpdg)->Mass();
1218     fRecV0Info->UpdateKF(*esdvertex,
1219                          fGenV0Info->GetPlus().GetPdg(),
1220                          fGenV0Info->GetMinus().GetPdg(),
1221                          mass);
1222     fTreeCmpV0->Fill();
1223   }
1224   //
1225   // write fake v0s
1226   //
1227   Int_t nV0MIs = fEvent->GetNumberOfV0s();
1228   for (Int_t i=0;i<nV0MIs;i++){
1229     if (fSignedV0[i]==0){
1230       AliV0 *v0MI  = (AliV0*)fEvent->GetV0(i);
1231       if (!v0MI) continue;
1232       fRecV0Info->Reset();  //reset all variables
1233       //
1234       new (fRecV0Info->fV0rec) AliV0(*v0MI);
1235       fRecV0Info->fV0Status  =-10;
1236       fRecV0Info->fRecStatus =-2;
1237       //
1238       AliESDtrack * trackn = fEvent->GetTrack((v0MI->GetNindex()));
1239       AliESDtrack * trackp = fEvent->GetTrack((v0MI->GetPindex()));
1240       Int_t vlabeln = (trackn==0) ? -1 : trackn->GetLabel();
1241       Int_t vlabelp = (trackp==0) ? -1 : trackp->GetLabel();
1242       fRecV0Info->fLab[0]=TMath::Abs(vlabelp);
1243       fRecV0Info->fLab[1]=TMath::Abs(vlabeln);
1244       if (TMath::Abs(fRecV0Info->fLab[0] - fRecV0Info->fLab[1])<2) continue;
1245       AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(vlabeln));
1246       AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(vlabelp));
1247       if (fRecInfo1 && fRecInfo2){
1248         fRecV0Info->fT1 = (*fRecInfo1);
1249         fRecV0Info->fT2 = (*fRecInfo2);
1250         fRecV0Info->fRecStatus =-1;
1251       }
1252       fRecV0Info->Update(vertex);
1253       fRecV0Info->UpdateKF(*esdvertex,211,211,0.49767);
1254       fTreeCmpV0->Fill();
1255     }
1256   }
1257
1258
1259
1260   fTreeCmpV0->AutoSave();
1261   printf("Time spended in BuilV0Info Loop\n");
1262   timer.Print();
1263   if (fDebug > 2) cerr<<"end of BuildV0Info Loop"<<endl;
1264   return 0;
1265 }
1266 ////////////////////////////////////////////////////////////////////////
1267 ////////////////////////////////////////////////////////////////////////
1268
1269