PWG2/FLOW/other/runFlowAnalysis.C

   1 #include "TStopwatch.h"
   2 #include "TObjArray"
   3 #include "Riostream.h"
   4
   5 //RUN SETTINGS
   6 //flow analysis method can be: (set to kTRUE or kFALSE)
   7 Bool_t SP    = kTRUE;
   8 Bool_t LYZ1  = kTRUE;
   9 Bool_t LYZ2  = kFALSE;
  10 Bool_t LYZEP = kFALSE;
  11 Bool_t GFC   = kTRUE;
  12 Bool_t QC    = kTRUE;
  13 Bool_t FQD   = kTRUE;
  14 Bool_t MCEP  = kFALSE; //does not work yet 24/12/08
  15
  16 Int_t offset = 0;
  17
  18 //CUT SETTINGS
  19 //integrated selection
  20 Double_t ptMaxInt  = 10.;
  21 Double_t ptMinInt  = 0.;
  22 Double_t etaMaxInt = 1.;
  23 Double_t etaMinInt = -1.;
  24 Double_t phiMaxInt = 7.5;
  25 Double_t phiMinInt = 0.;
  26 Int_t PIDInt       = 211;
  27
  28 //differential selection
  29 Double_t ptMaxDiff  = 10.;
  30 Double_t ptMinDiff  = 0.;
  31 Double_t etaMaxDiff = 1.;
  32 Double_t etaMinDiff = -1.;
  33 Double_t phiMaxDiff = 7.5;
  34 Double_t phiMinDiff = 0.;
  35 Int_t PIDDiff       = 211;
  36
  37 int runFlowAnalysis(Int_t aRuns = 100, const char*
  38                           //                      dir="/data/alice1/kolk/KineOnly3/")
  39                           dir="/Users/snelling/alice_data/KineOnly3/")
  40 {
  41   TStopwatch timer;
  42   timer.Start();
  43
  44   if (LYZ1 && LYZ2) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(); }
  45   if (LYZ2 && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
  46   if (LYZ1 && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
  47
  48
  49   cout<<endl;
  50   cout<<" ---- BEGIN ANALYSIS ---- "<<endl;
  51   cout<<endl;
  52
  53   //  gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  54   gSystem->AddIncludePath("-I$ROOTSYS/include");
  55
  56   // constants
  57   gROOT->LoadMacro("code/AliFlowCommonConstants.cxx+");
  58   gROOT->LoadMacro("code/AliFlowLYZConstants.cxx+");
  59   gROOT->LoadMacro("code/AliFlowCumuConstants.cxx+");
  60
  61   // flow event
  62   gROOT->LoadMacro("code/AliFlowVector.cxx+");
  63   gROOT->LoadMacro("code/AliFlowTrackSimple.cxx+");
  64   gROOT->LoadMacro("code/AliFlowEventSimple.cxx+");
  65
  66   // cuts
  67   gROOT->LoadMacro("code/AliFlowTrackSimpleCuts.cxx+");
  68
  69   // output histosgrams
  70   gROOT->LoadMacro("code/AliFlowCommonHist.cxx+");
  71   gROOT->LoadMacro("code/AliFlowCommonHistResults.cxx+");
  72   gROOT->LoadMacro("code/AliFlowLYZHist1.cxx+");
  73   gROOT->LoadMacro("code/AliFlowLYZHist2.cxx+");
  74
  75   // functions needed for various methods
  76   gROOT->LoadMacro("code/AliCumulantsFunctions.cxx+");
  77   gROOT->LoadMacro("code/AliQCumulantsFunctions.cxx+");
  78   gROOT->LoadMacro("code/AliFittingFunctionsForQDistribution.cxx+");
  79   gROOT->LoadMacro("code/AliFlowLYZEventPlane.cxx+");
  80
  81   // Flow Analysis code for various methods
  82   gROOT->LoadMacro("code/AliFlowAnalysisWithMCEventPlane.cxx+");
  83   gROOT->LoadMacro("code/AliFlowAnalysisWithScalarProduct.cxx+");
  84   gROOT->LoadMacro("code/AliFlowAnalysisWithLYZEventPlane.cxx+");
  85   gROOT->LoadMacro("code/AliFlowAnalysisWithLeeYangZeros.cxx+");
  86   gROOT->LoadMacro("code/AliFlowAnalysisWithCumulants.cxx+");
  87   gROOT->LoadMacro("code/AliFlowAnalysisWithQCumulants.cxx+");
  88   gROOT->LoadMacro("code/AliFittingQDistribution.cxx+");
  89
  90   // Method to fill the FlowEvent
  91   gROOT->LoadMacro("code/FlowEventSimpleMaker.cxx+");
  92
  93   //load needed libraries
  94   gSystem->Load("libTree.so");
  95
  96   cout << "finished loading macros!" << endl;
  97   //  gSystem->Load("libPWG2flow.so");
  98
  99   //------------------------------------------------------------------------
 100   //cuts
 101   AliFlowTrackSimpleCuts* cutsInt = new AliFlowTrackSimpleCuts();
 102   cutsInt->SetPtMax(ptMaxInt);
 103   cutsInt->SetPtMin(ptMinInt);
 104   cutsInt->SetEtaMax(etaMaxInt);
 105   cutsInt->SetEtaMin(etaMinInt);
 106   cutsInt->SetPhiMax(phiMaxInt);
 107   cutsInt->SetPhiMin(phiMinInt);
 108   cutsInt->SetPID(PIDInt);
 109
 110   AliFlowTrackSimpleCuts* cutsDiff = new AliFlowTrackSimpleCuts();
 111   cutsDiff->SetPtMax(ptMaxDiff);
 112   cutsDiff->SetPtMin(ptMinDiff);
 113   cutsDiff->SetEtaMax(etaMaxDiff);
 114   cutsDiff->SetEtaMin(etaMinDiff);
 115   cutsDiff->SetPhiMax(phiMaxDiff);
 116   cutsDiff->SetPhiMin(phiMinDiff);
 117   cutsDiff->SetPID(PIDDiff);
 118
 119   //------------------------------------------------------------------------
 120   //flow event
 121   FlowEventSimpleMaker* fEventMaker = new FlowEventSimpleMaker();
 122
 123   AliFlowAnalysisWithQCumulants    *qc    = NULL;
 124   AliFlowAnalysisWithCumulants     *gfc   = NULL;
 125   AliFittingQDistribution          *fqd   = NULL;
 126   AliFlowAnalysisWithLeeYangZeros  *lyz1  = NULL;
 127   AliFlowAnalysisWithLeeYangZeros  *lyz2  = NULL;
 128   AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
 129   AliFlowAnalysisWithScalarProduct *sp    = NULL;
 130   AliFlowAnalysisWithMCEventPlane  *mcep  = NULL;
 131
 132   //flow methods:
 133
 134   //MCEP = monte carlo event plane
 135   if (MCEP) {
 136     AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
 137     mcep->Init();
 138   }
 139
 140   //QC = Q-cumulants
 141   if(QC) {
 142     AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
 143     qc->CreateOutputObjects();
 144   }
 145
 146   //GFC = Generating Function Cumulants
 147   if(GFC) {
 148     AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
 149     gfc->CreateOutputObjects();
 150   }
 151
 152   //FQD = Fitting q-distribution
 153   if(FQD) {
 154     AliFittingQDistribution* fqd = new AliFittingQDistribution();
 155     fqd->CreateOutputObjects();
 156   }
 157   //LYZ1 = Lee-Yang Zeroes first run
 158   if(LYZ1) {
 159     AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
 160     lyz1->SetFirstRun(kTRUE);
 161     lyz1->SetUseSum(kTRUE);
 162     lyz1->Init();
 163   }
 164
 165   //LYZ2 = Lee-Yang Zeroes second run
 166
 167   if(LYZ2) {
 168     AliFlowAnalysisWithLeeYangZeros* lyz2 = new AliFlowAnalysisWithLeeYangZeros();
 169     // read the input file from the first run
 170     TString inputFileNameLYZ2 = "outputLYZ1analysis.root" ;
 171     TFile* inputFileLYZ2 = new TFile(inputFileNameLYZ2.Data(),"READ");
 172     if(!inputFileLYZ2 || inputFileLYZ2->IsZombie()) {
 173       cerr << " ERROR: NO First Run file... " << endl ;
 174     }
 175     else {
 176       TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");
 177       if (!inputListLYZ2) {cout<<"list is NULL pointer!"<<endl;}
 178       else {
 179         cout<<"LYZ2 input file/list read..."<<endl;
 180         lyz2->SetFirstRunList(inputListLYZ2);
 181         lyz2->SetFirstRun(kFALSE);
 182         lyz2->SetUseSum(kTRUE);
 183         lyz2->Init();
 184       }
 185     }
 186   }
 187
 188  //LYZEP = Lee-Yang Zeroes event plane
 189   if(LYZEP) {
 190     AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
 191     AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
 192     // read the input file from the second lyz run
 193     TString inputFileNameLYZEP = "outputLYZ2analysis.root" ;
 194     TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
 195     if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
 196       cerr << " ERROR: NO Second Run file... " << endl ; }
 197     else {
 198       TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");
 199       if (!inputListLYZEP) {cout<<"list is NULL pointer!"<<endl;}
 200       else {
 201         cout<<"LYZEP input file/list read..."<<endl;
 202         ep   ->SetSecondRunList(inputListLYZEP);
 203         lyzep->SetSecondRunList(inputListLYZEP);
 204         ep   ->Init();
 205         lyzep->Init();
 206       }
 207     }
 208   }
 209
 210   //SP = Scalar Product
 211   if(SP) {
 212     AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
 213     sp->Init();
 214   }
 215   //------------------------------------------------------------------------
 216
 217   //standard code
 218   Int_t fCount = 0;
 219   TString execDir(gSystem->pwd());
 220   TString targetDir(dir);
 221   TSystemDirectory* baseDir = new TSystemDirectory(".", dir);
 222   TList* dirList           = baseDir->GetListOfFiles();
 223   if (!dirList) {
 224     cout << endl << "No input files in: " << targetDir.Data() << endl;
 225     break;
 226   }
 227   Int_t nDirs              = dirList->GetEntries();
 228   cout<<endl;
 229   cout<<"Int_t nDirs = "<<nDirs<<endl;
 230   gSystem->cd(execDir);
 231
 232   for(Int_t iDir=0;iDir<nDirs;++iDir)
 233     {
 234       TSystemFile* presentDir = (TSystemFile*)dirList->At(iDir);
 235       if(!presentDir || !presentDir->IsDirectory() ||
 236          strcmp(presentDir->GetName(), ".") == 0 ||
 237          strcmp(presentDir->GetName(), "..") == 0)
 238         {
 239           cout << endl;
 240           cout << "Directory (" << iDir << "): " << presentDir->GetName() <<
 241             " - Skipping ... " << endl;
 242           continue ;
 243         }
 244
 245       if(offset > 0) { --offset ; continue ; }
 246       if((aRuns > 0) && (fCount >= aRuns)) { break ; }
 247
 248       TString presentDirName(dir); // aDataDir
 249       presentDirName += presentDir->GetName();
 250       presentDirName += "/";
 251       //cerr<<" presentDirName = "<<presentDirName<<endl;
 252
 253       TString fileName = presentDirName;
 254       fileName += "galice.root";
 255       Long_t *id, *size, *flags, *modtime;
 256       if(gSystem->GetPathInfo(fileName.Data(),id,size,flags,modtime))
 257         {
 258           cout << " File : " << fileName << " does NOT exist ! - Skipping ... "
 259                << endl;
 260           continue;
 261         }
 262       //      cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName << " ... " << endl;
 263
 264       //loop (simulations in the present dir)
 265       TSystemDirectory* evtsDir = new TSystemDirectory(".",
 266                                                        presentDirName.Data());
 267       TList* fileList       = evtsDir->GetListOfFiles();
 268       Int_t nFiles                  = fileList->GetEntries();
 269       //cout<<" Int_t nFiles = "<<nFiles<<endl;
 270       gSystem->cd(execDir);
 271       for(Int_t iFiles=0; iFiles<nFiles; ++iFiles)
 272         {
 273           TSystemFile* presentFile = (TSystemFile*) fileList->At(iFiles);
 274           TString presentFileName(presentDirName);
 275           presentFileName += presentFile->GetName();
 276
 277           if(!(presentFileName.Contains("Kinematics") &&
 278                presentFileName.Contains("root"))) { continue ; }
 279
 280           //cout << " found: " << presentFileName.Data() << endl;
 281
 282           TFile* kineFile = new TFile(presentFileName.Data(), "READ");
 283           // kineFile->ls();
 284           Int_t nEvts = kineFile->GetNkeys() ;
 285           //cout << "  . found: " << nEvts << " KineTree(s) in " << presentFileName.Data() << endl;
 286           TList* kineEventsList = (TList*)kineFile->GetListOfKeys();
 287           TTree* kTree;
 288           TIter next(kineEventsList);
 289           TKey* key;
 290
 291           //loop over the events
 292           while( key=(TKey *)next() )
 293             {
 294               TDirectory* tDir = (TDirectory*)key->ReadObj();
 295               if(!tDir) break;
 296
 297               TString evtDir(tDir->GetName());
 298               //cout << "  . . found: " << tDir->GetName() << endl;
 299
 300               kTree = (TTree *)tDir->Get("TreeK");
 301               if(!kTree) break;
 302
 303               Int_t nPart = kTree->GetEntries();
 304               //cout << "  . . . kTree " << fCount << " has " << nPart << " particles " << endl;
 305
 306               //-----------------------------------------------------------
 307               //fill and save the flow event
 308
 309               AliFlowEventSimple* fEvent = fEventMaker->FillTracks(kTree, cutsInt, cutsDiff);
 310
 311               //pass the flow event to flow methods for analysis
 312               if (MCEP) //mcep->Make(fEvent,fEP);  //fix fEP
 313               if(QC) qc->Make(fEvent);
 314               if(GFC) gfc->Make(fEvent);
 315               if(FQD) fqd->Make(fEvent);
 316               if(LYZ1) lyz1->Make(fEvent);
 317               if(LYZ2) lyz2->Make(fEvent);
 318               if(LYZEP) lyzep->Make(fEvent,ep);
 319               if(SP) sp->Make(fEvent);
 320               //-----------------------------------------------------------
 321
 322               fCount++;
 323               cout << "# " << fCount << " events processed" << endl;
 324               delete kTree;
 325               delete fEvent;
 326             }
 327           delete kineFile ;
 328         }
 329       delete evtsDir ;
 330     }
 331
 332   //--------------------------------------------------------------
 333   //calculating and storing the final results of flow analysis
 334   if (MCEP) {mcep->Finish(); mcep->WriteHistograms("outputMCEPanalysis.root");}
 335   if(QC) {qc->Finish(); qc->WriteHistograms("outputQCanalysis.root");}
 336   if(GFC) {gfc->Finish(); gfc->WriteHistograms("outputGFCanalysis.root");}
 337   if(FQD) {fqd->Finish(); fqd->WriteHistograms("outputFQDanalysis.root");}
 338   if(LYZ1) {lyz1->Finish(); lyz1->WriteHistograms("outputLYZ1analysis.root");}
 339   if(LYZ2) {lyz2->Finish(); lyz2->WriteHistograms("outputLYZ2analysis.root");}
 340   if(LYZEP) {lyzep->Finish(); lyzep->WriteHistograms("outputLYZEPanalysis.root");}
 341   if(SP) {sp->Finish(); sp->WriteHistograms("outputSPanalysis.root");}
 342   //--------------------------------------------------------------
 343
 344   cout << endl;
 345   cout << " Finished ... " << endl;
 346   cout << endl;
 347
 348   timer.Stop();
 349   cout << endl;
 350   timer.Print();
 351 }