[u/mrichter/AliRoot.git] / PWGCF / Correlations / JCORRAN / AliJEfficiencyScanner.cxx

//
//#include <Riostream.h>
// #include <TChain.h>
// #include <TVectorT.h> 
// #include <TVector3.h> 
// #include <TFile.h>
// #include <TH1.h> 
// #include <TObjArray.h>
// #include <TObjString.h>
// #include <TFormula.h>
// #include <TString.h>
// #include <TRefArray.h>
// #include <TNtuple.h>
// #include <TArrayF.h>

#include <TClonesArray.h>
#include <TH1D.h>
#include <TH2D.h>

#include "AliJEfficiencyScanner.h" 
#include "AliJTrack.h"
#include "AliJMCTrack.h"
#include "AliJPhoton.h"
#include "AliJEventHeader.h"
#include "AliJRunHeader.h"
#include <TFile.h>
#include <TRandom.h>


ClassImp(AliJEfficiencyScanner);

//______________________________________________________________________________
AliJEfficiencyScanner::AliJEfficiencyScanner() :   
    TNamed(),
	fMBTriggMask(0),
    fTrackList(0),
    fMCTrackList(0x0),
    fEventHeader(0x0),
    fRunHeader(0x0),
    fhVertexZMC(0x0),
    fhVertexZTrigg(0x0),
    fhVertexZTriggVtx(0x0),
    fhVZRawMC(0x0),
    fhVZRecMC(0x0),
    fhVZRecAccMC(0x0),
    fh2VtxCent(0x0),
    fhL0Input(0x0),
    fhTriggerAlice(0x0),
    fhZVtxMCAll(0x0),
    fhZVtxMCTrigg(0x0),
    fhZVtxMCTriggVtx(0x0),
    fhZVtxRecAll(0x0),
    fhZVtxRecTrigg(0x0),
    fhZVtxRecTriggVtx(0x0),
    fVtxReFunc(0x0),
    fVtxMCFunc(0x0),
    fVtxRatioFunc(0x0),
    fVtxRatioMax(1)
{
    //Default constructor
    for( int ivtx=0;ivtx<kNVtxBin;ivtx++ ){
        for( int icent=0;icent<kNCentBin;icent++ ){
            fhChargedPtMCTriggVtx[ivtx][icent] = 0x0;
            for( int ifilter=0;ifilter<AliJTrackCut::kJNTrackCuts;ifilter++ ){
                for( int ipri=0;ipri<kNPrimaryStatus;ipri++ ){
                    for( int itt=0;itt<kNTrackType;itt++ ){
                        fhChargedPtMCRecoCentVtx[ivtx][icent][ifilter][ipri][itt] = 0x0;
                    }
                }
            }
        }
    }

}

//______________________________________________________________________________
AliJEfficiencyScanner::AliJEfficiencyScanner(const char *name):
    TNamed(name,name), 
	fMBTriggMask(0),
    fTrackList(0),
    fMCTrackList(0x0),
    fEventHeader(0x0),
    fRunHeader(0x0),
    fhVertexZMC(0x0),
    fhVertexZTrigg(0x0),
    fhVertexZTriggVtx(0x0),
    fhVZRawMC(0x0),
    fhVZRecMC(0x0),
    fhVZRecAccMC(0x0),
    fh2VtxCent(0x0),
    fhL0Input(0x0),
    fhTriggerAlice(0x0),
    fhZVtxMCAll(0x0),
    fhZVtxMCTrigg(0x0),
    fhZVtxMCTriggVtx(0x0),
    fhZVtxRecAll(0x0),
    fhZVtxRecTrigg(0x0),
    fhZVtxRecTriggVtx(0x0),
    fVtxReFunc(0x0),
    fVtxMCFunc(0x0),
    fVtxRatioFunc(0x0),
    fVtxRatioMax(1)
{
    // Constructor
    for( int ivtx=0;ivtx<kNVtxBin;ivtx++ ){
        for( int icent=0;icent<kNCentBin;icent++ ){
            fhChargedPtMCTriggVtx[ivtx][icent] = 0x0;
            for( int ifilter=0;ifilter<AliJTrackCut::kJNTrackCuts;ifilter++ ){
                for( int ipri=0;ipri<kNPrimaryStatus;ipri++ ){
                    for( int itt=0;itt<kNTrackType;itt++ ){
                        fhChargedPtMCRecoCentVtx[ivtx][icent][ifilter][ipri][itt] = 0x0;
                    }
                }
            }
        }
    }
}

//____________________________________________________________________________
AliJEfficiencyScanner::AliJEfficiencyScanner(const AliJEfficiencyScanner& ap) :
    TNamed(ap.GetName(), ap.GetTitle()),
	fMBTriggMask(ap.fMBTriggMask),
    fTrackList(ap.fTrackList),
    fMCTrackList(ap.fMCTrackList),
    fEventHeader( ap.fEventHeader ),
    fRunHeader(ap.fRunHeader),
    fhVertexZMC(ap.fhVertexZMC),
    fhVertexZTrigg(ap.fhVertexZTrigg),
    fhVertexZTriggVtx(ap.fhVertexZTriggVtx),
    fhVZRawMC(ap.fhVZRawMC),
    fhVZRecMC(ap.fhVZRecMC),
    fhVZRecAccMC(ap.fhVZRecAccMC),
    fh2VtxCent(ap.fh2VtxCent),
    fhL0Input(ap.fhL0Input),
    fhTriggerAlice(ap.fhTriggerAlice),
    fhZVtxMCAll(ap.fhZVtxMCAll),
    fhZVtxMCTrigg(ap.fhZVtxMCTrigg),
    fhZVtxMCTriggVtx(ap.fhZVtxMCTriggVtx),
    fhZVtxRecAll(ap.fhZVtxRecAll),
    fhZVtxRecTrigg(ap.fhZVtxRecTrigg),
    fhZVtxRecTriggVtx(ap.fhZVtxRecTriggVtx),
    fVtxReFunc(ap.fVtxReFunc),
    fVtxMCFunc(ap.fVtxMCFunc),
    fVtxRatioFunc(ap.fVtxRatioFunc),
    fVtxRatioMax(ap.fVtxRatioMax)
{ 
    // cpy ctor
    for( int ivtx=0;ivtx<kNVtxBin;ivtx++ ){
        for( int icent=0;icent<kNCentBin;icent++ ){
            fhChargedPtMCTriggVtx[ivtx][icent] = ap.fhChargedPtMCTriggVtx[ivtx][icent];
            for( int ifilter=0;ifilter<AliJTrackCut::kJNTrackCuts;ifilter++ ){
                for( int ipri=0;ipri<kNPrimaryStatus;ipri++ ){
                    for( int itt=0;itt<kNTrackType;itt++ ){
                        fhChargedPtMCRecoCentVtx[ivtx][icent][ifilter][ipri][itt] = ap.fhChargedPtMCRecoCentVtx[ivtx][icent][ifilter][ipri][itt];
                    }
                }
            }
        }
    }

}

//_____________________________________________________________________________
AliJEfficiencyScanner& AliJEfficiencyScanner::operator = (const AliJEfficiencyScanner& ap)
{
    // assignment operator

    this->~AliJEfficiencyScanner();
    new(this) AliJEfficiencyScanner(ap);
    return *this;
}

//______________________________________________________________________________
AliJEfficiencyScanner::~AliJEfficiencyScanner()
{
    // destructor 
}

//________________________________________________________________________

