//ROOT includes
#include <TFile.h>
#include <TLorentzVector.h>
#include <TStopwatch.h>
#include <TChain.h>
#include <TTree.h>
#include <TChainIndex.h>
#include <TTreeIndex.h>
//std-lib
#include <iostream>

//collaborating classes
#include "MatchingTuple.h"
#include "FopiPidHub.h"
#include "FopiPidInfo.h"
#include "TpcEventIdentifier.h"
#include "TpcRefGFTrkResCalc_Alignment.h"
//FOPI2ROOT
#include "CdcTrack.h"
#include "CdcHit.h"
#include "RpcTrack.h"
#include "BarrelTrack.h"
#include "FopiEvent.h"
//Genfit
#include <GFRaveConverters.h>
#include <GFRaveTrackParameters.h>
#include <GFRaveVertex.h>
#include <GFRaveVertexFactory.h>
#include <GFException.h>
#include <GFTrack.h>
#include <GFAbsTrackRep.h>

#include <GFTrack.h>
#include <GFAbsTrackRep.h>
#include <RKTrackRep.h>
//header
#include "CdcTrackChain.h"
#include "../base/TpcdEdx.h"

CdcTrackChain::CdcTrackChain(TChain *chain_, TChain *friendChain_, TString outFile_, bool GFMode):
    fopiEventBranchName("FopiEventBr"),
    tpcEventIdentifierBranchName("TpcEventIdentifier"),
    outFileName(outFile_),
    fVerbose(1),
    fFopiIndex(NULL),
    fGFMode(GFMode)
{
    chain=chain_;
    //friendChain_->Print();
    friendChain=friendChain_;
    //frendChain_->GetEntries();
    cdcHitBranchName="CdcHitBr";
    cdcHitArray=NULL;
    nEvents=0;
    if(GFMode){
        cdcArray=NULL;//new TClonesArray("GFTrack");
        cdcBranchName="TpcTrackPostFit";
        tpcEventIdentifiererArray=NULL;//new TClonesArray("TpcEventIdentifier");
    }else{
        cdcArray=NULL;//new TClonesArray("CdcTrack");
        cdcBranchName="CdcTrackBr";
    }
    arrays.push_back(fopiEventArray);
    fopiEventArray=NULL;//new TClonesArray("FopiEvent");
    arrays.push_back(cdcArray);
    
}


void CdcTrackChain::Init(){
    std::cout<<"CdcTrackChain::Init()"<<std::endl;
    //std::cout<<chain->GetEntries()<<std::endl;
    std::cout<<friendChain->GetEntries()<<std::endl;
    //chain->Print();
    //friendChain->Print();
    chain->SetBranchStatus("*.*",0);
    chain->SetBranchStatus(cdcBranchName+".*",1);
    
    std::cout<<"chain->SetBranchAddress("+cdcBranchName+",&cdcArray);"<<std::endl;
    std::cout<<chain->SetBranchAddress(cdcBranchName,&cdcArray)<<std::endl;
    if(!fGFMode){
        chain->SetBranchStatus(fopiEventBranchName+".*",1);
        chain->SetBranchStatus(cdcHitBranchName+".*",1);
        std::cout<<"chain->SetBranchAddress(fopiEventBranchName,&cfopiEventArray);"<<std::endl;
        std::cout<<chain->SetBranchAddress(fopiEventBranchName,&fopiEventArray)<<std::endl;
        std::cout<<"chain->SetBranchAddress("+cdcHitBranchName+",&cdcHitArray)"<<std::endl;
        std::cout<<chain->SetBranchAddress(cdcHitBranchName,&cdcHitArray)<<std::endl;
    }else{
        chain->SetBranchStatus(tpcEventIdentifierBranchName+".*",1);
        friendChain->SetBranchStatus("*.*",0);
        friendChain->SetBranchStatus(fopiEventBranchName+".*",1);
        std::cout<<"chain->SetBranchAddress(fopiEventBranchName,&cfopiEventArray);"<<std::endl;
        std::cout<<friendChain->SetBranchAddress(fopiEventBranchName,&fopiEventArray)<<std::endl;
        //build tree Index for fast access later on:
        std::cout<<"FopiDataReaderTask::Init() Building tree index ... "
                <<std::endl;
        TString major = fopiEventBranchName;
        TString minor = major;
        //please let us never change the accessor names
        major += ".GetSpillNum()";
        minor += ".GetEvtNum()";
        int nE = friendChain->BuildIndex(major, minor);
        std::cout<<"   ... Index has "<<nE<<" entries"<<std::endl;
        fFopiIndex = (TTreeIndex*) friendChain->GetTreeIndex();  //Tree has ownership
        //fFopiIndex->Print();
        std::cout<<"chain->SetBranchAddress(tpcEventIdentifierBranchName,&tpcEventIdentifiererArray);"<<std::endl;
        std::cout<<chain->SetBranchAddress(tpcEventIdentifierBranchName,&tpcEventIdentifiererArray)<<std::endl;
    }
    
    
    InitHistos();
    
}


