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