void AliJEfficiencyScanner::UserCreateOutputObjects()
{  
    //=== create the jcorran outputs objects
    cout<<"DEBUG Start AliJEfficiencyScanner::UserCreateOutputObjects() "<<"\t"<<gDirectory<<endl;

    double ptbin[300] = {0};
    double pt = 0;
    int i = 0;
    for( i=0;i<300;i++ ){
        ptbin[i] = pt;
        if( pt > 50 ) break;
        if( pt < 3 ) pt+= 0.05;
        else if( pt < 5 ) pt+= 0.1;
        else if( pt < 10 ) pt+= 1;
        else pt+= 1;
    }
    cout<<"n Ptbin = "<<i<<endl;
    int nPtBin = i-1;

    fh2VtxCent = new TH2D("h2VtxCent","h2VtxCent",100,0,10,110,0,110 );
    fh2VtxCent->SetDirectory(gDirectory);

    int nVtxBin = kNVtxBin;
    int nCentBin = kNCentBin;
    if( fRunHeader->IsPP() ) nCentBin = 1;

    TString name ="";
    for( int ivtx=0;ivtx<nVtxBin;ivtx++ ){
        for( int icent=0;icent<nCentBin;icent++ ){
            fhChargedPtMC[ivtx][icent] = AddTH1D( Form("hChargedPtMC%02d%02d",ivtx,icent), 
                    new TH1D("","", nPtBin, ptbin) );
            fhChargedPtMCTriggVtx[ivtx][icent]=AddTH1D( 
                    Form("hChargedPtMCTriggVtx%02d%02d", ivtx, icent ),
                    new TH1D("", "", nPtBin, ptbin) );
            fhChargedPtMCTrigg[ivtx][icent]=AddTH1D( 
                    Form("hChargedPtMCTrigg%02d%02d", ivtx, icent ),
                    new TH1D("", "", nPtBin, ptbin) );
            //for any MC track filled with MC pt in triggered event
            name = Form("h2DChargedPtTrigg%02d%02d",ivtx,icent);
            fh2DChargedPtTrigg[ivtx][icent]=AddTH2D(name, new TH2D(name,name, nPtBin, ptbin, 20, -0.8, 0.8));
            //for any MC track filled with MC pt in triggered event with rec vertex
            name = Form("h2DChargedPtTriggVtx%02d%02d",ivtx,icent);
            fh2DChargedPtTriggVtx[ivtx][icent]=AddTH2D(name, new TH2D(name,name, nPtBin, ptbin, 20, -0.8, 0.8));

        }
    }
    fhTriggerAlice = AddTH1D("hTriggerAlice",new TH1D("","",32,0,32));
    fhL0Input = AddTH1D("hL0Input",new TH1D("","",32,0,32));
    // DCA bin
    double dcaBin[1000];
    double dbin=-50-1;
    int ndbin = 0;
    dcaBin[0] = dbin;
    double tol = 1e-5;
    while(dbin < 50-1){
        if( fabs(dbin) < 2-tol ) dbin+=0.01;
        else if( fabs(dbin) < 5-tol  ) dbin+= 0.05;
        else if( fabs(dbin) <= 10-tol ) dbin+= 0.1;
        else dbin += 1;
        if( fabs(dbin) < tol ) dbin=0;
        dcaBin[ndbin++] = dbin;
    }

    for( int ivtx=0;ivtx<nVtxBin;ivtx++ ){
        for( int icent=0;icent<nCentBin;icent++ ){
            for( int ifilter=0;ifilter<AliJTrackCut::kJNTrackCuts;ifilter++ ){

                name = Form("h2DChargedPtAll%02d%02d%02d",ivtx,icent,ifilter);
                fh2DChargedPtAll[ivtx][icent][ifilter]=AddTH2D(name, new TH2D(name,name, nPtBin, ptbin, 20, -0.8, 0.8));
                name = Form("h2DChargedPtRec%02d%02d%02d",ivtx,icent,ifilter);
                fh2DChargedPtRec[ivtx][icent][ifilter]=AddTH2D(name, new TH2D(name,name, nPtBin, ptbin, 20, -0.8, 0.8));
                for( int ipri=0;ipri<kNPrimaryStatus;ipri++ ){
                    for( int itt=0;itt<kNTrackType;itt++ ){
                        name = Form("hChargedPtMCRecoCentVtx%02d%02d%02d%02d%02d", ivtx, icent, ifilter, ipri, itt );
                        fhChargedPtMCRecoCentVtx[ivtx][icent][ifilter][ipri][itt]=AddTH1D(name, new TH1D( name, name, nPtBin, ptbin));
                    }
                }
                name = Form("hDCA2VertexXY%02d%02d%02d", ivtx, icent, ifilter );
                fhDCA2VertexXY[ivtx][icent][ifilter] = AddTH1D(name, new TH1D("","",ndbin-1, dcaBin ));
                name = Form("hDCA2VertexZ%02d%02d%02d", ivtx, icent, ifilter );
                fhDCA2VertexZ[ivtx][icent][ifilter] = AddTH1D( name, new TH1D("","",ndbin-1, dcaBin ));

            }
        }
    }
//  TH2D * fh2MultGenRawPrimary[AliJTrackCut::kJNTrackCuts];
//  TH2D * fh2MultGenRawAll[AliJTrackCut::kJNTrackCuts];

//  TH1D * fhDCA2VertexXY[kNVtxBin][kNCentBin][AliJTrackCut::kJNTrackCuts];
    for( int ifilter=0;ifilter<AliJTrackCut::kJNTrackCuts;ifilter++ ){
        name = Form("h2MultGenRawPrimary%02d", ifilter);
        fh2MultGenRawPrimary[ifilter] = AddTH2D( name, new TH2D(name, name, 100,0,100,100,0,100 ));
        name = Form("h2MultGenRawAll%02d", ifilter);
        fh2MultGenRawAll[ifilter] = AddTH2D( name, new TH2D(name, name, 100,0,100,100,0,100 ));
    }
    double   binsVertexMult[] = {0,1,2,3,4,5,10000};
    int   NbinsVertexMult  = sizeof(binsVertexMult)/sizeof(double)-1;
    double binsVertexZ[]    = {-10,-6,-3,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,3,6,10};
    int   NbinsVertexZ   = sizeof(binsVertexZ)/sizeof(double)-1;
    fhVertexZMC = AddTH2D("hVertexZMC",new TH2D("hVertexZMC","hVertexZMC", NbinsVertexMult, binsVertexMult, NbinsVertexZ, binsVertexZ));
    fhVertexZTrigg = AddTH2D("hVertexZTrigg",new TH2D("hVertexZTrigg","hVertexZTrigg",  NbinsVertexMult, binsVertexMult, NbinsVertexZ, binsVertexZ));
    fhVertexZTriggVtx = AddTH2D("hVertexZTriggVtx",new TH2D("hVertexZTriggVtx","hVertexZTriggVtx",  NbinsVertexMult, binsVertexMult, NbinsVertexZ, binsVertexZ));
    fhVZRawMC = AddTH1D("hVtxZMC", new TH1D("hVtxZMC","VertexZ in MC ",200,-50,50));
    fhVZRecMC = AddTH1D("hVtxZRec", new TH1D("hVtxZRec","VertexZ Rec in MC",200,-50,50));
    fhVZRecAccMC = AddTH1D("hVtxZRecAcc", new TH1D("hVtxZRecAcc","VertexZ Rec in MC and acc",200,-50,50));

	int NvtxBins = 400;
	double vtxedge = 50;
    fhZVtxMCAll = AddTH1D("hZVtxMCAll", new TH1D("hZVtxMCAll","VetrexZ in MC for all event",NvtxBins, -vtxedge, vtxedge ));
    fhZVtxMCTrigg = AddTH1D("hZVtxMCTrigg", new TH1D("hZVtxMCTrigg","VetrexZ in MC for triggered event",NvtxBins, -vtxedge, vtxedge ));
    fhZVtxMCTriggVtx = AddTH1D("hZVtxMCTriggVtx", new TH1D("hZVtxMCTriggVtx","VetrexZ in MC for triggered and vtx reconstructed event",NvtxBins, -vtxedge, vtxedge ));

    fhZVtxRecAll = AddTH1D("hZVtxRecAll", new TH1D("hZVtxRecAll","VetrexZ in Rec for all event",NvtxBins, -vtxedge, vtxedge ));
    fhZVtxRecTrigg = AddTH1D("hZVtxRecTrigg", new TH1D("hZVtxRecTrigg","VetrexZ in Rec for triggered event",NvtxBins, -vtxedge, vtxedge ));
    fhZVtxRecTriggVtx = AddTH1D("hZVtxRecTriggVtx", new TH1D("hZVtxRecTriggVtx","VetrexZ in Rec for triggered and vtx reconstructed event",NvtxBins, -vtxedge, vtxedge ));

    double v0 = -20;
    double v1 = 20;
    fVtxReFunc = new TF1("VtxReFunc", "gaus",v0,v1);
    fVtxReFunc->SetParameters(1 ,  -7.14076e-01, 6.27110 );
    fVtxMCFunc = new TF1("VtxMCFunc", "gaus",v0,v1);
    fVtxMCFunc->SetParameters(1, -4.53674e-01,  5.27088e+00 );
    fVtxRatioFunc = new TF1("VtxRatioFunc", "VtxReFunc/VtxMCFunc",v0,v1);
    fVtxRatioMax = fVtxRatioFunc->GetMaximum();


    cout << "Add(fAliJRunHeader) in UserCreateObject() ======= " << endl;
    cout<<"DEBUG END AliJEfficiencyScanner::UserCreateOutputObjects() "<<endl;

}