void CdcTrackChain::Exec(int cutOnValidVertex, double momCutLow,double momCutHigh){
    std::cout<<"CdcTrackChain::Exec()"<<std::endl;
    
    int nEntries = chain->GetEntries();
    std::cout<<nEntries<<std::endl;
    if(nEvents==0){
        nEvents=nEntries;
    }else if(nEvents>nEntries){
        nEvents=nEntries;
    }
    std::cout<<"Chain with "<<nEntries<<" events loaded, looping over "<<nEvents<<" events"<<std::endl;
    
    TStopwatch clock;
    
    clock.Start();
    int nFailedEvents=0;
    for(unsigned int i_ev=0;i_ev<nEvents;++i_ev){
        if(cdcArray!=NULL){
            cdcArray->Delete();
        }
        if(fVerbose){
            if(i_ev%5000==0){
                std::cout<<"."<<std::flush;
            }
            if(i_ev%100000==0){
                std::cout<<"------------------- "<<i_ev<<"/"<<nEvents<<" ("<<(i_ev+0.0)/(nEvents+0.0)*100<<" %)"
                        <<"---------------------------------"<<std::endl;
                clock.Print();
                clock.Continue();
            }
        }
        
        chain->GetEvent(i_ev);
        
        //prelim cuts:
        
        //Radius
        double r_test=20;
        
        if(!fGFMode){
            FopiEvent* fopiEv=(FopiEvent*) fopiEventArray->At(0);
            if(fopiEv==NULL){
                continue;
            }
            //fopiEv->GetVertex().Print();
            //http://mathworld.wolfram.com/Circle-CircleIntersection.html
            unsigned int nTracks=cdcArray->GetEntries();
            for(unsigned i_tr=0;i_tr<nTracks;++i_tr){
                CdcTrack* theTrack = (CdcTrack*)cdcArray->UncheckedAt(i_tr);
                double  mom = fabs(theTrack->GetMom());
                double theta=theTrack->GetTheta();
                double pt=mom*TMath::Sin(theta);
                int charge = theTrack->GetCharge();
                //std::cout<<charge<<std::endl;
                TVector3 trackCenter(theTrack->GetMx(),theTrack->GetMy(),0);
                
                double phi0=trackCenter.Phi();
                double d=trackCenter.Mag();
                double r_tr=fabs(theTrack->GetRadius());
                if((r_tr+r_test)<d){
                    continue;
                }
                if((r_test+d)<r_tr){
                    continue;
                }
                
                double r_test2=5;
                if((r_tr+r_test2)<d){
                    continue;
                }
                if((r_test2+d)<r_tr){
                    continue;
                }
                double z0=theTrack->GetZ0();
                if(z0>20 or z0<-20){
                    continue;
                }
                double x=(d*d-r_tr*r_tr+r_test*r_test)/(2*d);
                double y=sqrt(r_test*r_test-x*x);
                
                TVector3 trackPos(x,y,0);
                trackPos.RotateZ(phi0);
                TVector3 cent_to_pos=trackPos-trackCenter;
                //TVector3 cent_to_pos=trackCenter-trackPos;
                
                double phiIncrement=acos(x/r_test);
                
                double phi=phi0+phiIncrement*charge/abs(charge);
                
                cent_to_pos.RotateZ(-TMath::PiOver2()*charge/abs(charge));
                double phiMom=cent_to_pos.Phi();
                
                //phiMom += TMath::PiOver2()*charge/abs(charge);
                phiMom  = phiMom*TMath::RadToDeg();
                phi     = phi*TMath::RadToDeg();
                if(phiMom>180){
                    phiMom-=360;
                }
                if(phiMom<-180){
                    phiMom+=360;
                }
                if(phi>180){
                    phi=phi-360;
                }
                if(phi<-180){
                    phi=phi+360;
                }
               
                double dPhi=TpcRefGFTrkResCalc_Alignment::phiDifference(phiMom*TMath::DegToRad(),phi*TMath::DegToRad());
                
                
                if(phi!=phi){
                    std::cout<<r_tr<<", d:"<<d<<", x: "<<x<<std::endl; 
                    trackCenter.Print();
                }
                std::set<unsigned int> sectorIds;
                std::set<int> sectorIdsLeftRight;
                int nHits=theTrack->GetNpoint();
                for(int iHit=0;iHit<nHits;++iHit){
                    CdcHit* hit= (CdcHit*)cdcHitArray->At(theTrack->getNthHitID(iHit));
                    
                    if(hit!=NULL){
                        TVector3 wirePos=hit->GetWirePosA();
                        double drift=hit->GetHitDrift();
                        unsigned int sectorId=getCdcSector(wirePos);
                        int sectorIdLeftRight=sectorId;
                        if(drift<0){
                            sectorIdLeftRight*=-1;
                            if(sectorIdLeftRight==0){
                                sectorIdLeftRight=-999;
                            }
                        }
                        
                        
                        sectorIds.insert(sectorId);
                        sectorIdsLeftRight.insert(sectorIdLeftRight);
                        //std::cout<<sectorId<<std::endl;
                    }else{
                        std::cout<<"hit not found!!!!"<<std::endl;
                    }
                }
                int nSectors=sectorIds.size();
                int nSectorsLR=sectorIdsLeftRight.size();
                for(int vtxCut =-1;vtxCut<2;vtxCut+=2){
                    //todo replace == 
                    if((fabs(fopiEv->GetVertex().X())<1e-20&&vtxCut>0)){
                        continue;
                    }
                    if((fabs(fopiEv->GetVertex().X())>=1e-20
                        &&vtxCut<0)){
                        continue;
                    }
                    histos["histPhi"]->Fill(phi);
                    histos["histDPhiVPhi"]->Fill(phi,dPhi);
                    if(charge>0){
                        histos["histDPhiVPhiPlus"]->Fill(phi,dPhi);
                        histos["histPhiPlus"]->Fill(phi);
                    }else{
                        histos["histDPhiVPhiMinus"]->Fill(phi,dPhi);
                        histos["histPhiMinus"]->Fill(phi);
                    }
                    {
                        TString postfix = TString::Format("vtxCut_%d",vtxCut);
                        histos["histPhi"+postfix]->Fill(phi);
                        histos["histDPhiVPhi"+postfix]->Fill(phi,dPhi);
                        if(charge>0){
                            histos["histDPhiVPhiPlus"+postfix]->Fill(phi,dPhi);
                            histos["histPhiPlus"+postfix]->Fill(phi);
                        }else{
                            histos["histDPhiVPhiMinus"+postfix]->Fill(phi,dPhi);
                            histos["histPhiMinus"+postfix]->Fill(phi);
                        }
                    }
                    for(double cutMom=0.200;cutMom<1.9;cutMom+=0.200){
                        if(pt>cutMom or pt<(cutMom-0.200)){
                            continue;
                        }
                        TString postfix = TString::Format("vtxCut_%d_mom_%f_%f",vtxCut, cutMom-0.2,cutMom);
                        
                        histos["histPhi"+postfix]->Fill(phi);
                        histos["histDPhiVPhi"+postfix]->Fill(phi,dPhi);
                        histos["nSectors"+postfix]->Fill(nSectors);
                        histos["nSectorsLR"+postfix]->Fill(nSectorsLR);
                        if(charge>0){
                            histos["histDPhiVPhiPlus"+postfix]->Fill(phi,dPhi);
                            histos["histPhiPlus"+postfix]->Fill(phi);
                            histos["nSectorsPlus"+postfix]->Fill(nSectors);
                            histos["nSectorsLRPlus"+postfix]->Fill(nSectorsLR);
                        }else{
                            histos["histDPhiVPhiMinus"+postfix]->Fill(phi,dPhi);
                            histos["histPhiMinus"+postfix]->Fill(phi);
                            histos["nSectorsMinus"+postfix]->Fill(nSectors);
                            histos["nSectorsLRMinus"+postfix]->Fill(nSectorsLR);
                        }
                    }
                    
                    
                }//end loop over vtx cut
                
            }//end loop over tracks       
        }else{
            TpcEventIdentifier* tpcid=(TpcEventIdentifier*) tpcEventIdentifiererArray->At(0);
            unsigned int tpcEvent = tpcid->getEventInSpill();
            unsigned int tpcSpill = tpcid->getSpill();
            
            //fFopiIndex->Print();
            int treeEv = fFopiIndex->GetEntryNumberWithIndex((int)tpcSpill,
                                                             (int)tpcEvent);
            if(treeEv<0) {
                std::cout<<"error looking for entry with"<<"\n"
                        <<"   Spill: "<<tpcSpill<<"\n"
                       <<"   Event: "<<tpcEvent<<"\n"
                      <<"# "<<++nFailedEvents<<"/ "<<i_ev<<" "<<(double) nFailedEvents/(double) i_ev<<"\n";
//                Warning("FopiDataReaderTask::Exec()",
//                        "No matching Entry in Fopi Tree could be found");
                continue;
            }
            
            //            std::cout<<"###########################################################"<<std::endl;
            //            std::cout<<"TpcEventInChain: "<<i_ev<<", spill: "<<tpcSpill<<", ev: "<<tpcEvent<<std::endl;
            friendChain->GetEvent(treeEv);
            FopiEvent* fopiEv=(FopiEvent*) fopiEventArray->At(0);
            //            std::cout<<"FopiEventInChain: "<<treeEv<<", spill: "<<fopiEv->GetSpillNum()
            //                    <<", ev: "<<fopiEv->GetEvtNum()<<std::endl;
            //            std::cout<<"###########################################################"<<std::endl;
            
            if(fopiEv==NULL){
                
                //std::cout<<"no FOPI event!!!"<<treeEv<<" ,"<<friendChain->GetEntries()<<std::endl;
                continue;
            }
            
            
            unsigned int nTracks=cdcArray->GetEntries();
            //            if(cutOnValidVertex==1&&nTracks!=1){
            //                continue;
            //            }if(cutOnValidVertex==-1&&nTracks==1){
            //                continue;
            //            }
            for(unsigned i_tr=0;i_tr<nTracks;++i_tr){
                GFTrack* theTrack = (GFTrack*) cdcArray->At(i_tr);
                GFAbsTrackRep* theRep= theTrack->getCardinalRep();
                TVector3 pos,mom;
                if(theRep->getStatusFlag()!=0){
                    continue;
                }
                               
                GFDetPlane targetPlane(TVector3(0,0,0),TVector3(0,0,1));
                
                TVector3 testPos;
                try{
                    testPos=theRep->getPos(targetPlane);
                    theRep->extrapolateToCylinder(20,pos,mom);
                }catch(GFException exp){
                    continue;
                }
                int charge=theRep->getCharge();
                histos["vertex2d"]->Fill(testPos.X(),testPos.Y());
                if(charge>0){
                    histos["vertex2dPlus"]->Fill(testPos.X(),testPos.Y());
                }else{
                    histos["vertex2dMinus"]->Fill(testPos.X(),testPos.Y());
                }
                if(pos.Z()>50 or pos.Z()<-50 or testPos.Perp()>5){
                    continue;
                }
                
                double phi=pos.Phi();
                double momPhi=mom.Phi();
                double momTheta=mom.Theta()*TMath::RadToDeg();
                if(momTheta<5){
                    continue;
                }
                //takes input in radian returns in degree, arg1-arg2
                TVector3 mom2d(mom.X(),mom.Y(),0);
                TVector3 pos2d(pos.X(),pos.Y(),0);
                double dPhi=TpcRefGFTrkResCalc_Alignment::phiDifference(momPhi,phi);
                //double dPhi=mom2d.Angle(pos2d);
                //dPhi=dPhi*TMath::RadToDeg();
                phi=phi*TMath::RadToDeg();
                
                for(int vtxCut =-1;vtxCut<2;vtxCut+=2){
                    //todo replace == 
                    if((fabs(fopiEv->GetVertex().X())<1e-20&&vtxCut>0)){
                        continue;
                    }
                    if((fabs(fopiEv->GetVertex().X())>=1e-20
                        &&vtxCut<0)){
                        continue;
                    }
                    histos["histPhi"]->Fill(phi);
                    histos["histDPhiVPhi"]->Fill(phi,dPhi);
                    if(charge>0){
                        histos["histDPhiVPhiPlus"]->Fill(phi,dPhi);
                        histos["histPhiPlus"]->Fill(phi);
                    }else{
                        histos["histDPhiVPhiMinus"]->Fill(phi,dPhi);
                        histos["histPhiMinus"]->Fill(phi);
                    }
                    {
                        TString postfix = TString::Format("vtxCut_%d",vtxCut);
                        histos["histPhi"+postfix]->Fill(phi);
                        histos["histDPhiVPhi"+postfix]->Fill(phi,dPhi);
                        if(charge>0){
                            histos["histDPhiVPhiPlus"+postfix]->Fill(phi,dPhi);
                            histos["histPhiPlus"+postfix]->Fill(phi);
                        }else{
                            histos["histDPhiVPhiMinus"+postfix]->Fill(phi,dPhi);
                            histos["histPhiMinus"+postfix]->Fill(phi);
                        }
                    }
                    for(double cutMom=0.200;cutMom<1.9;cutMom+=0.200){
                        if(mom.Perp()>cutMom or mom.Perp()<(cutMom-0.200)){
                            continue;
                        }
                        TString postfix = TString::Format("vtxCut_%d_mom_%f_%f",vtxCut, cutMom-0.2,cutMom);
                        
                        histos["histPhi"+postfix]->Fill(phi);
                        histos["histDPhiVPhi"+postfix]->Fill(phi,dPhi);
                        if(charge>0){
                            histos["histDPhiVPhiPlus"+postfix]->Fill(phi,dPhi);
                            histos["histPhiPlus"+postfix]->Fill(phi);
                        }else{
                            histos["histDPhiVPhiMinus"+postfix]->Fill(phi,dPhi);
                            histos["histPhiMinus"+postfix]->Fill(phi);
                        }
                    }//end loop over cutMom
                }//end loop over cutVtx
            }//end if genfit
        }
        
    }//end loop over events
    
}

