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