//______________________________________________________________________________

void AliJEfficiencyScanner::UserExec(Option_t *option)
{

    // Main loop
    // Called for each event
    JUNUSED(option);

    double zVtxCut = 10;
    double etaCut  = 0.8;
    double nCentBin = 20;

    AliJEventHeader * eventHeader = GetJEventHeader();

    //== TRIGGER( MB ) TODO is this?
    Bool_t triggeredEventMB  = eventHeader->GetTriggerMaskJCorran() & fMBTriggMask; // 1:kMB 8:kINT7 //TODO
    UInt_t trigAlice = eventHeader->GetTriggerMaskAlice();
/*
cout<<"Trigger"<<endl;
	for( int i=0;i<32 ;i++ ) cout<<Form("%02d",i)<<" ";
	cout<<endl;
	for( int i=0;i<32 ;i++ ) cout<<(eventHeader->GetTriggerMaskJCorran()&(1<<i)?1:0)<<"  ";
	cout<<endl;
	for( int i=0;i<32 ;i++ ) cout<<(eventHeader->GetTriggerMaskAlice()&(1<<i)?1:0)<<"  ";
	cout<<endl;
*/
    /*  
    bool isMB1   = trigAlice & (1UL<<0);
    bool isMBBG1 = trigAlice & (1UL<<3);
    bool isGFO = trigAlice & (1UL<<14);
    bool isMBBG3 = trigAlice & (1UL<<1); 
    bool isV0A   = trigAlice & (1UL<<10);
    bool isV0C   = trigAlice & (1UL<11);
    */

    // L0Input
    UInt_t l0Input = eventHeader->GetL0TriggerInputs();
    /* 
    bool isL0V0A  = l0Input & (1UL<<9);
    bool isL0V0C  = l0Input & (1UL<<10);
    bool isL0V0AND  = l0Input & (1UL<<4);
    bool isL0V0OR   = l0Input & (1UL<<1);
    bool isL0V0BeamGas   = l0Input & (1UL<<2);
    bool isL0SMB   = l0Input & (1UL<<3);
    bool isL0T0X   = l0Input & (1UL<<0);
    */


    for( UInt_t i=0;i<32;i++ ){
        if( l0Input  & 1UL<<i ) fhL0Input->Fill(i);
        if( trigAlice & 1UL<<i ) fhTriggerAlice->Fill(i);
/*
        if( isL0V0A != isV0A ) cout<<"J_WARN 0 : L0:V0A != Class:V0A "<<isL0V0A<<"\t"<<isV0A<<endl;
        if( isL0V0C != isV0C ) cout<<"J_WARN 1 : L0:V0C != Class:V0C "<<isL0V0C<<"\t"<<isV0C<<endl;
        if( isL0SMB != isGFO ) cout<<"J_WARN 2 : L0:SMB != Class:GFO "<<isL0SMB<<"\t"<<isGFO<<endl;
        if( isL0V0OR != isMBBG1 ) cout<<"J_WARN 3 : L0:V0OR != Class:MBBG1 "<<isL0V0OR<<"\t"<<isMBBG1<<endl;
        if( (isL0V0A||isL0V0C) != isMBBG1 ) cout<<"J_WARN 4 : L0:V0A|V0C != Class:MBBG1 "<<(isL0V0A||isL0V0C)<<"\t"<<isMBBG1<<endl;
        if( ( isL0T0X || isL0V0OR || isL0SMB && !isL0V0BeamGas ) != isMB1 ) cout<<"J_WARN 5 : L0:MB1 != Class:MB1 "<<isMB1<<endl;
        if( (isL0V0OR || isL0SMB && !isL0V0BeamGas ) != (( isMBBG1 || isGFO ) && isMB1 )) cout<<"J_WARN 6 : L0:MB != Class:MB "<<(isL0V0OR || isL0SMB && !isL0V0BeamGas)<<"\t"<<((isMBBG1 || isGFO ) && isMB1)<<endl; 
        if( (isL0V0OR || isL0SMB ) != (( isMBBG1 || isGFO ) )) cout<<"J_WARN 6 : L0:MB != Class:MB "<<(isL0V0OR || isL0SMB && !isL0V0BeamGas)<<"\t"<<((isMBBG1 || isGFO ) && isMB1)<<endl; 
*/

    }
    //* FOR TEST */
    //Bool_t triggeredEventMB  = ( isMBBG1 || isGFO ) && isMB1;
    //Bool_t triggeredEventMB  = isMB1;
    //Bool_t triggeredEventMB  = ( isMBBG1 || isGFO );
    //
    //== REAL MB with L0 INput
    //triggeredEventMB  = ( (isL0V0OR || isL0SMB) );
    //== CENTRALITY
    double centrality = eventHeader->GetCentrality();
    int iCent = int( centrality/(100./nCentBin));
    if( iCent< 0 || iCent>20 ) return;
    //if( centrality < 0 || centrality > 100 ) return;
    if( fRunHeader->IsPP() ) {iCent = 0; centrality=0;}

    //== Vertex
    //==== Reco
    int ncontributors = eventHeader->GetVtxMult();
    double zVtxRec = eventHeader->GetZVertex();
    Bool_t goodRecVertex = (ncontributors>0) && (fabs(zVtxRec)<=zVtxCut);
    int iVtx = 0;
    //==== MC
    double zVtxMC = eventHeader->GetZVertexMC();
    //==== vtx sampling
    /*
    if( zVtxMC > fVtxRatioFunc->GetXmin() && zVtxMC < fVtxRatioFunc->GetXmax() ){
        double vtxRatio = fVtxRatioFunc->Eval(zVtxMC)/fVtxRatioMax;
        double random = gRandom->Uniform(0,1);
        if( random > vtxRatio ) return;
    }
    */

    fhVZRawMC->Fill(zVtxMC);
    fhVertexZMC->Fill(ncontributors,zVtxMC);
    // for crosscheck MC input events
    fhZVtxMCAll->Fill(zVtxMC);
    fhZVtxRecAll->Fill(zVtxRec);

    if(triggeredEventMB){
        fhVertexZTrigg->Fill(ncontributors,zVtxMC);
       // for crosscheck passing trigger condition
        fhZVtxMCTrigg->Fill(zVtxMC);
        fhZVtxRecTrigg->Fill(zVtxRec);

        if(ncontributors>0){//reconstructed vertex
            fhVZRecMC->Fill(zVtxRec);
       		// for crosscheck passing trigger condition + reconstruted vtx
            fhZVtxMCTriggVtx->Fill(zVtxMC);
            fhZVtxRecTriggVtx->Fill(zVtxRec);
        }
        if(goodRecVertex){
            fhVertexZTriggVtx->Fill(ncontributors,zVtxMC);
            fhVZRecAccMC->Fill(zVtxRec); 
        }
    }

    fh2VtxCent->Fill( zVtxRec, centrality );

    // TODO ? MC or REC?

    int nRawMultPri = 0;
    int nGenMultPri[AliJTrackCut::kJNTrackCuts] = {0};
    int nGenMultAll[AliJTrackCut::kJNTrackCuts] = {0};

    //==============================================
    //  LOOP over MC
    //==============================================

    //if( fabs(zVtxMC) > zVtxCut ) return;
    int nMCTrks = 0;
    if( fRunHeader->IsMC() )
        nMCTrks = GetJMCTracks()->GetEntriesFast();
    for( int it=0; it< nMCTrks; it++ ){
        AliJMCTrack * track = GetJMCTrack( it );//Always IsPhysicalPrimary
        if( !track ) continue;
        if( !track->IsTrue( AliJMCTrack::kPrimary ) ) continue;
        // TODO need? if( ! track->IsFinal() ) continue;
        double eta = track->Eta();
        double pt = track->Pt();
        if( fabs(eta) > etaCut ) continue;
        if( ! track->IsCharged() ) continue;


        nRawMultPri++;
        fhChargedPtMC[iVtx][iCent]->Fill(pt);     //ALL charged physical primaries
        if(triggeredEventMB){//triggered
            fhChargedPtMCTrigg[iVtx][iCent]->Fill(pt);
            fh2DChargedPtTrigg[iVtx][iCent]->Fill(pt, eta); 
            if(goodRecVertex){ //triggered+vertex  
                fhChargedPtMCTriggVtx[iVtx][iCent]->Fill(pt);
                fh2DChargedPtTriggVtx[iVtx][iCent]->Fill(pt, eta);
            }
        }
    }

    //--------------------------------------------------------------  
    //---- Reconstracted TRACKS 
    //--------------------------------------------------------------  
    const int nTrkCut = AliJTrackCut::kJNTrackCuts;
    if( ! triggeredEventMB ) return;
    if( ! goodRecVertex ) return;
    int nTrks   = GetJTracks()->GetEntriesFast();
    for(Int_t it = 0; it < nTrks; it++) {
        AliJTrack * track = GetJTrack(it);
        if( !track ) continue;
        bool isTriggered = track->GetFilterMap();
        if( !isTriggered ) continue;


        //== init Vars of Track
        double eta   = track->Eta();
        double etaTPC = eta;
        double etaGCG = eta;
        double ptRec = track->Pt();
        double ptTPC = track->PtTPC();
        if( fabs(ptTPC) < 1e-2 ) ptTPC=ptRec;
        else{
            TVector3 v(track->GetTPCTrack());
            etaTPC = v.Eta();
        }
        double ptGCG = track->PtGCG();
        if( fabs(ptGCG) < 1e-2 ) ptGCG=ptRec;
        else{
            TVector3 v(track->GetGCGTrack());
            etaGCG = v.Eta();
        }
        double ptMC  = -1;
        int iPrimary = kJFake;

        //== Skip tracks with out of Eta
        if( fabs(eta) > etaCut && fabs(etaTPC)>etaCut && fabs(etaGCG)>etaCut) continue;

        //== Find MC Info
        for( int imc=0;imc<nMCTrks;imc++ ){
            AliJMCTrack * mcTrack = GetJMCTrack(imc);
            if( !mcTrack || !mcTrack->IsTrue( AliJMCTrack::kPrimary ) ) continue;
            if( mcTrack && (TMath::Abs(track->GetLabel()) == TMath::Abs(mcTrack->GetLabel())) ){
                iPrimary = kJPhysicsPrimary;
                ptMC = mcTrack->Pt();
                break;
            }
        }


        //== FILL HIST
        for( int icut=0;icut<nTrkCut;icut++ ){
            if( IsSelected( track ,icut) ){
                if( fabs(eta) < etaCut ){
                    fh2DChargedPtAll[iVtx][iCent][icut]->Fill(ptRec, eta);
                    fh2DChargedPtRec[iVtx][iCent][icut]->Fill(ptMC, eta);
                    fhChargedPtMCRecoCentVtx[iVtx][iCent][icut][iPrimary][kJGlobal]->Fill( ptRec );
                    if( iPrimary == kJPhysicsPrimary ) nGenMultPri[icut]++;
                    nGenMultAll[icut]++;
                    //fhDCA2VertexXY[iVtx][iCent][icut]->Fill( track->GetDCAtoVertexXY() );
                    //fhDCA2VertexZ[iVtx][iCent][icut]->Fill( track->GetDCAtoVertexZ() );
                }
                if( fabs(etaTPC) < etaCut )
                    fhChargedPtMCRecoCentVtx[iVtx][iCent][icut][iPrimary][kJTPCOnly]->Fill( ptTPC );
                if( fabs(etaGCG) < etaCut )
                    fhChargedPtMCRecoCentVtx[iVtx][iCent][icut][iPrimary][kJGCG]->Fill( ptGCG );
                if( (icut==0 && fabs(etaTPC) < etaCut) || (icut!=0&&fabs(eta)<etaCut) )
                    fhChargedPtMCRecoCentVtx[iVtx][iCent][icut][iPrimary][kJMCTrack]->Fill( ptMC );
            }
        }
    }

    for( int icut=0;icut<nTrkCut;icut++ ){
        fh2MultGenRawPrimary[icut]->Fill( nRawMultPri, nGenMultPri[icut] );
        fh2MultGenRawAll[icut]->Fill( nRawMultPri, nGenMultAll[icut] );
    }
    //cout<<"DEBUG END AliJEfficiencyScanner::UserExec"<<endl;
}