void CdcTrackChain::InitHistos(){
    TH2D* vertex2d = new TH2D("vertex2d","vertex2d",200,-20,20,200,-20,20);
    histos["vertex2d"]=vertex2d;
    TH2D* vertex2dPlus = new TH2D("vertex2dPlus","vertex2dPlus",200,-20,20,200,-20,20);
    histos["vertex2dPlus"]=vertex2dPlus;
    TH2D* vertex2dMinus = new TH2D("vertex2dMinus","vertex2dMinus",200,-20,20,200,-20,20);
    histos["vertex2dMinus"]=vertex2dMinus;
    for(int vtxCut =-1;vtxCut<2;vtxCut+=2){
        for(double cutMom=0.200;cutMom<1.9;cutMom+=0.200){
            TString postfix = TString::Format("vtxCut_%d_mom_%f_%f",vtxCut, cutMom-0.2,cutMom);
            TString prettyPostfix = TString::Format(" vtxCut%d %f< mom <%f",vtxCut,cutMom-0.2,cutMom);
            TH1D* histPhi=new TH1D("phiHist"+postfix,"#phi"+prettyPostfix,360,-180,180);
            histos["histPhi"+postfix]=histPhi;
            TH1D* histPhiPlus=new TH1D("phiHistPlus"+postfix,"#phi plus"+prettyPostfix,360,-180,180);
            histos["histPhiPlus"+postfix]=histPhiPlus;
            TH1D* histPhiMinus=new TH1D("phiHistMins"+postfix,"#phi mins"+prettyPostfix,360,-180,180);
            histos["histPhiMinus"+postfix]=histPhiMinus;
            
            TH2D* histDPhiVPhi=new TH2D("dphiVphiHist"+postfix,"d#phi vs #phi"+prettyPostfix,360,-180,180,360,-50,50);
            histos["histDPhiVPhi"+postfix]=histDPhiVPhi;
            TH2D* histDPhiVPhiPlus=new TH2D("dphiVphiHistPlus"+postfix,"d#phi vs #phi plus"+prettyPostfix,360,-180,180,360,-50,50);
            histos["histDPhiVPhiPlus"+postfix]=histDPhiVPhiPlus;
            TH2D* histDPhiVPhiMinus=new TH2D("dphiVphiHistMins"+postfix,"d#phi vs #phi mins"+prettyPostfix,360,-180,180,360,-50,50);
            histos["histDPhiVPhiMinus"+postfix]=histDPhiVPhiMinus;
            
            
            
            TH1D* histNSectorsTmp=new TH1D("+nSectors"+postfix,"nSectors"+prettyPostfix,16,0,16);
            TH1D* histNSectorsLRTmp=new TH1D("nSectorsLR"+postfix,"nSectors Left/right"+prettyPostfix,16,0,16);
            TH1D* histNSectorsTmpPlus=new TH1D("nSectorsPlus"+postfix,"nSectors Plus"+prettyPostfix,16,0,16);
            TH1D* histNSectorsLRTmpPlus=new TH1D("nSectorsLRPlus"+postfix,"nSectors Left/right Plus"+prettyPostfix,16,0,16);
            TH1D* histNSectorsTmpMinus=new TH1D("nSectorsMinus"+postfix,"nSectors Minus"+prettyPostfix,16,0,16);
            TH1D* histNSectorsLRTmpMinus=new TH1D("nSectorsLRMinus"+postfix,"nSectors Left/right Minus"+prettyPostfix,16,0,16);
            histos["nSectors"+postfix]=histNSectorsTmp;
            histos["nSectorsLR"+postfix]=histNSectorsLRTmp;
            histos["nSectorsPlus"+postfix]=histNSectorsTmpPlus;
            histos["nSectorsLRPlus"+postfix]=histNSectorsLRTmpPlus;
            histos["nSectorsMinus"+postfix]=histNSectorsTmpMinus;
            histos["nSectorsLRMinus"+postfix]=histNSectorsLRTmpMinus;
            
        }
    }
    TH1D* histPhi=new TH1D("phiHist","#phi",360,-180,180);
    histos["histPhi"]=histPhi;
    TH1D* histPhiPlus=new TH1D("phiHistPlus","#phi plus",360,-180,180);
    histos["histPhiPlus"]=histPhiPlus;
    TH1D* histPhiMinus=new TH1D("phiHistMins","#phi mins",360,-180,180);
    histos["histPhiMinus"]=histPhiMinus;
    {
        TH2D* histDPhiVPhi=new TH2D("dphiVphiHist","d#phi vs #phi",360,-180,180,360,-50,50);
        histos["histDPhiVPhi"]=histDPhiVPhi;
        TH2D* histDPhiVPhiPlus=new TH2D("dphiVphiHistPlus","d#phi vs #phi plus",360,-180,180,360,-50,50);
        histos["histDPhiVPhiPlus"]=histDPhiVPhiPlus;
        TH2D* histDPhiVPhiMinus=new TH2D("dphiVphiHistMins","d#phi vs #phi mins",360,-180,180,360,-50,50);
        histos["histDPhiVPhiMinus"]=histDPhiVPhiMinus;
    }
    {
        TString postfix = TString::Format("vtxCut_%d",-1);
        TString prettyPostfix = TString::Format(" vtxCut%d ",-1);
        TH1D* histPhiTmp=new TH1D("phiHist"+postfix,"#phi"+prettyPostfix,360,-180,180);
        histos["histPhi"+postfix]=histPhiTmp;
        TH1D* histPhiPlusTmp=new TH1D("phiHistPlus"+postfix,"#phi plus"+prettyPostfix,360,-180,180);
        histos["histPhiPlus"+postfix]=histPhiPlusTmp;
        TH1D* histPhiMinusTmp=new TH1D("phiHistMins"+postfix,"#phi mins"+prettyPostfix,360,-180,180);
        histos["histPhiMinus"+postfix]=histPhiMinusTmp;
        
        TH2D* histDPhiVPhi=new TH2D("dphiVphiHist"+postfix,"d#phi vs #phi"+prettyPostfix,360,-180,180,360,-50,50);
        histos["histDPhiVPhi"+postfix]=histDPhiVPhi;
        TH2D* histDPhiVPhiPlus=new TH2D("dphiVphiHistPlus"+postfix,"d#phi vs #phi plus"+prettyPostfix,360,-180,180,360,-50,50);
        histos["histDPhiVPhiPlus"+postfix]=histDPhiVPhiPlus;
        TH2D* histDPhiVPhiMinus=new TH2D("dphiVphiHistMins"+postfix,"d#phi vs #phi mins"+prettyPostfix,360,-180,180,360,-50,50);
        histos["histDPhiVPhiMinus"+postfix]=histDPhiVPhiMinus;
    }
    {
        TString postfix = TString::Format("vtxCut_%d",1);
        TString prettyPostfix = TString::Format(" vtxCut%d ",1);
        TH1D* histPhiTmp=new TH1D("phiHist"+postfix,"#phi"+prettyPostfix,360,-180,180);
        histos["histPhi"+postfix]=histPhiTmp;
        TH1D* histPhiPlusTmp=new TH1D("phiHistPlus"+postfix,"#phi plus"+prettyPostfix,360,-180,180);
        histos["histPhiPlus"+postfix]=histPhiPlusTmp;
        TH1D* histPhiMinusTmp=new TH1D("phiHistMins"+postfix,"#phi mins"+prettyPostfix,360,-180,180);
        histos["histPhiMinus"+postfix]=histPhiMinusTmp;
        
        TH2D* histDPhiVPhi=new TH2D("dphiVphiHist"+postfix,"d#phi vs #phi"+prettyPostfix,360,-180,180,360,-50,50);
        histos["histDPhiVPhi"+postfix]=histDPhiVPhi;
        TH2D* histDPhiVPhiPlus=new TH2D("dphiVphiHistPlus"+postfix,"d#phi vs #phi plus"+prettyPostfix,360,-180,180,360,-50,50);
        histos["histDPhiVPhiPlus"+postfix]=histDPhiVPhiPlus;
        TH2D* histDPhiVPhiMinus=new TH2D("dphiVphiHistMins"+postfix,"d#phi vs #phi mins"+prettyPostfix,360,-180,180,360,-50,50);
        histos["histDPhiVPhiMinus"+postfix]=histDPhiVPhiMinus;
    }
    

}
void CdcTrackChain::writeHistos(){
    TFile* outFile= new TFile(outFileName,"RECREATE");
    outFile->cd();
    outFile->mkdir("plus");
    outFile->mkdir("minus");
    std::map<TString,TH1*>::iterator it;
    for(it=histos.begin();it!=histos.end();++it){
        outFile->cd();
        if((*it).first.Contains("Minus")){
            outFile->cd("minus");
            (*it).second->Write();
            delete (*it).second;
        }else if ((*it).first.Contains("Plus")){
            outFile->cd("plus");
            (*it).second->Write();
            delete (*it).second;
        }else{
            (*it).second->Write();
            delete (*it).second;
        }
        
    }
    outFile->Close();
    
}
int CdcTrackChain::getCdcSector(TVector3 wirePos){
    double phi = wirePos.Phi()*TMath::RadToDeg();
    //std::cout<<"Phi: "<<phi;
    for(unsigned int i=0;i<16;++i){
        if(phi<(-155+i*22.5)){
          //  std::cout<<", sectorId: "<<i<<std::endl;
            return i;
        }
    }
    return -1;
}
int CdcTrackChain::getCdcSector(double phi){
    //std::cout<<"Phi: "<<phi;
    for(unsigned int i=0;i<16;++i){
        if(phi<(-155+i*22.5)){
          //  std::cout<<", sectorId: "<<i<<std::endl;
            return i;
        }
    }
    return -1;
}    
    


