#ifndef __MAKECINT__ 
 #ifndef __CINT__
  #include "alles.h"
  #include <TPDGCode.h>
  #include "AliRunAnalysis.h"
  #include "AliHBTAnalysis.h"
  #include "AliReader.h"
  #include "AliReaderESD.h"
  #include "AliReaderKineTree.h"
  #include "AliParticleCut.h"
  #include "AliAOD.h"
  #include "AliAODPairCut.h"
  #include "AliAODParticleCut.h"
  #include "AliHBTQResolutionFctns.h"
  #include "AliHBTQDistributionFctns.h"
  #include "AliHBTMonDistributionFctns.h"
  #include "AliHBTTwoTrackEffFctn.h"
  #include "AliHBTCorrelFctn.h"
  #include "AliHBTWeightFctn.h"
  #include "AliHBTWeightTheorFctn.h"
  #include "AliHBTWeights.h"
  #include "AliHBTLLWeights.h"
  #include "TSystem.h"
  #include "TObjString.h"
  #include "TString.h"
  #include "AliPDG.h"
  #include "AliHBTMonDistributionFctns.h"
  #include "AliHBTMonResolutionFctns.h"
 #endif 
#endif

//This macro can not be executed from 

AliHBTAnalysis* buildanalysis(Int_t pid1, Int_t pid2, Float_t qmax, Bool_t threeDcuts, Bool_t weights, const char* anasuffix = "",const char* titlesfx = "");


void nonidhbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles", 
                Int_t first = -1,Int_t last = -1, 
                char *prefname = "", char *datname = "", 
                Int_t readfrom = 0, Int_t readto = 0,const char* dir = 0, const char* run = 0)
 {

::Info("nonidanalysis","New ANALYSIS schema");
AliHBTAnalysis* buildanalysis(Int_t pid1, Int_t pid2, Float_t qmax, Bool_t threeDcuts, Bool_t weights, const char* anasuffix = "",const char* titlesfx = "");

//HBT Anlysis Macro
//Anlyzes TPC recontructed tracks and simulated particles that corresponds to them

//datatype defines type of data to be read
//  Kine  - analyzes Kine Tree: simulated particles
//  TPC   - analyzes TPC   tracking + particles corresponding to these tracks
//  ITSv1 - analyzes ITSv1 ----------------------//--------------------------
//  ITSv2 - analyzes ITSv2 ----------------------//--------------------------

//processopt defines option passed to AliHBTAnlysis::Process method
// default: TracksAndParticles - process both recontructed tracks and sim. particles corresponding to them 
//          Tracks - process only recontructed tracks
//          Particles - process only simulated particles

//Reads data from diroctories from first to last(including)
// For examples if first=3 and last=5 it reads from
//  ./3/
//  ./4/
//  ./5/
//if first or last is negative (or both), it reads from current directory
//
//these names I use when analysis is done directly from CASTOR files via RFIO
//This switches control to analyse:  
  Bool_t ident = kFALSE;//identical particles 
  Bool_t nonident = kTRUE;//non identical analysis
  Bool_t weights = kTRUE;//non identical analysis
  
//  AliHBTParticle::SetDebug(0);
  
  char* directory[200];
  const char* basedir="rfio:/afs/cern.ch/user/s/skowron";
  const char* field="wrk";
  const char* serie ="dane";
  
  if (dir == 0)
   { 
     sprintf(directory,"%s/%s/%s",basedir,field,serie);
   }
  else
   {
     sprintf(directory,"%s",dir);
   } 
//  const char* basedir=".";
//  const char* serie=".";
//  const char* field = ".";
 
//Here are specified particles to be analysed
//Last should be 0 - analysis relies on it
//User can add any number of particles
// Dynamically link some shared libs                    
//  gROOT->LoadMacro("loadlibs.C");
//  loadlibs();
  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libANALYSIS");
  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libHBTAN");

  Int_t PID[10];

  PID[0]=kKPlus;
  PID[1]=kKMinus;
  PID[2]=0;//Last must be 0!!!!!!!!!!!!!!!!!!
  
//  PID[0]=kPiMinus;
//  PID[0]=kKPlus;
//  PID[1]=kKMinus;
//  PID[2]=kProton;
//  PID[1]=0;//Last must be 0!!!!!!!!!!!!!!!!!!
  
//  PID[0]=kProton;
//  PID[1]=kProtonBar;
//  PID[2]=kKPlus;
//  PID[3]=kKMinus;
//  PID[4]=kPiPlus;
//  PID[5]=kPiMinus;
//  PID[6]=0;//Last must be 0!!!!!!!!!!!!!!!!!!
  Float_t PIDprob[10];
  PIDprob[0] = 0.5;
  PIDprob[1] = 0.5;
  PIDprob[2] = 0.0;
  PIDprob[3] = 0.0;
  PIDprob[4] = 0.0;
  PIDprob[5] = 0.0;
  PIDprob[6] = 0.0;
  PIDprob[7] = 0.0;
  PIDprob[8] = 0.0;
  PIDprob[9] = 0.0;
 
  
  /***********************************************************/
  //Create analysis and reader
  
  TString dataname;
  AliReader* reader;
  Int_t kine = strcmp(datatype,"Kine");
  Int_t ESD = strcmp(datatype,"ESD");
  Int_t intern = strcmp(datatype,"Intern");
  
  if(!kine)
   {
    reader = new AliReaderKineTree();
    processopt="Particles"; //this reader by definition reads only simulated particles
    dataname = "Kine";
   }
  else if(!ESD)
   {
    AliReaderESD* esdreader = new AliReaderESD();
//    esdreader->SetNumberOfTrackPoints(5,30.);//5 points every 30 cm
//    esdreader->SetITSTrackPoints(kTRUE);
//    esdreader->SetClusterMap();


    if (strstr(processopt,"Particles"))
      esdreader->ReadSimulatedData(kTRUE);
    else
      esdreader->ReadSimulatedData(kFALSE);

    reader = esdreader;
    dataname = "ESD";
//    esdreader->MustTPC(kTRUE);
    
//    esdreader->SetChi2Range(0.0,3.0);
//    esdreader->SetC22Range(0.0,3e-6);
//    esdreader->SetC33Range(0.0,10e-6);
//    esdreader->SetC44Range(0.0,0.0012);
    
   }
  else if(!intern)
   {
    dataname = datname;
    cerr<<"Option "<<datatype<<"  AOD Reader not yet implemented. Exiting"<<endl;
    return;
//    reader = new AliHBTReaderInternal(dataname+".root");
   }
  else
   {
    cerr<<"Option "<<datatype<<"  not recognized. Exiting"<<endl;
    return;
   }

  cout<<"hbtanalysis.C: Reader Set, "<<reader->Class()->GetName()<<" dataname: "<<dataname<<endl;
  cout<<"hbtanalysis.C: Setting dirs\n";
  TObjArray* dirs=0;
  if ( ((first >= 0) && (last >= 0) ) && ( (last-first)>=0 ) )
  {//read from many dirs dirs
    char buff[50];
    dirs = new TObjArray(last-first+1);
    for (Int_t i = first; i<=last; i++)
     {
//     sprintf(buff,"%s/%s/%s/%05.5d",basedir,field,serie,i);
       printf("%s/%d\n",directory,i);
       sprintf(buff,"%s/%d",directory,i);
       TObjString *odir= new TObjString(buff);
       dirs->Add(odir);
     }
   }
  else
   {
     dirs = new TObjArray(1);
     TObjString *odir= new TObjString("rfio:.");
     dirs->Add(odir);
   } 
  reader->SetDirs(dirs);
  reader->SetName("READER");
  reader->SetTitle("Reader Used in HBT Analysis");
  cout<<"hbtanalysis.C: Dirs set\n";
  //reader->SetBlend(kTRUE);
  reader->ReadEventsFromTo(readfrom,readto);
  /***********************************************************/    
  /*****   R E A D E R ***************************************/    
  /***********************************************************/    
  
  //we want have only low pt pi+ so set a cut to reader
  AliAODParticleCut* readerpartcut= new AliAODParticleCut();
  readerpartcut->SetPtRange(0.0,100000.0);
  readerpartcut->SetVxRange(-0.004,0.004);//select particles from vertex
  readerpartcut->SetVyRange(-0.004,0.004);//select particles from vertex
  readerpartcut->SetVzRange(-0.004,0.004);//select particles from vertex
//  readerpartcut->SetThetaRange(0.25*TMath::Pi(), 0.75*TMath::Pi());//45->135 deg.
  
//  readerpartcut->SetYRange(-1.0,1.0);//!!!!!!!!!!!!!!!!!!!!
  Int_t l = 0;
  while (PID[l] != 0)
   {
     cout<<"hbtanalysis.C: Adding PID  = "<<PID[l]<<" l = "<<l<<endl;
     readerpartcut->SetPID(PID[l]);
     
     AliAODParticleCut * pcut = (AliAODParticleCut*)readerpartcut->Clone();
     AliAODPIDCut* pidcut = new AliAODPIDCut(PID[l],PIDprob[l]);
     pcut->AddBasePartCut(pidcut);
     reader->AddParticleCut(pcut);//read this particle type with this cut
     delete pcut;
     
     l++;
   }
  cout<<"hbtanalysis.C: PID Cut set in Reader "<<endl;
  //*******************
  //R A N D O M I Z E R
  //*******************
  AliHBTPositionRandomizer* rndmreader = new AliHBTPositionRandomizer(reader);
  rndmreader->SetGaussianBall(8.0);
  reader = rndmreader;
  
  TString filename = prefname + dataname + ".anal.";
  TString fn;
  TParticlePDG* pdg1;  
  TParticlePDG* pdg2;  
  TString p2name;  
  TString p1name ;
  
  AliHBTAnalysis* hbtanalysis;

  /************************************************************/
  /****   P R O C E S S                 ***********************/
  /************************************************************/

  cout<<"Process 1 "<<endl;

  l = 0;
  Int_t k;
  while (PID[l] != 0)
   {
     cout<<"Process l "<<l<<endl;
     
     cout<<"TDatabasePDG::Instance()->GetParticle(PID[l])\n";
     pdg1 = TDatabasePDG::Instance()->GetParticle(PID[l]);
     cout<<"TDatabasePDG::Instance()->GetParticle(PID[l]) Done\n";
     
     if (pdg1 == 0x0)
       {
         cerr<<"ERROR: hbtanalysis.C: Can not find particle type "<<PID[l]<<" in PDG DataBase\n";
         l++;
         continue;
       }
     cout<<"pdg1->GetName()\n";
     p1name= pdg1->GetName();
     cout<<"pdg1->GetName() Done\n";
     k=l;
     while (PID[k] != 0)
      {
        if ( (k==l)&&(ident==kFALSE) )  { k++;continue;}
        if ( (k!=l)&&(nonident==kFALSE) ) { k++;continue; }
        

        Int_t pd1 = PID[l];
        Int_t pd2 = PID[k];
        
        Bool_t w =  weights;//weight code does not accpept pi- kaon combinations
        if ((pd1 == kPiMinus) && (TMath::Abs(pd2) == kKPlus)) w = kFALSE; 
        if ((pd1 == kPiMinus) && (pd2 == kPiMinus)) w = kFALSE; 
        if ((pd1 == kKMinus)  && (pd2 == kKMinus))  w = kFALSE; 
        
        
        cout<<"TDatabasePDG::Instance()->GetParticle(PID[k])\n";
        pdg2 = TDatabasePDG::Instance()->GetParticle(PID[k]);
        cout<<"TDatabasePDG::Instance()->GetParticle(PID[k]) Done\n";
        if (pdg2 == 0x0)
         {
           cerr<<"ERROR: hbtanalysis.C: Can not find particle type "<<PID[k]<<" in PDG DataBase\n";
           k++;
           continue;
         }
         
        cout<<"pdg2->GetName()\n";
        p2name = pdg2->GetName();
        cout<<"pdg2->GetName() Done\n";
        
        cout<<"hbtanalysis.C: l = "<<l<<" k = "<<k<<endl;
        cout<<"hbtanalysis.C: buildanalysis for "<<PID[l]<<"("<<p1name<<")  "<<PID[k]<<"("<<p2name<<")\n";
        fn = filename + p1name + p2name + ".root";
        cout<<"hbtanalysis.C: File Name is: "<<fn<<endl;
        TFile* histoOutput = TFile::Open(fn,"recreate");

        AliRunAnalysis* analysis = new AliRunAnalysis();
        cout<<"hbtanalysis.C: Set reader\n";
        analysis->SetReader(reader);

        cout<<"hbtanalysis.C: buildanalysis "<<endl;
        hbtanalysis = buildanalysis(pd1,pd2,.2,kFALSE,w,"","");
        cout<<"hbtanalysis.C: buildanalysis Done"<<endl;
        hbtanalysis->SetProcessOption(AliHBTAnalysis::kSimulatedAndReconstructed);
        hbtanalysis->SetDisplayInfo(1000);

        if (hbtanalysis == 0x0)
         {
           cout<<"ERROR: hbtanalysis.C: buildanalysis function returned NULL pointer\n";
           cout<<"ERROR: Deleting Run Analysis ..."<<endl;
           delete analysis;
           cout<<"ERROR: Deleting Run Analysis Done"<<endl;
           continue; 
         }

        analysis->Add(hbtanalysis);
        
         
 
        if (w)
         {
           AliHBTLLWeights::Instance()->SetParticlesTypes(pd1,pd2);
           
           if (pd1 == pd2)
            {
              AliHBTLLWeights::Instance()->SetColoumb(kFALSE);
              AliHBTLLWeights::Instance()->SetStrongInterSwitch(kFALSE);
              AliHBTLLWeights::Instance()->SetQuantumStatistics(kTRUE);
            }
           else
            {
              AliHBTLLWeights::Instance()->SetColoumb(kTRUE);
              AliHBTLLWeights::Instance()->SetStrongInterSwitch(kTRUE);
              AliHBTLLWeights::Instance()->SetQuantumStatistics(kFALSE);
            } 
           
           AliHBTLLWeights::Instance()->SetRandomPosition(kFALSE);
//           AliHBTLLWeights::Instance()->SetR1dw(8.0);
           
           AliHBTLLWeights::Instance()->Init();
           AliHBTLLWeights::Instance()->Set();
         }

        
        histoOutput->cd();
        cout<<"hbtanalysis.C: Run\n";
        analysis->Run();
        cout<<"hbtanalysis.C: Run Done\n";
        
        cout<<"hbtanalysis.C: Write Reader ... \n";
        histoOutput->cd();
        reader->Write();
        reader->WriteTrackCounter();
        cout<<"hbtanalysis.C: Write Reader Done \n";
        
        cout<<"hbtanalysis.C: Delete analysis ...\n";
        analysis->SetOwner(kTRUE);
        delete analysis;
        analysis = 0x0;
        cout<<"hbtanalysis.C:                     DONE\n";


        cout<<"hbtanalysis.C: Closing File ...\n";
        histoOutput->Close();
        delete histoOutput;
        cout<<"hbtanalysis.C: Closing File Done\n";
        
        if (run != 0)
         {
           TString cmd("rfcp ");
           cmd+=fn;
           cmd+=" /castor/cern.ch/user/c/chajecki/results/";
           cmd+=run;
           cmd+="/";
           cout<<"Executing command"<< cmd <<endl;
           Int_t cmdret = gSystem->Exec(cmd);
           cout<<"Command completed. Returned "<<cmdret<<endl;
         }
        k++;
      }
     l++;
   }
  
  cout<<"hbtanalysis.C: Delete readerpartcut ...\n";
  delete readerpartcut;
  cout<<"hbtanalysis.C: Delete readerpartcut ... Done\n";
  
  if (dirs) 
   {
    cout<<"hbtanalysis.C: Delete dirs ...\n";
    dirs->SetOwner();
    delete dirs;
    cout<<"hbtanalysis.C: Delete dirs ...Done\n";
   }
   
  cout<<"hbtanalysis.C: This is the very end!!!!! \n";
  ::Info("hbtanalysis.C"," This is the very end!!!!!");
//    ::Fatal("hbtanalysis.C","Maybe works???");
 }