//______________________________________________________________________________
void AliJEfficiencyScanner::Init()
{
    // Intialisation of parameters

}

//______________________________________________________________________________
void AliJEfficiencyScanner::Terminate(Option_t *)
{
    // termination
}
Commit	Line	Data
9dc4f671	1	//
	2	//#include <Riostream.h>
	3	// #include <TChain.h>
	4	// #include <TVectorT.h>
	5	// #include <TVector3.h>
	6	// #include <TFile.h>
	7	// #include <TH1.h>
	8	// #include <TObjArray.h>
	9	// #include <TObjString.h>
	10	// #include <TFormula.h>
	11	// #include <TString.h>
	12	// #include <TRefArray.h>
	13	// #include <TNtuple.h>
	14	// #include <TArrayF.h>
	15
	16	#include <TClonesArray.h>
	17	#include <TH1D.h>
	18	#include <TH2D.h>
	19
	20	#include "AliJEfficiencyScanner.h"
	21	#include "AliJTrack.h"
	22	#include "AliJMCTrack.h"
	23	#include "AliJPhoton.h"
	24	#include "AliJEventHeader.h"
	25	#include "AliJRunHeader.h"
	26	#include <TFile.h>
	27	#include <TRandom.h>
	28
	29
	30	ClassImp(AliJEfficiencyScanner);
	31
	32	//______________________________________________________________________________
	33	AliJEfficiencyScanner::AliJEfficiencyScanner() :
	34	TNamed(),
