PWG1/AliRecInfoMaker.cxx

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