AliHBTAnalysis* buildanalysis(Int_t pid1, Int_t pid2, Float_t qmax, Bool_t threeDcuts, Bool_t weights, const char* anasuffix, const char* titlesfx)
 {
  /************************************************************/
  /****    Build Analysises            ************************/
  /************************************************************/
  
  //Ranges of the cuts on Ous-Side-Long Functions are made to extract the best 3D signal,
  
  
  Int_t binstoscale = 200;
  Int_t nbins = 1000;
  Float_t qmin = 0.0;

  Float_t ptnbins = 20;
  Float_t ptmax = 0.2;
  
  Float_t thetanbins = 60;
  Float_t thetamax = 0.06;
  
  Float_t phinbins = 60;
  Float_t phimax = 3.15;
  
  
  AliHBTAnalysis * analysis = new AliHBTAnalysis();
  TString anasfx(anasuffix);
  TString tsfx(titlesfx);
  analysis->SetName("hbtan"+anasfx);
  analysis->SetTitle(tsfx);
  analysis->SetBufferSize(2);
//  analysis->SetCutsOnParticles();
  analysis->SetCutsOnRec();
  AliAODPairCut *paircut = new AliAODPairCut();
  
  paircut->SetQInvRange(qmin,1.5*qmax);
//  paircut->SetAvSeparationRange(2.); //AntiMerging
//  paircut->SetClusterOverlapRange(-.5,.6);//AntiSplitting

//  AliHBTSeparationCut* entrsepar = new AliHBTSeparationCut(5, 10e5, 0);
//  paircut->AddBasePairCut(entrsepar);
  
//  paircut->SetITSSeparation(0,0.025,0.08);//rphi 250, z 800 microns
//  paircut->SetITSSeparation(1,0.025,0.08);//rphi 250, z 800 microns
//  paircut->SetITSSeparation(2,0.025,0.08);//rphi 250, z 800 microns
//  paircut->SetITSSeparation(3,0.025,0.08);//rphi 250, z 800 microns
//  paircut->SetITSSeparation(4,0.025,0.1);//rphi 250, z 800 microns
//  paircut->SetITSSeparation(5,0.025,0.1);//rphi 250, z 800 microns
 
//  AliHBTDeltaPhiCut* deltaphicut = new AliHBTDeltaPhiCut(0.04,10e5);
//  AliHBTDeltaThetaCut* deltathetacut = new AliHBTDeltaThetaCut(0.04,10e5);
//  AliHBTOrPairCut* thetaphipaircut = new AliHBTOrPairCut(deltaphicut,deltathetacut);
//  paircut->AddBasePairCut(thetaphipaircut);
  
  
   
  TParticlePDG* pdg1 = TDatabasePDG::Instance()->GetParticle(pid1);
  TString system = pdg1->GetName();
  system+=" ";
  if (pid2)
   {
    cout<<"Setting cuts on particles: "<<pid1<<"  "<<pid2<<endl;
    AliAODParticleCut* partcut = new AliAODParticleCut();
    partcut->SetPID(pid1);
    paircut->SetFirstPartCut(partcut);
    partcut->SetPID(pid2);
    paircut->SetSecondPartCut(partcut);
    delete partcut;
    TParticlePDG* pdg2 = TDatabasePDG::Instance()->GetParticle(pid2);
    system+=pdg2->GetName();
   }
  else
   {
    system+=pdg1->GetName();
   }
  analysis->SetGlobalPairCut(paircut);
//  analysis->SetAntiMergingCut(2.);//in cm
  
  paircut->SetName("AnalysisGlobalPairCut");
  paircut->Write();

  AliAODPairCut *theorFpaircut = new AliAODPairCut();
  AliAODParticleCut* thFpartcut = new AliAODParticleCut();
  thFpartcut->SetPID(pid1);
  theorFpaircut->SetFirstPartCut(thFpartcut);
  thFpartcut->SetPID(pid2);
  theorFpaircut->SetSecondPartCut(thFpartcut);
  delete thFpartcut;


  AliAODPairCut* PosSidePairCut = new AliAODPairCut();
  PosSidePairCut->SetQSideLCMSRange(0.,1e5);

  AliAODPairCut* NegSidePairCut = new AliAODPairCut();
  NegSidePairCut->SetQSideLCMSRange(-1e5.,0.);

  AliAODPairCut* PosOutPairCut = new AliAODPairCut();
  PosOutPairCut->SetQOutLCMSRange(0.,1e5);

  AliAODPairCut* NegOutPairCut = new AliAODPairCut();
  NegOutPairCut->SetQOutLCMSRange(-1e5.,0.);
  
  AliAODPairCut* PosLongPairCut = new AliAODPairCut();
  PosLongPairCut->SetQLongLCMSRange(0.,1e5);

  AliAODPairCut* NegLongPairCut = new AliAODPairCut();
  NegLongPairCut->SetQLongLCMSRange(-1e5.,0.);
  

  AliAODPairCut* paircutpt1 = new AliAODPairCut();
  paircutpt1->SetKtRange(0.1,0.5);
  paircutpt1->SetQSideLCMSRange(0.0,1e5);
//  AliAODParticleCut* partcut1 = new AliAODParticleCut();
//  partcut1->SetPtRange(0.1,0.3);
//  paircutpt1->SetPartCut(partcut1);
  AliAODPairCut* paircutpt2 = new AliAODPairCut();
  paircutpt2->SetKtRange(0.5,1.0);
  paircutpt2->SetQSideLCMSRange(0.0,1e5);
//  AliAODParticleCut* partcut2 = new AliAODParticleCut();
//  partcut2->SetPtRange(0.3,0.6);
//  paircutpt2->SetPartCut(partcut2);
  

  AliAODPairCut* paircutpt3 = new AliAODPairCut();
  paircutpt3->SetKtRange(1.0,10e6);
  paircutpt3->SetQSideLCMSRange(0.0,1e5);
//  AliAODParticleCut* partcut3 = new AliAODParticleCut();
//  partcut3->SetPtRange(0.6,10e6);
//  paircutpt3->SetPartCut(partcut3);
  
  
 
  
  /************************************************************/
  /****   Q INV Correlation Function   ************************/
  /************************************************************/
  
  AliHBTQInvCorrelFctn * qinvcfP = new AliHBTQInvCorrelFctn(nbins,qmax);
  qinvcfP->SetNumberOfBinsToScale(binstoscale);
  qinvcfP->Rename("qinvcfP"+anasfx,"Particle (simulated) Q_{inv} Correlation Function "+system);
  analysis->AddParticleFunction(qinvcfP);

  AliHBTQInvCorrelFctn * qinvcfT = new AliHBTQInvCorrelFctn(nbins,qmax);
  qinvcfT->SetNumberOfBinsToScale(binstoscale);
  qinvcfT->Rename("qinvcfT"+anasfx,"Track (recontructed) Q_{inv} Correlation Function "+system);
  analysis->AddTrackFunction(qinvcfT);

  AliHBTTwoKStarCorrelFctn* kstarP = new AliHBTTwoKStarCorrelFctn(nbins,qmax);
  kstarP->SetNumberOfBinsToScale(binstoscale);
  kstarP->Rename("kstarP"+anasfx,"Particle (simulated) 2K^{*} Correlation Function "+system);
  analysis->AddParticleFunction(kstarP);

  AliHBTTwoKStarCorrelFctn* kstarT = new AliHBTTwoKStarCorrelFctn(nbins,qmax);
  kstarT->SetNumberOfBinsToScale(binstoscale);
  kstarT->Rename("kstarT"+anasfx,"Track (recontructed) 2K^{*} Correlation Function "+system);
  analysis->AddTrackFunction(kstarT);


  

  //LEDNICKY
  if (weights)
   {
   
     AliHBTWeightTheorQInvFctn  *wqinvcfP = new AliHBTWeightTheorQInvFctn(nbins,qmax);
     wqinvcfP->SetNumberOfBinsToScale(binstoscale);
     wqinvcfP->Rename("wqinvcfP"+anasfx,"Lednicky Q_{inv} Theoretical Correlation Function "+system);
     wqinvcfP->SetPairCut(theorFpaircut);
     analysis->AddParticleFunction(wqinvcfP);
     
     AliHBTWeightQInvFctn  *wqinvcf = new AliHBTWeightQInvFctn(nbins,qmax);
     wqinvcf->SetNumberOfBinsToScale(binstoscale);
     wqinvcf->Rename("wqinvcf"+anasfx,"Lednicky Q_{inv} Correlation Function "+system);
     wqinvcf->SetPairCut(theorFpaircut);
     analysis->AddParticleAndTrackFunction(wqinvcf);
     
     
     AliHBTWeightTwoKStarFctn* wkstarP = new AliHBTWeightTwoKStarFctn(nbins,qmax);
     wkstarP->SetNumberOfBinsToScale(binstoscale);
     wkstarP->Rename("wkstar"+anasfx,"Lednicky 2K^{*} Correlation Function "+system);
     analysis->AddParticleAndTrackFunction(wkstarP);


     AliHBTWeightTwoKStarFctn* wkstarNegS = new AliHBTWeightTwoKStarFctn(nbins,qmax);
     wkstarNegS->SetNumberOfBinsToScale(binstoscale);
     wkstarNegS->Rename("wkstarNegS"+anasfx,"Lednicky 2K^{*} Correlation Function Neg Side"+system);
     wkstarNegS->SetPairCut(NegSidePairCut);
     analysis->AddParticleAndTrackFunction(wkstarNegS);
     
     AliHBTWeightTwoKStarFctn* wkstarPosS = new AliHBTWeightTwoKStarFctn(nbins,qmax);
     wkstarPosS->SetNumberOfBinsToScale(binstoscale);
     wkstarPosS->Rename("wkstarPosS"+anasfx,"Lednicky 2K^{*} Correlation Function Neg Side "+system);
     wkstarPosS->SetPairCut(PosSidePairCut);
     analysis->AddParticleAndTrackFunction(wkstarPosS);
     

     AliHBTWeightTwoKStarFctn* wkstarNegO = new AliHBTWeightTwoKStarFctn(nbins,qmax);
     wkstarNegO->SetNumberOfBinsToScale(binstoscale);
     wkstarNegO->Rename("wkstarNegO"+anasfx,"Lednicky 2K^{*} Correlation Function Neg Out"+system);
     wkstarNegO->SetPairCut(NegOutPairCut);
     analysis->AddParticleAndTrackFunction(wkstarNegO);
     
     AliHBTWeightTwoKStarFctn* wkstarPosO = new AliHBTWeightTwoKStarFctn(nbins,qmax);
     wkstarPosO->SetNumberOfBinsToScale(binstoscale);
     wkstarPosO->Rename("wkstarPosO"+anasfx,"Lednicky 2K^{*} Correlation Function Pos Out"+system);
     wkstarPosO->SetPairCut(PosOutPairCut);
     analysis->AddParticleAndTrackFunction(wkstarPosO);

     AliHBTWeightTwoKStarFctn* wkstarNegL = new AliHBTWeightTwoKStarFctn(nbins,qmax);
     wkstarNegL->SetNumberOfBinsToScale(binstoscale);
     wkstarNegL->Rename("wkstarNegL"+anasfx,"Lednicky 2K^{*} Correlation Function Neg Long"+system);
     wkstarNegL->SetPairCut(NegLongPairCut);
     analysis->AddParticleAndTrackFunction(wkstarNegL);
     
     AliHBTWeightTwoKStarFctn* wkstarPosL = new AliHBTWeightTwoKStarFctn(nbins,qmax);
     wkstarPosL->SetNumberOfBinsToScale(binstoscale);
     wkstarPosL->Rename("wkstarPosL"+anasfx,"Lednicky 2K^{*} Correlation Function Pos Long "+system);
     wkstarPosL->SetPairCut(PosLongPairCut);
     analysis->AddParticleAndTrackFunction(wkstarPosL);



     
   }


  /************************************************************/
  /****   Q OUT Correlation Function   ************************/
  /************************************************************/
  
  AliHBTQOutLCMSCorrelFctn* qoutP = new AliHBTQOutLCMSCorrelFctn(nbins,qmax,-qmax);
  qoutP->Rename("qoutP"+anasfx,"Particle (simulated) Q_{out} CF "+system);
  qoutP->SetNumberOfBinsToScale(binstoscale);
  
  AliHBTQOutLCMSCorrelFctn* qoutT = new AliHBTQOutLCMSCorrelFctn(nbins,qmax,-qmax);
  qoutT->Rename("qoutT"+anasfx,"Track (recontructed) Q_{out} CF "+system);
  qoutT->SetNumberOfBinsToScale(binstoscale);

  AliAODPairCut *outPairCut = new AliAODPairCut();

  outPairCut->SetQSideLCMSRange(-0.015,0.015);
  outPairCut->SetQLongLCMSRange(-0.015,0.015);
  
  qoutP->SetPairCut(outPairCut);
  qoutT->SetPairCut(outPairCut);
  
  analysis->AddParticleFunction(qoutP);
  analysis->AddTrackFunction(qoutT);
  
  if (weights)
   {
     AliHBTWeightQOutFctn* wqoutP = new AliHBTWeightQOutFctn(nbins,qmax,-qmax);
     wqoutP->SetNumberOfBinsToScale(binstoscale);
     wqoutP->Rename("wqout"+anasfx,"Lednicky Q_{out} CF "+system);
     wqoutP->SetPairCut(outPairCut);
     analysis->AddParticleAndTrackFunction(wqoutP);
   }
  /************************************************************/
  /****   Q SIDE Correlation Function   ***********************/
  /************************************************************/
  
  AliHBTQSideLCMSCorrelFctn* qsideP = new AliHBTQSideLCMSCorrelFctn(nbins,qmax,-qmax);
  qsideP->Rename("qsideP"+anasfx,"Particle (simulated) Q_{side} CF "+system);
  qsideP->SetNumberOfBinsToScale(binstoscale);
  
  AliHBTQSideLCMSCorrelFctn* qsideT = new AliHBTQSideLCMSCorrelFctn(nbins,qmax,-qmax);
  qsideT->Rename("qsideT"+anasfx,"Track (recontructed) Q_{side} CF "+system);
  qsideT->SetNumberOfBinsToScale(binstoscale);

  AliAODPairCut *sidePairCut = new AliAODPairCut();
  sidePairCut->SetQOutLCMSRange(-0.02,0.02);
  sidePairCut->SetQLongLCMSRange(-0.015,0.015);
  
  qsideP->SetPairCut(sidePairCut);
  qsideT->SetPairCut(sidePairCut);
  
  analysis->AddParticleFunction(qsideP);
  analysis->AddTrackFunction(qsideT);
  
  if (weights)
   {
     AliHBTWeightQSideFctn* wqside = new AliHBTWeightQSideFctn(nbins,qmax,-qmax);
     wqside->SetNumberOfBinsToScale(binstoscale);
     wqside->Rename("wqside"+anasfx,"Lednicky Q_{side} CF "+system);
     wqside->SetPairCut(sidePairCut);
     analysis->AddParticleAndTrackFunction(wqside);
   }
  /************************************************************/
  /****   Q LONG Correlation Function   ***********************/
  /************************************************************/
    
  AliHBTQLongLCMSCorrelFctn* qlongP = new AliHBTQLongLCMSCorrelFctn(nbins,qmax,-qmax);
  qlongP->Rename("qlongP"+anasfx,"Particle (simulated) Q_{long} CF "+system);
  qlongP->SetNumberOfBinsToScale(binstoscale);

  AliHBTQLongLCMSCorrelFctn* qlongT = new AliHBTQLongLCMSCorrelFctn(nbins,qmax,-qmax);
  qlongT->Rename("qlongT"+anasfx,"Track (recontructed) Q_{long} CF "+system);
  qlongT->SetNumberOfBinsToScale(binstoscale);

  AliAODPairCut *longPairCut = new AliAODPairCut();
  longPairCut->SetQOutLCMSRange(-0.02,0.02);
  longPairCut->SetQSideLCMSRange(-0.015,0.015);
  longPairCut->SetQLongLCMSRange(-qmax,qmax);
  
  qlongP->SetPairCut(longPairCut);
  qlongT->SetPairCut(longPairCut);
  
  analysis->AddParticleFunction(qlongP);
  analysis->AddTrackFunction(qlongT);
  
  if (weights)
   {
     AliHBTWeightQLongFctn* wqlong = new AliHBTWeightQLongFctn(nbins,qmax,-qmax);
     wqlong->SetNumberOfBinsToScale(binstoscale);
     wqlong->Rename("wqlong"+anasfx,"Lednicky Q_{long} CF "+system);
     wqlong->SetPairCut(longPairCut);
     analysis->AddParticleAndTrackFunction(wqlong);
     
   }
   
  /************************************************************/
  /****   Q OUT - Q SIDE Correlation Function   ***************/
  /************************************************************/
   
  AliAODPairCut *outsidePairCut = new AliAODPairCut();
  outsidePairCut->SetQLongLCMSRange(-0.02,0.02);
  
  if (weights)
   {
     AliHBTWeightQOutQSideFctn* wqoutqside = new AliHBTWeightQOutQSideFctn(200,qmax,-qmax,200,qmax,-qmax);
     wqoutqside->SetNumberOfBinsToScale(40,40);
     wqoutqside->Rename("wqoutqside"+anasfx,"Lednicky Weight Q_{out}-Q_{side} CF "+system);
     wqoutqside->SetPairCut(outsidePairCut);
     analysis->AddParticleAndTrackFunction(wqoutqside);
     
   }
   

  /************************************************************/
  /****   Q OUT - Q LONG Correlation Function   ***************/
  /************************************************************/
   
  AliAODPairCut *outlongPairCut = new AliAODPairCut();
  outlongPairCut->SetQSideLCMSRange(-0.02,0.02);
  
  if (weights)
   {
     AliHBTWeightQOutQLongFctn* wqoutqlong = new AliHBTWeightQOutQLongFctn(200,qmax,-qmax,200,qmax,-qmax);
     wqoutqlong->SetNumberOfBinsToScale(40,40);
     wqoutqlong->Rename("wqoutqlong"+anasfx,"Lednicky Weight Q_{out}-Q_{long} CF "+system);
     wqoutqlong->SetPairCut(outlongPairCut);
     analysis->AddParticleAndTrackFunction(wqoutqlong);
     
   }

  /************************************************************/
  /****   Q SIDE - Q LONG Correlation Function   ***************/
  /************************************************************/
   
  AliAODPairCut *sidelongPairCut = new AliAODPairCut();
  sidelongPairCut->SetQOutLCMSRange(-0.02,0.02);
  
  if (weights)
   {
     AliHBTWeightQSideQLongFctn* wqsideqlong  = new AliHBTWeightQSideQLongFctn(200,qmax,-qmax,200,qmax,-qmax);
     wqsideqlong->SetNumberOfBinsToScale(40,40);
     wqsideqlong->Rename("wqsideqlong"+anasfx,"Lednicky Weight Q_{side}-Q_{long} CF "+system);
     wqsideqlong->SetPairCut(sidelongPairCut);
     analysis->AddParticleAndTrackFunction(wqsideqlong);
     
   }
  
  /************************************************************/
  /************************************************************/
  /************************************************************/
  /****    M O N I T O R  R E S O L U T I O N S   *************/
  /************************************************************/
  /************************************************************/
  /************************************************************/

  AliHBTMonPxResolutionVsPtFctn* PxVsPtResF = new AliHBTMonPxResolutionVsPtFctn(200,2.0,0.0,200,0.06,-0.06);
  analysis->AddParticleAndTrackMonitorFunction(PxVsPtResF);

  AliHBTMonPyResolutionVsPtFctn* PyVsPtResF = new AliHBTMonPyResolutionVsPtFctn(200,2.0,0.0,200,0.06,-0.06);
  analysis->AddParticleAndTrackMonitorFunction(PyVsPtResF);

  AliHBTMonPzResolutionVsPtFctn* PzVsPtResF = new AliHBTMonPzResolutionVsPtFctn(200,2.0,0.0,200,0.06,-0.06);
  analysis->AddParticleAndTrackMonitorFunction(PzVsPtResF);

  AliHBTMonPResolutionVsPtFctn* PVsPtResF = new AliHBTMonPResolutionVsPtFctn(200,2.0,0.0,200,0.06,-0.06);
  analysis->AddParticleAndTrackMonitorFunction(PVsPtResF);

  AliHBTMonPtResolutionVsPtFctn* PtVsPtResF = new AliHBTMonPtResolutionVsPtFctn(200,2.0,0.0,200,0.06,-0.06);
  analysis->AddParticleAndTrackMonitorFunction(PtVsPtResF);

  AliHBTMonPhiResolutionVsPtFctn* PhiVsPtResF = new AliHBTMonPhiResolutionVsPtFctn(200,2.0,0.0,100,0.04,-0.04);
  analysis->AddParticleAndTrackMonitorFunction(PhiVsPtResF);

  AliHBTMonThetaResolutionVsPtFctn* ThetaVsPtResF = new AliHBTMonThetaResolutionVsPtFctn(200,2.0,0.0,100,0.04,-0.04);
  analysis->AddParticleAndTrackMonitorFunction(ThetaVsPtResF);

  AliHBTMonPhiResolutionVsPhiFctn* PhiVsPhiResF = 
       new AliHBTMonPhiResolutionVsPhiFctn(200,TMath::TwoPi(),0.0,200,0.04,-0.04);
  analysis->AddParticleAndTrackMonitorFunction(PhiVsPhiResF);

  AliHBTMonThetaResolutionVsThetaFctn* ThetaVsThetaResF = 
        new AliHBTMonThetaResolutionVsThetaFctn(200,0.75*TMath::Pi(),0.25*TMath::Pi(),200,0.04,-0.04);
  analysis->AddParticleAndTrackMonitorFunction(ThetaVsThetaResF);


  AliHBTMonPxDistributionFctn* pxdistr = new AliHBTMonPxDistributionFctn();
  AliHBTMonPyDistributionFctn* pydistr = new AliHBTMonPyDistributionFctn();
  AliHBTMonPzDistributionFctn* pzdistr = new AliHBTMonPzDistributionFctn();
  AliHBTMonPDistributionFctn*  pdistr  = new AliHBTMonPDistributionFctn();
  AliHBTMonPtDistributionFctn* ptdistr = new AliHBTMonPtDistributionFctn();
  
  analysis->AddTrackMonitorFunction(pxdistr);
  analysis->AddTrackMonitorFunction(pydistr);
  analysis->AddTrackMonitorFunction(pzdistr);
  analysis->AddTrackMonitorFunction(pdistr);
  analysis->AddTrackMonitorFunction(ptdistr);

  AliHBTMonPhiDistributionVsPtFctn*  phidistr  = new AliHBTMonPhiDistributionVsPtFctn();
  AliHBTMonThetaDistributionVsPtFctn* thetadistr = new AliHBTMonThetaDistributionVsPtFctn();
  analysis->AddParticleMonitorFunction(phidistr);
  analysis->AddParticleMonitorFunction(thetadistr);
  
  AliHBTMonVxDistributionFctn* vxdistr = new AliHBTMonVxDistributionFctn(200,1.0e-1,-1.0e-1);
  AliHBTMonVyDistributionFctn* vydistr = new AliHBTMonVyDistributionFctn(200,1.0e-1,-1.0e-1);
  AliHBTMonVzDistributionFctn* vzdistr = new AliHBTMonVzDistributionFctn(200,1.0e-1,-1.0e-1);
  AliHBTMonRDistributionFctn*  vrdistr = new AliHBTMonRDistributionFctn (200,1.0e-1,-1.0e-1);
  analysis->AddTrackMonitorFunction(vxdistr);
  analysis->AddTrackMonitorFunction(vydistr);
  analysis->AddTrackMonitorFunction(vzdistr);
  analysis->AddTrackMonitorFunction(vrdistr);



  /************************************************************/
  /************************************************************/
  /************************************************************/
  /****   D I S T R B U T I O N S         ***********************/
  /************************************************************/
  /************************************************************/
  /************************************************************/

  cout<<"Check 7"<<endl;


  /************************************************************/
  /****   OUT VS KT                 ***********************/
  /************************************************************/
  AliHBTQOutDistributionVsKtFctn* qoutVsKtDistP = new AliHBTQOutDistributionVsKtFctn(200,4.0,0.0,100,qmax,-qmax);
  qoutVsKtDistP->Rename("qOutVsKtDistP",qoutVsKtDistP->GetTitle());
  AliHBTQOutDistributionVsKtFctn* qoutVsKtDistT = new AliHBTQOutDistributionVsKtFctn(200,4.0,0.0,100,qmax,-qmax);
  qoutVsKtDistT->Rename("qOutVsKtDistT",qoutVsKtDistT->GetTitle());
  if (threeDcuts)  qoutVsKtDistP->SetPairCut(outPairCut);
  if (threeDcuts)  qoutVsKtDistT->SetPairCut(outPairCut);
  analysis->AddParticleFunction(qoutVsKtDistP);
  analysis->AddTrackFunction(qoutVsKtDistT);
  /************************************************************/
  /****   Side VS KT                 ***********************/
  /************************************************************/
  AliHBTQSideDistributionVsKtFctn* qSideVsKtDistP = new AliHBTQSideDistributionVsKtFctn(200,4.0,0.0,100,qmax,-qmax);
  qSideVsKtDistP->Rename("qSideVsKtDistP",qSideVsKtDistP->GetTitle());
  AliHBTQSideDistributionVsKtFctn* qSideVsKtDistT = new AliHBTQSideDistributionVsKtFctn(200,4.0,0.0,100,qmax,-qmax);
  qSideVsKtDistT->Rename("qSideVsKtDistT",qSideVsKtDistT->GetTitle());
  if (threeDcuts)  qSideVsKtDistP->SetPairCut(sidePairCut);
  if (threeDcuts)  qSideVsKtDistT->SetPairCut(sidePairCut);
  analysis->AddParticleFunction(qSideVsKtDistP);
  analysis->AddTrackFunction(qSideVsKtDistT);
  /************************************************************/
  /****   Long VS KT                 ***********************/
  /************************************************************/
  AliHBTQLongDistributionVsKtFctn* qLongVsKtDistP = new AliHBTQLongDistributionVsKtFctn(200,4.0,0.0,100,qmax,-qmax);
  qLongVsKtDistP->Rename("qLongVsKtDistP",qLongVsKtDistP->GetTitle());
  AliHBTQLongDistributionVsKtFctn* qLongVsKtDistT = new AliHBTQLongDistributionVsKtFctn(200,4.0,0.0,100,qmax,-qmax);
  qLongVsKtDistT->Rename("qLongVsKtDistT",qLongVsKtDistT->GetTitle());
  if (threeDcuts)  qLongVsKtDistP->SetPairCut(longPairCut);
  if (threeDcuts)  qLongVsKtDistT->SetPairCut(longPairCut);
  analysis->AddParticleFunction(qLongVsKtDistP);
  analysis->AddTrackFunction(qLongVsKtDistT);


  /************************************************************/
  /****   OUT VS QINV                 ***********************/
  /************************************************************/
  AliHBTQOutDistributionVsQInvFctn* qoutVsQInvDistP = new AliHBTQOutDistributionVsQInvFctn(100,qmax,qmin,100,qmax,-qmax);
  qoutVsQInvDistP->Rename("qOutVsQInvDistP",qoutVsQInvDistP->GetTitle());
  AliHBTQOutDistributionVsQInvFctn* qoutVsQInvDistT = new AliHBTQOutDistributionVsQInvFctn(100,qmax,qmin,100,qmax,-qmax);
  qoutVsQInvDistT->Rename("qOutVsQInvDistT",qoutVsQInvDistT->GetTitle());
  if (threeDcuts)  qoutVsQInvDistP->SetPairCut(outPairCut);  if (threeDcuts)  qoutVsQInvDistT->SetPairCut(outPairCut);
  analysis->AddParticleFunction(qoutVsQInvDistP);
  analysis->AddTrackFunction(qoutVsQInvDistT);
  /************************************************************/
  /****   Side VS QINV                 ***********************/
  /************************************************************/
  AliHBTQSideDistributionVsQInvFctn* qSideVsQInvDistP = new AliHBTQSideDistributionVsQInvFctn(100,qmax,qmin,100,qmax,-qmax);
  qSideVsQInvDistP->Rename("qSideVsQInvDistP",qSideVsQInvDistP->GetTitle());
  AliHBTQSideDistributionVsQInvFctn* qSideVsQInvDistT = new AliHBTQSideDistributionVsQInvFctn(100,qmax,qmin,100,qmax,-qmax);
  qSideVsQInvDistT->Rename("qSideVsQInvDistT",qSideVsQInvDistT->GetTitle());
  if (threeDcuts)  qSideVsQInvDistP->SetPairCut(sidePairCut);
  if (threeDcuts)  qSideVsQInvDistT->SetPairCut(sidePairCut);
  analysis->AddParticleFunction(qSideVsQInvDistP);
  analysis->AddTrackFunction(qSideVsQInvDistT);
  /************************************************************/
  /****   Long VS QINV                 ***********************/
  /************************************************************/
  AliHBTQLongDistributionVsQInvFctn* qLongVsQInvDistP = new AliHBTQLongDistributionVsQInvFctn(100,qmax,qmin,100,qmax,-qmax);
  qLongVsQInvDistP->Rename("qLongVsQInvDistP",qLongVsQInvDistP->GetTitle());
  AliHBTQLongDistributionVsQInvFctn* qLongVsQInvDistT = new AliHBTQLongDistributionVsQInvFctn(100,qmax,qmin,100,qmax,-qmax);
  qLongVsQInvDistT->Rename("qLongVsQInvDistT",qLongVsQInvDistT->GetTitle());
  if (threeDcuts)  qLongVsQInvDistP->SetPairCut(longPairCut);
  if (threeDcuts)  qLongVsQInvDistT->SetPairCut(longPairCut);
  analysis->AddParticleFunction(qLongVsQInvDistP);
  analysis->AddTrackFunction(qLongVsQInvDistT);




  delete outPairCut;
  delete sidePairCut;
  delete longPairCut;

  delete paircut;
  
  delete paircutpt1;
//  delete partcut1;

  delete paircutpt2;
//  delete partcut2;

  delete paircutpt3;
//  delete partcut3;
  
  return analysis;


 }