5f9b8c57	35	fMBTriggMask(0),
9dc4f671	36	fTrackList(0),
	37	fMCTrackList(0x0),
	38	fEventHeader(0x0),
	39	fRunHeader(0x0),
	40	fhVertexZMC(0x0),
	41	fhVertexZTrigg(0x0),
	42	fhVertexZTriggVtx(0x0),
	43	fhVZRawMC(0x0),
	44	fhVZRecMC(0x0),
	45	fhVZRecAccMC(0x0),
	46	fh2VtxCent(0x0),
	47	fhL0Input(0x0),
	48	fhTriggerAlice(0x0),
	49	fhZVtxMCAll(0x0),
	50	fhZVtxMCTrigg(0x0),
	51	fhZVtxMCTriggVtx(0x0),
	52	fhZVtxRecAll(0x0),
	53	fhZVtxRecTrigg(0x0),
	54	fhZVtxRecTriggVtx(0x0),
	55	fVtxReFunc(0x0),
	56	fVtxMCFunc(0x0),
	57	fVtxRatioFunc(0x0),
	58	fVtxRatioMax(1)
	59	{
	60	//Default constructor
	61	for( int ivtx=0;ivtx<kNVtxBin;ivtx++ ){
	62	for( int icent=0;icent<kNCentBin;icent++ ){
	63	fhChargedPtMCTriggVtx[ivtx][icent] = 0x0;
	64	for( int ifilter=0;ifilter<AliJTrackCut::kJNTrackCuts;ifilter++ ){
	65	for( int ipri=0;ipri<kNPrimaryStatus;ipri++ ){
	66	for( int itt=0;itt<kNTrackType;itt++ ){
	67	fhChargedPtMCRecoCentVtx[ivtx][icent][ifilter][ipri][itt] = 0x0;
	68	}
	69	}
	70	}
	71	}
	72	}
	73
	74	}
	75
	76	//______________________________________________________________________________
	77	AliJEfficiencyScanner::AliJEfficiencyScanner(const char *name):
	78	TNamed(name,name),
5f9b8c57	79	fMBTriggMask(0),
9dc4f671	80	fTrackList(0),
	81	fMCTrackList(0x0),
	82	fEventHeader(0x0),
	83	fRunHeader(0x0),
	84	fhVertexZMC(0x0),
	85	fhVertexZTrigg(0x0),
	86	fhVertexZTriggVtx(0x0),
	87	fhVZRawMC(0x0),
	88	fhVZRecMC(0x0),
	89	fhVZRecAccMC(0x0),
	90	fh2VtxCent(0x0),
	91	fhL0Input(0x0),
	92	fhTriggerAlice(0x0),
	93	fhZVtxMCAll(0x0),
	94	fhZVtxMCTrigg(0x0),
	95	fhZVtxMCTriggVtx(0x0),
	96	fhZVtxRecAll(0x0),
	97	fhZVtxRecTrigg(0x0),
	98	fhZVtxRecTriggVtx(0x0),
	99	fVtxReFunc(0x0),
	100	fVtxMCFunc(0x0),
	101	fVtxRatioFunc(0x0),
	102	fVtxRatioMax(1)
	103	{
	104	// Constructor
	105	for( int ivtx=0;ivtx<kNVtxBin;ivtx++ ){
	106	for( int icent=0;icent<kNCentBin;icent++ ){
	107	fhChargedPtMCTriggVtx[ivtx][icent] = 0x0;
	108	for( int ifilter=0;ifilter<AliJTrackCut::kJNTrackCuts;ifilter++ ){
	109	for( int ipri=0;ipri<kNPrimaryStatus;ipri++ ){
	110	for( int itt=0;itt<kNTrackType;itt++ ){
	111	fhChargedPtMCRecoCentVtx[ivtx][icent][ifilter][ipri][itt] = 0x0;
	112	}
	113	}
	114	}
	115	}
	116	}
	117	}
	118
	119	//____________________________________________________________________________
	120	AliJEfficiencyScanner::AliJEfficiencyScanner(const AliJEfficiencyScanner& ap) :
	121	TNamed(ap.GetName(), ap.GetTitle()),
