void runLumi3myTrackFit(const int nEvents=10, TString storePath="tmpOutput", const int verboseLevel=0, double dv=0.5, bool no4d=false)
{
  // ----  Load libraries   -------------------------------------------------
  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
  gROOT->LoadMacro("line3Dfit.C");
  // ------------------------------------------------------------------------

  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------

  TChain t("cbmsim");
  t.Add((storePath+"/Lumi_reco.root"));

  TClonesArray* reco_array=new TClonesArray("PndMvdHit");
  t.SetBranchAddress("MVDHitsStrip",&reco_array);  //Points for Tracks

  TFile *f;
  if(no4d){
    f = new TFile(storePath+"/Lumi_track_no4d.root","RECREATE");
  }else{
    f = new TFile(storePath+"/Lumi_track.root","RECREATE");
  }

  Double_t theta;
  Double_t accuracy;
  Double_t direction[3];
  Double_t startpoint[3];

  TTree *tree = new TTree("tracks","reconstructed Tracks");
  tree->Branch("theta",&theta,"theta/D");
  tree->Branch("accuracy",&accuracy,"accuracy/D");
  tree->Branch("direction",&direction,"x/D:y/D:z/D");
  tree->Branch("startpoint",&startpoint,"x/D:y/D:z/D");

  for (Int_t j=0; j<nEvents; j++)
  {
    if(verboseLevel>0) cout << "Event: "<< j <<endl;
    t.GetEntry(j);

//fit pairs of first and second disc events
    std::vector<int> hitsd1, hitsd2, hitsd3, hitsd4; //hit'ids splitted by detectorplane
    std::vector<double> trackStartx, trackStarty, trackStartz, trackVecx, trackVecy, trackVecz;//pseudo tracks
    std::vector<double> trackStartdx, trackStartdy, trackStartdz, trackVecdx, trackVecdy, trackVecdz;//save errors 
    int nHits = reco_array->GetEntriesFast();

    if(nHits<3)
      continue;
    
    double z;
    //check whitch disc was hit
    for (Int_t i=0; i<nHits; i++)
    {
      PndMvdHit *hit=(PndMvdHit*)reco_array->At(i);
      z = hit->GetZ();
      if( (z<1071)&&(z>1069) )  //TODO: just use z koordinate without knowing where the planes are
        hitsd1.push_back(i);
      else if( (z<1121)&&(z>1119) )
        hitsd2.push_back(i);
      else if( (z<1171)&&(z>1169) )
        hitsd3.push_back(i);
      else if( (z<1221)&&(z>1219) )
        hitsd4.push_back(i);
    }//end of hits

    if(verboseLevel>1) cout << "Hits: "<< nHits <<" DiscHits: "<< hitsd1.size()+hitsd2.size()+hitsd3.size()+hitsd4.size() <<endl;

    //iterate first discs, fit with all seconds, save pseudo-tracks
    TVector3 start, tmp, vec, dstart, dvec; //temp-vars
    for (Int_t i=0; i<hitsd1.size(); i++)
    {
      PndMvdHit *hit1=(PndMvdHit*)reco_array->At(hitsd1.at(i));
      start.SetXYZ(hit1->GetX(), hit1->GetY(), hit1->GetZ());
      
      for (Int_t k=0; k<hitsd2.size(); k++)
      {
        PndMvdHit *hit2=(PndMvdHit*)reco_array->At(hitsd2.at(k));
        tmp.SetXYZ(hit2->GetX(), hit2->GetY(), hit2->GetZ());
        vec = tmp - start;   //calc direction vector for FINDING
        trackStartx.push_back(start.x());
        trackStarty.push_back(start.y());
        trackStartz.push_back(start.z());
        trackStartdx.push_back(hit1->GetDx());
        trackStartdy.push_back(hit1->GetDy());
        trackStartdz.push_back(hit1->GetDz());
        trackVecx.push_back(vec.x());        //save vector from start to second
        trackVecy.push_back(vec.y());
        trackVecz.push_back(vec.z());
        trackVecdx.push_back(hit2->GetDx()); //save error of second point for FIT
        trackVecdy.push_back(hit2->GetDy()); //NOT error of direction vector
        trackVecdz.push_back(hit2->GetDz());
      }//end of disc 2 hits
    }//end of disc 1 hits

    if(verboseLevel>1) cout << "Pseudos: "<< trackStartx.size() <<endl;

    //check if disc3 or disc4 has hits in track, add points for fitting
    TVector3 disc3hit, disc4hit, disc3hitD, disc4hitD; //temp-vars
    std::vector<double> fitpointsx, fitpointsy, fitpointsz;//points for fit
    std::vector<double> fitpointsdx, fitpointsdy, fitpointsdz;//errors for fit
    std::vector<int> numPtsTrack;//if more then one track is saven in fit-array
    fitpointsx.clear(); fitpointsy.clear(); fitpointsz.clear();
    for (Int_t i=0; i<trackStartx.size(); i++)
    {
      int pntcnt=2;
      start.SetXYZ(trackStartx.at(i), trackStarty.at(i), trackStartz.at(i));
      dstart.SetXYZ(trackStartdx.at(i), trackStartdy.at(i), trackStartdz.at(i));
      vec.SetXYZ(trackVecx.at(i), trackVecy.at(i), trackVecz.at(i));
      dvec.SetXYZ(trackVecdx.at(i), trackVecdy.at(i), trackVecdz.at(i));

      bool d3hit=false, d4hit=false;
      for (Int_t k=0; k<hitsd3.size(); k++)
      {
        PndMvdHit *hit3=(PndMvdHit*)reco_array->At(hitsd3.at(k));
        disc3hit.SetXYZ(hit3->GetX(), hit3->GetY(), hit3->GetZ());
        disc3hitD.SetXYZ(hit3->GetDx(), hit3->GetDy(), hit3->GetDz());
        tmp = start + vec + vec; //find hit in disc3
        if( abs(tmp.x()-disc3hit.x())<dv && abs(tmp.y()-disc3hit.y())<dv){
          d3hit=true;
          pntcnt++;
          k=hitsd3.size();
        }
      }//end of disc 3 hits

      if(verboseLevel>2 && d3hit) cout << "Disc 3 Hit! "<<endl;

      for (Int_t k=0; k<hitsd4.size(); k++)
      {
        PndMvdHit *hit4=(PndMvdHit*)reco_array->At(hitsd4.at(k));
        disc4hit.SetXYZ(hit4->GetX(), hit4->GetY(), hit4->GetZ());
        disc4hitD.SetXYZ(hit4->GetDx(), hit4->GetDy(), hit4->GetDz());
        tmp = start + vec + vec + vec; //find hit in disc4
        if( abs(tmp.x()-disc4hit.x())<(2*dv) && abs(tmp.y()-disc4hit.y())<(2*dv)){
          d4hit=true;
          pntcnt++;
          k=hitsd4.size();
        }
      }//end of disc 4 hits

      //for testing just 3 planes
      if(no4d && d4hit){
        d4hit=false;
        pntcnt--;
      }  

      if(verboseLevel>2 && d4hit) cout << "Disc 4 Hit! "<<endl;
      
      //save track für fitting
      if(d3hit || d4hit){ //third hit found <=> track found
        //first point with errors
        fitpointsx.push_back(start.x()); fitpointsy.push_back(start.y()); fitpointsz.push_back(start.z());
        fitpointsdx.push_back(dstart.x()); fitpointsdy.push_back(dstart.y()); fitpointsdz.push_back(dstart.z());
        //second point (start+vec because vec:=second-start) with perviously saved errors 
        fitpointsx.push_back(start.x()+vec.x()); fitpointsy.push_back(start.y()+vec.y()); fitpointsz.push_back(start.z()+vec.z());
        fitpointsdx.push_back(dvec.x()); fitpointsdy.push_back(dvec.y()); fitpointsdz.push_back(dvec.z());
        if(d3hit){
          fitpointsx.push_back(disc3hit.x()); fitpointsy.push_back(disc3hit.y()); fitpointsz.push_back(disc3hit.z());
          fitpointsdx.push_back(disc3hitD.x()); fitpointsdy.push_back(disc3hitD.y()); fitpointsdz.push_back(disc3hitD.z());
        }
        if(d4hit){
          fitpointsx.push_back(disc4hit.x()); fitpointsy.push_back(disc4hit.y()); fitpointsz.push_back(disc4hit.z());
          fitpointsdx.push_back(disc4hitD.x()); fitpointsdy.push_back(disc4hitD.y()); fitpointsdz.push_back(disc4hitD.z());
        }
        numPtsTrack.push_back(pntcnt);//save number of points in this track
      }
    }//end of pseudo-tracks

    if(verboseLevel>2) cout << "VecSizes: "<< fitpointsx.size()<< " " << fitpointsy.size()<< " " << fitpointsz.size() <<endl;
    if(verboseLevel>2) cout << "DVecSizes: "<< fitpointsdx.size()<< " " << fitpointsdy.size()<< " " << fitpointsdz.size() <<endl;
    if(verboseLevel>1) cout << "Points Taken: "<< fitpointsx.size()<< " of #Tracks: "<< numPtsTrack.size() <<endl;

    //straight line fit
    TGraph2DErrors* fitme;
    int pt=0;//point counter, loop over tracks, not points!
    for(int track=0; track<numPtsTrack.size(); track++)
    {
      const int numPts = numPtsTrack.at(track); //read how many points in this track
      if(verboseLevel>2) cout << "Track: "<< track<< " Points: "<< numPts <<endl;
      if(numPts > fitpointsx.size()){
        cout << "Error: To less points for tracks!!" << endl;
        track=numPtsTrack.size();
      }

      fitme = new TGraph2DErrors(numPts); //new graph for fitting  

      fitme->SetPoint(0, fitpointsx.at(pt), fitpointsy.at(pt), fitpointsz.at(pt));
      fitme->SetPointError(0, fitpointsdx.at(pt), fitpointsdy.at(pt), fitpointsdz.at(pt));
      fitme->SetPoint(1, fitpointsx.at(pt+1), fitpointsy.at(pt+1), fitpointsz.at(pt+1));
      fitme->SetPointError(1, fitpointsdx.at(pt+1), fitpointsdy.at(pt+1), fitpointsdz.at(pt+1));
      fitme->SetPoint(2, fitpointsx.at(pt+2), fitpointsy.at(pt+2), fitpointsz.at(pt+2));
      fitme->SetPointError(2, fitpointsdx.at(pt+2), fitpointsdy.at(pt+2), fitpointsdz.at(pt+2));
      if(numPts==4){
        fitme->SetPoint(3, fitpointsx.at(pt+3), fitpointsy.at(pt+3), fitpointsz.at(pt+3));
        fitme->SetPointError(3, fitpointsdx.at(pt+3), fitpointsdy.at(pt+3), fitpointsdz.at(pt+3));
      }
      pt+=numPts;

      Double_t parFit[4]; //fit-parameter
      accuracy = line3Dfit(numPts, fitme, parFit);

      //---- save Track as startpoint & direction, also theta and Chi² as error
      startpoint[0]=parFit[0]; startpoint[1]=parFit[2]; startpoint[2]=0.;//save fit as startpoint at z=0
      direction[0]=parFit[1]; direction[1]=parFit[3]; direction[2]=1.;   //and a direction vector
      theta = fabs(atan2(sqrt(direction[0]*direction[0]+direction[1]*direction[1]),fabs(direction[2])));
      if(verboseLevel>0) cout << "Tracked Theta /mrad " << theta*1000 << endl;

      tree->Fill(); //save results
      delete fitme;
    }

    if(verboseLevel>1) cout<<endl;

  }//end of Event's

  tree->Write(tree->GetName(),TObject::kOverwrite);
  f->Close();

  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  cout << endl << endl;
  cout << "Macro finished succesfully." << endl;
  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
  cout << endl;
}