//                    if((fabs(fopiEv->GetVertex().X())<1e-20&&vtxCut>0)){
//                        continue;
//                    }
//                    if((fabs(fopiEv->GetVertex().X())>=1e-20
//                        &&vtxCut<0)){
//                        continue;
//                    }

//                    histos["histPhi"]->Fill(phi);
//                    if(charge>0){
//                        histos["histPhiPlus"]->Fill(phi);
//                    }else{
//                        histos["histPhiMinus"]->Fill(phi);
//                    }
//                    {
//                        TString postfix = TString::Format("vtxCut_%d",vtxCut);
//                        histos["histPhi"+postfix]->Fill(phi);
//                        if(charge>0){
//                            histos["histPhiPlus"+postfix]->Fill(phi);
//                        }else{
//                            histos["histPhiMinus"+postfix]->Fill(phi);
//                        }
//                    }
//                    for(double cutMom=0.200;cutMom<1.9;cutMom+=0.200){
//                        if(mom>cutMom or mom<(cutMom-0.200)){
//                            continue;
//                        }
//                        TString postfix = TString::Format("vtxCut_%d_mom_%f_%f",vtxCut, cutMom-0.2,cutMom);
//                        histos["histPhi"+postfix]->Fill(phi);
//                        if(charge>0){
//                            histos["histPhiPlus"+postfix]->Fill(phi);
//                        }else{
//                            histos["histPhiMinus"+postfix]->Fill(phi);
//                        }
//                    }


ClassImp(CdcTrackChain)