5f9b8c57	122	fMBTriggMask(ap.fMBTriggMask),
9dc4f671	123	fTrackList(ap.fTrackList),
	124	fMCTrackList(ap.fMCTrackList),
	125	fEventHeader( ap.fEventHeader ),
	126	fRunHeader(ap.fRunHeader),
	127	fhVertexZMC(ap.fhVertexZMC),
	128	fhVertexZTrigg(ap.fhVertexZTrigg),
	129	fhVertexZTriggVtx(ap.fhVertexZTriggVtx),
	130	fhVZRawMC(ap.fhVZRawMC),
	131	fhVZRecMC(ap.fhVZRecMC),
	132	fhVZRecAccMC(ap.fhVZRecAccMC),
	133	fh2VtxCent(ap.fh2VtxCent),
	134	fhL0Input(ap.fhL0Input),
	135	fhTriggerAlice(ap.fhTriggerAlice),
	136	fhZVtxMCAll(ap.fhZVtxMCAll),
	137	fhZVtxMCTrigg(ap.fhZVtxMCTrigg),
	138	fhZVtxMCTriggVtx(ap.fhZVtxMCTriggVtx),
	139	fhZVtxRecAll(ap.fhZVtxRecAll),
	140	fhZVtxRecTrigg(ap.fhZVtxRecTrigg),
	141	fhZVtxRecTriggVtx(ap.fhZVtxRecTriggVtx),
	142	fVtxReFunc(ap.fVtxReFunc),
	143	fVtxMCFunc(ap.fVtxMCFunc),
	144	fVtxRatioFunc(ap.fVtxRatioFunc),
	145	fVtxRatioMax(ap.fVtxRatioMax)
	146	{
	147	// cpy ctor
	148	for( int ivtx=0;ivtx<kNVtxBin;ivtx++ ){
	149	for( int icent=0;icent<kNCentBin;icent++ ){
	150	fhChargedPtMCTriggVtx[ivtx][icent] = ap.fhChargedPtMCTriggVtx[ivtx][icent];
	151	for( int ifilter=0;ifilter<AliJTrackCut::kJNTrackCuts;ifilter++ ){
	152	for( int ipri=0;ipri<kNPrimaryStatus;ipri++ ){
	153	for( int itt=0;itt<kNTrackType;itt++ ){
	154	fhChargedPtMCRecoCentVtx[ivtx][icent][ifilter][ipri][itt] = ap.fhChargedPtMCRecoCentVtx[ivtx][icent][ifilter][ipri][itt];
	155	}
	156	}
	157	}
	158	}
	159	}
	160
	161	}
	162
	163	//_____________________________________________________________________________
	164	AliJEfficiencyScanner& AliJEfficiencyScanner::operator = (const AliJEfficiencyScanner& ap)
	165	{
	166	// assignment operator
	167
	168	this->~AliJEfficiencyScanner();
	169	new(this) AliJEfficiencyScanner(ap);
	170	return *this;
	171	}
	172
	173	//______________________________________________________________________________
	174	AliJEfficiencyScanner::~AliJEfficiencyScanner()
	175	{
	176	// destructor
	177	}
	178
	179	//________________________________________________________________________
	180
	181	void AliJEfficiencyScanner::UserCreateOutputObjects()
	182	{
	183	//=== create the jcorran outputs objects
	184	cout<<"DEBUG Start AliJEfficiencyScanner::UserCreateOutputObjects() "<<"\t"<<gDirectory<<endl;
	185
	186	double ptbin[300] = {0};
187	double pt = 0;
188	int i = 0;
189	for( i=0;i<300;i++ ){
190	ptbin[i] = pt;
191	if( pt > 50 ) break;
192	if( pt < 3 ) pt+= 0.05;
193	else if( pt < 5 ) pt+= 0.1;
194	else if( pt < 10 ) pt+= 1;
195	else pt+= 1;
196	}
197	cout<<"n Ptbin = "<<i<<endl;
198	int nPtBin = i-1;
199
200	fh2VtxCent = new TH2D("h2VtxCent","h2VtxCent",100,0,10,110,0,110 );
201	fh2VtxCent->SetDirectory(gDirectory);
202
203	int nVtxBin = kNVtxBin;
204	int nCentBin = kNCentBin;
205	if( fRunHeader->IsPP() ) nCentBin = 1;
206
207	TString name ="";
208	for( int ivtx=0;ivtx<nVtxBin;ivtx++ ){
209	for( int icent=0;icent<nCentBin;icent++ ){
210	fhChargedPtMC[ivtx][icent] = AddTH1D( Form("hChargedPtMC%02d%02d",ivtx,icent),
211	new TH1D("","", nPtBin, ptbin) );
212	fhChargedPtMCTriggVtx[ivtx][icent]=AddTH1D(
213	Form("hChargedPtMCTriggVtx%02d%02d", ivtx, icent ),
214	new TH1D("", "", nPtBin, ptbin) );
215	fhChargedPtMCTrigg[ivtx][icent]=AddTH1D(
216	Form("hChargedPtMCTrigg%02d%02d", ivtx, icent ),
217	new TH1D("", "", nPtBin, ptbin) );
218	//for any MC track filled with MC pt in triggered event
219	name = Form("h2DChargedPtTrigg%02d%02d",ivtx,icent);
220	fh2DChargedPtTrigg[ivtx][icent]=AddTH2D(name, new TH2D(name,name, nPtBin, ptbin, 20, -0.8, 0.8));
221	//for any MC track filled with MC pt in triggered event with rec vertex
222	name = Form("h2DChargedPtTriggVtx%02d%02d",ivtx,icent);
223	fh2DChargedPtTriggVtx[ivtx][icent]=AddTH2D(name, new TH2D(name,name, nPtBin, ptbin, 20, -0.8, 0.8));
224
225	}
226	}
227	fhTriggerAlice = AddTH1D("hTriggerAlice",new TH1D("","",32,0,32));
228	fhL0Input = AddTH1D("hL0Input",new TH1D("","",32,0,32));
229	// DCA bin
230	double dcaBin[1000];
231	double dbin=-50-1;
232	int ndbin = 0;
233	dcaBin[0] = dbin;
234	double tol = 1e-5;
235	while(dbin < 50-1){
236	if( fabs(dbin) < 2-tol ) dbin+=0.01;
237	else if( fabs(dbin) < 5-tol ) dbin+= 0.05;
238	else if( fabs(dbin) <= 10-tol ) dbin+= 0.1;
239	else dbin += 1;
240	if( fabs(dbin) < tol ) dbin=0;
241	dcaBin[ndbin++] = dbin;
242	}
243
244	for( int ivtx=0;ivtx<nVtxBin;ivtx++ ){
245	for( int icent=0;icent<nCentBin;icent++ ){
246	for( int ifilter=0;ifilter<AliJTrackCut::kJNTrackCuts;ifilter++ ){
247
248	name = Form("h2DChargedPtAll%02d%02d%02d",ivtx,icent,ifilter);
249	fh2DChargedPtAll[ivtx][icent][ifilter]=AddTH2D(name, new TH2D(name,name, nPtBin, ptbin, 20, -0.8, 0.8));
250	name = Form("h2DChargedPtRec%02d%02d%02d",ivtx,icent,ifilter);
251	fh2DChargedPtRec[ivtx][icent][ifilter]=AddTH2D(name, new TH2D(name,name, nPtBin, ptbin, 20, -0.8, 0.8));
252	for( int ipri=0;ipri<kNPrimaryStatus;ipri++ ){
253	for( int itt=0;itt<kNTrackType;itt++ ){
254	name = Form("hChargedPtMCRecoCentVtx%02d%02d%02d%02d%02d", ivtx, icent, ifilter, ipri, itt );
255	fhChargedPtMCRecoCentVtx[ivtx][icent][ifilter][ipri][itt]=AddTH1D(name, new TH1D( name, name, nPtBin, ptbin));
256	}
257	}
258	name = Form("hDCA2VertexXY%02d%02d%02d", ivtx, icent, ifilter );
259	fhDCA2VertexXY[ivtx][icent][ifilter] = AddTH1D(name, new TH1D("","",ndbin-1, dcaBin ));
260	name = Form("hDCA2VertexZ%02d%02d%02d", ivtx, icent, ifilter );
261	fhDCA2VertexZ[ivtx][icent][ifilter] = AddTH1D( name, new TH1D("","",ndbin-1, dcaBin ));
262
263	}
264	}
265	}
266	// TH2D * fh2MultGenRawPrimary[AliJTrackCut::kJNTrackCuts];
267	// TH2D * fh2MultGenRawAll[AliJTrackCut::kJNTrackCuts];
268
269	// TH1D * fhDCA2VertexXY[kNVtxBin][kNCentBin][AliJTrackCut::kJNTrackCuts];
270	for( int ifilter=0;ifilter<AliJTrackCut::kJNTrackCuts;ifilter++ ){
271	name = Form("h2MultGenRawPrimary%02d", ifilter);
272	fh2MultGenRawPrimary[ifilter] = AddTH2D( name, new TH2D(name, name, 100,0,100,100,0,100 ));
273	name = Form("h2MultGenRawAll%02d", ifilter);
274	fh2MultGenRawAll[ifilter] = AddTH2D( name, new TH2D(name, name, 100,0,100,100,0,100 ));
275	}
276	double binsVertexMult[] = {0,1,2,3,4,5,10000};
277	int NbinsVertexMult = sizeof(binsVertexMult)/sizeof(double)-1;
278	double binsVertexZ[] = {-10,-6,-3,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,3,6,10};
279	int NbinsVertexZ = sizeof(binsVertexZ)/sizeof(double)-1;
280	fhVertexZMC = AddTH2D("hVertexZMC",new TH2D("hVertexZMC","hVertexZMC", NbinsVertexMult, binsVertexMult, NbinsVertexZ, binsVertexZ));
281	fhVertexZTrigg = AddTH2D("hVertexZTrigg",new TH2D("hVertexZTrigg","hVertexZTrigg", NbinsVertexMult, binsVertexMult, NbinsVertexZ, binsVertexZ));
282	fhVertexZTriggVtx = AddTH2D("hVertexZTriggVtx",new TH2D("hVertexZTriggVtx","hVertexZTriggVtx", NbinsVertexMult, binsVertexMult, NbinsVertexZ, binsVertexZ));
283	fhVZRawMC = AddTH1D("hVtxZMC", new TH1D("hVtxZMC","VertexZ in MC ",200,-50,50));
284	fhVZRecMC = AddTH1D("hVtxZRec", new TH1D("hVtxZRec","VertexZ Rec in MC",200,-50,50));
285	fhVZRecAccMC = AddTH1D("hVtxZRecAcc", new TH1D("hVtxZRecAcc","VertexZ Rec in MC and acc",200,-50,50));
286
287	int NvtxBins = 400;
288	double vtxedge = 50;
289	fhZVtxMCAll = AddTH1D("hZVtxMCAll", new TH1D("hZVtxMCAll","VetrexZ in MC for all event",NvtxBins, -vtxedge, vtxedge ));
290	fhZVtxMCTrigg = AddTH1D("hZVtxMCTrigg", new TH1D("hZVtxMCTrigg","VetrexZ in MC for triggered event",NvtxBins, -vtxedge, vtxedge ));
291	fhZVtxMCTriggVtx = AddTH1D("hZVtxMCTriggVtx", new TH1D("hZVtxMCTriggVtx","VetrexZ in MC for triggered and vtx reconstructed event",NvtxBins, -vtxedge, vtxedge ));
292
293	fhZVtxRecAll = AddTH1D("hZVtxRecAll", new TH1D("hZVtxRecAll","VetrexZ in Rec for all event",NvtxBins, -vtxedge, vtxedge ));
294	fhZVtxRecTrigg = AddTH1D("hZVtxRecTrigg", new TH1D("hZVtxRecTrigg","VetrexZ in Rec for triggered event",NvtxBins, -vtxedge, vtxedge ));
295	fhZVtxRecTriggVtx = AddTH1D("hZVtxRecTriggVtx", new TH1D("hZVtxRecTriggVtx","VetrexZ in Rec for triggered and vtx reconstructed event",NvtxBins, -vtxedge, vtxedge ));
296
297	double v0 = -20;
298	double v1 = 20;
299	fVtxReFunc = new TF1("VtxReFunc", "gaus",v0,v1);
300	fVtxReFunc->SetParameters(1 , -7.14076e-01, 6.27110 );
301	fVtxMCFunc = new TF1("VtxMCFunc", "gaus",v0,v1);
302	fVtxMCFunc->SetParameters(1, -4.53674e-01, 5.27088e+00 );
303	fVtxRatioFunc = new TF1("VtxRatioFunc", "VtxReFunc/VtxMCFunc",v0,v1);
304	fVtxRatioMax = fVtxRatioFunc->GetMaximum();
305
306
307
308	cout << "Add(fAliJRunHeader) in UserCreateObject() ======= " << endl;
309	cout<<"DEBUG END AliJEfficiencyScanner::UserCreateOutputObjects() "<<endl;
310
311	}
312
313	//______________________________________________________________________________
314
315	void AliJEfficiencyScanner::UserExec(Option_t *option)
316	{
317
318	// Main loop
319	// Called for each event
320	JUNUSED(option);
321
322	double zVtxCut = 10;
323	double etaCut = 0.8;
324	double nCentBin = 20;
325
326	AliJEventHeader * eventHeader = GetJEventHeader();
327
328	//== TRIGGER( MB ) TODO is this?
5f9b8c57	329	Bool_t triggeredEventMB = eventHeader->GetTriggerMaskJCorran() & fMBTriggMask; // 1:kMB 8:kINT7 //TODO
9dc4f671	330	UInt_t trigAlice = eventHeader->GetTriggerMaskAlice();
	331	/*
	332	cout<<"Trigger"<<endl;
	333	for( int i=0;i<32 ;i++ ) cout<<Form("%02d",i)<<" ";
	334	cout<<endl;
	335	for( int i=0;i<32 ;i++ ) cout<<(eventHeader->GetTriggerMaskJCorran()&(1<<i)?1:0)<<" ";
	336	cout<<endl;
	337	for( int i=0;i<32 ;i++ ) cout<<(eventHeader->GetTriggerMaskAlice()&(1<<i)?1:0)<<" ";
	338	cout<<endl;
	339	*/
	340	/*
	341	bool isMB1 = trigAlice & (1UL<<0);
	342	bool isMBBG1 = trigAlice & (1UL<<3);
	343	bool isGFO = trigAlice & (1UL<<14);
	344	bool isMBBG3 = trigAlice & (1UL<<1);
	345	bool isV0A = trigAlice & (1UL<<10);
	346	bool isV0C = trigAlice & (1UL<11);
	347	*/
	348
	349	// L0Input
	350	UInt_t l0Input = eventHeader->GetL0TriggerInputs();
	351	/*
	352	bool isL0V0A = l0Input & (1UL<<9);
	353	bool isL0V0C = l0Input & (1UL<<10);
	354	bool isL0V0AND = l0Input & (1UL<<4);
	355	bool isL0V0OR = l0Input & (1UL<<1);
	356	bool isL0V0BeamGas = l0Input & (1UL<<2);
	357	bool isL0SMB = l0Input & (1UL<<3);
	358	bool isL0T0X = l0Input & (1UL<<0);
	359	*/
	360
	361
	362
	363	for( UInt_t i=0;i<32;i++ ){
	364	if( l0Input & 1UL<<i ) fhL0Input->Fill(i);
	365	if( trigAlice & 1UL<<i ) fhTriggerAlice->Fill(i);
	366	/*
	367	if( isL0V0A != isV0A ) cout<<"J_WARN 0 : L0:V0A != Class:V0A "<<isL0V0A<<"\t"<<isV0A<<endl;
	368	if( isL0V0C != isV0C ) cout<<"J_WARN 1 : L0:V0C != Class:V0C "<<isL0V0C<<"\t"<<isV0C<<endl;
	369	if( isL0SMB != isGFO ) cout<<"J_WARN 2 : L0:SMB != Class:GFO "<<isL0SMB<<"\t"<<isGFO<<endl;
	370	if( isL0V0OR != isMBBG1 ) cout<<"J_WARN 3 : L0:V0OR != Class:MBBG1 "<<isL0V0OR<<"\t"<<isMBBG1<<endl;
	371	if( (isL0V0A\|\|isL0V0C) != isMBBG1 ) cout<<"J_WARN 4 : L0:V0A\|V0C != Class:MBBG1 "<<(isL0V0A\|\|isL0V0C)<<"\t"<<isMBBG1<<endl;
	372	if( ( isL0T0X \|\| isL0V0OR \|\| isL0SMB && !isL0V0BeamGas ) != isMB1 ) cout<<"J_WARN 5 : L0:MB1 != Class:MB1 "<<isMB1<<endl;
	373	if( (isL0V0OR \|\| isL0SMB && !isL0V0BeamGas ) != (( isMBBG1 \|\| isGFO ) && isMB1 )) cout<<"J_WARN 6 : L0:MB != Class:MB "<<(isL0V0OR \|\| isL0SMB && !isL0V0BeamGas)<<"\t"<<((isMBBG1 \|\| isGFO ) && isMB1)<<endl;
	374	if( (isL0V0OR \|\| isL0SMB ) != (( isMBBG1 \|\| isGFO ) )) cout<<"J_WARN 6 : L0:MB != Class:MB "<<(isL0V0OR \|\| isL0SMB && !isL0V0BeamGas)<<"\t"<<((isMBBG1 \|\| isGFO ) && isMB1)<<endl;
	375	*/
	376
	377	}
	378	//* FOR TEST */
	379	//Bool_t triggeredEventMB = ( isMBBG1 \|\| isGFO ) && isMB1;
	380	//Bool_t triggeredEventMB = isMB1;
	381	//Bool_t triggeredEventMB = ( isMBBG1 \|\| isGFO );
	382	//
	383	//== REAL MB with L0 INput
	384	//triggeredEventMB = ( (isL0V0OR \|\| isL0SMB) );
	385	//== CENTRALITY
	386	double centrality = eventHeader->GetCentrality();
	387	int iCent = int( centrality/(100./nCentBin));
	388	if( iCent< 0 \|\| iCent>20 ) return;
	389	//if( centrality < 0 \|\| centrality > 100 ) return;
5f9b8c57	390	if( fRunHeader->IsPP() ) {iCent = 0; centrality=0;}
9dc4f671	391
	392	//== Vertex
	393	//==== Reco
	394	int ncontributors = eventHeader->GetVtxMult();
	395	double zVtxRec = eventHeader->GetZVertex();
	396	Bool_t goodRecVertex = (ncontributors>0) && (fabs(zVtxRec)<=zVtxCut);
	397	int iVtx = 0;
	398	//==== MC
	399	double zVtxMC = eventHeader->GetZVertexMC();
	400	//==== vtx sampling
	401	/*
	402	if( zVtxMC > fVtxRatioFunc->GetXmin() && zVtxMC < fVtxRatioFunc->GetXmax() ){
	403	double vtxRatio = fVtxRatioFunc->Eval(zVtxMC)/fVtxRatioMax;
	404	double random = gRandom->Uniform(0,1);
	405	if( random > vtxRatio ) return;
	406	}
	407	*/
	408
	409	fhVZRawMC->Fill(zVtxMC);
	410	fhVertexZMC->Fill(ncontributors,zVtxMC);
	411	// for crosscheck MC input events
	412	fhZVtxMCAll->Fill(zVtxMC);
	413	fhZVtxRecAll->Fill(zVtxRec);
	414
	415	if(triggeredEventMB){
	416	fhVertexZTrigg->Fill(ncontributors,zVtxMC);
	417	// for crosscheck passing trigger condition
	418	fhZVtxMCTrigg->Fill(zVtxMC);
	419	fhZVtxRecTrigg->Fill(zVtxRec);
	420
	421	if(ncontributors>0){//reconstructed vertex
	422	fhVZRecMC->Fill(zVtxRec);
	423	// for crosscheck passing trigger condition + reconstruted vtx
	424	fhZVtxMCTriggVtx->Fill(zVtxMC);
	425	fhZVtxRecTriggVtx->Fill(zVtxRec);
	426	}
	427	if(goodRecVertex){
	428	fhVertexZTriggVtx->Fill(ncontributors,zVtxMC);
	429	fhVZRecAccMC->Fill(zVtxRec);
	430	}
	431	}
	432
	433	fh2VtxCent->Fill( zVtxRec, centrality );
	434
	435	// TODO ? MC or REC?
	436
	437	int nRawMultPri = 0;
	438	int nGenMultPri[AliJTrackCut::kJNTrackCuts] = {0};
	439	int nGenMultAll[AliJTrackCut::kJNTrackCuts] = {0};
	440
	441	//==============================================
	442	// LOOP over MC
	443	//==============================================
	444
	445	//if( fabs(zVtxMC) > zVtxCut ) return;
	446	int nMCTrks = 0;
	447	if( fRunHeader->IsMC() )
	448	nMCTrks = GetJMCTracks()->GetEntriesFast();
	449	for( int it=0; it< nMCTrks; it++ ){
	450	AliJMCTrack * track = GetJMCTrack( it );//Always IsPhysicalPrimary
	451	if( !track ) continue;
	452	if( !track->IsTrue( AliJMCTrack::kPrimary ) ) continue;
	453	// TODO need? if( ! track->IsFinal() ) continue;
	454	double eta = track->Eta();
455	double pt = track->Pt();
456	if( fabs(eta) > etaCut ) continue;
457	if( ! track->IsCharged() ) continue;
458
459
460	nRawMultPri++;
461	fhChargedPtMC[iVtx][iCent]->Fill(pt); //ALL charged physical primaries
462	if(triggeredEventMB){//triggered
463	fhChargedPtMCTrigg[iVtx][iCent]->Fill(pt);
464	fh2DChargedPtTrigg[iVtx][iCent]->Fill(pt, eta);
465	if(goodRecVertex){ //triggered+vertex
466	fhChargedPtMCTriggVtx[iVtx][iCent]->Fill(pt);
467	fh2DChargedPtTriggVtx[iVtx][iCent]->Fill(pt, eta);
468	}
469	}
470	}
471
472	//--------------------------------------------------------------
473	//---- Reconstracted TRACKS
474	//--------------------------------------------------------------
475	const int nTrkCut = AliJTrackCut::kJNTrackCuts;
476	if( ! triggeredEventMB ) return;
477	if( ! goodRecVertex ) return;
478	int nTrks = GetJTracks()->GetEntriesFast();
479	for(Int_t it = 0; it < nTrks; it++) {
480	AliJTrack * track = GetJTrack(it);
481	if( !track ) continue;
482	bool isTriggered = track->GetFilterMap();
483	if( !isTriggered ) continue;
484
485
486	//== init Vars of Track
487	double eta = track->Eta();
488	double etaTPC = eta;
489	double etaGCG = eta;
490	double ptRec = track->Pt();
491	double ptTPC = track->PtTPC();
492	if( fabs(ptTPC) < 1e-2 ) ptTPC=ptRec;
493	else{
494	TVector3 v(track->GetTPCTrack());
495	etaTPC = v.Eta();
496	}
497	double ptGCG = track->PtGCG();
498	if( fabs(ptGCG) < 1e-2 ) ptGCG=ptRec;
499	else{
500	TVector3 v(track->GetGCGTrack());
501	etaGCG = v.Eta();
502	}
503	double ptMC = -1;
504	int iPrimary = kJFake;
505
506	//== Skip tracks with out of Eta
507	if( fabs(eta) > etaCut && fabs(etaTPC)>etaCut && fabs(etaGCG)>etaCut) continue;
508
509	//== Find MC Info
510	for( int imc=0;imc<nMCTrks;imc++ ){
511	AliJMCTrack * mcTrack = GetJMCTrack(imc);
512	if( !mcTrack \|\| !mcTrack->IsTrue( AliJMCTrack::kPrimary ) ) continue;
513	if( mcTrack && (TMath::Abs(track->GetLabel()) == TMath::Abs(mcTrack->GetLabel())) ){
514	iPrimary = kJPhysicsPrimary;
515	ptMC = mcTrack->Pt();
516	break;
517	}
518	}
519
520
521	//== FILL HIST
522	for( int icut=0;icut<nTrkCut;icut++ ){
523	if( IsSelected( track ,icut) ){
524	if( fabs(eta) < etaCut ){
525	fh2DChargedPtAll[iVtx][iCent][icut]->Fill(ptRec, eta);
526	fh2DChargedPtRec[iVtx][iCent][icut]->Fill(ptMC, eta);
527	fhChargedPtMCRecoCentVtx[iVtx][iCent][icut][iPrimary][kJGlobal]->Fill( ptRec );
528	if( iPrimary == kJPhysicsPrimary ) nGenMultPri[icut]++;
529	nGenMultAll[icut]++;
530	//fhDCA2VertexXY[iVtx][iCent][icut]->Fill( track->GetDCAtoVertexXY() );
531	//fhDCA2VertexZ[iVtx][iCent][icut]->Fill( track->GetDCAtoVertexZ() );
532	}
533	if( fabs(etaTPC) < etaCut )
534	fhChargedPtMCRecoCentVtx[iVtx][iCent][icut][iPrimary][kJTPCOnly]->Fill( ptTPC );
535	if( fabs(etaGCG) < etaCut )
536	fhChargedPtMCRecoCentVtx[iVtx][iCent][icut][iPrimary][kJGCG]->Fill( ptGCG );
537	if( (icut==0 && fabs(etaTPC) < etaCut) \|\| (icut!=0&&fabs(eta)<etaCut) )
538	fhChargedPtMCRecoCentVtx[iVtx][iCent][icut][iPrimary][kJMCTrack]->Fill( ptMC );
539	}
540	}
541	}
542
543	for( int icut=0;icut<nTrkCut;icut++ ){
544	fh2MultGenRawPrimary[icut]->Fill( nRawMultPri, nGenMultPri[icut] );
545	fh2MultGenRawAll[icut]->Fill( nRawMultPri, nGenMultAll[icut] );
546	}
547	//cout<<"DEBUG END AliJEfficiencyScanner::UserExec"<<endl;
548	}
549
550	//______________________________________________________________________________
551	void AliJEfficiencyScanner::Init()
552	{
553	// Intialisation of parameters
554
555	}
556
557	//______________________________________________________________________________
558	void AliJEfficiencyScanner::Terminate(Option_t *)
559	{
560	// termination
561	}
562
563