PWG3/muon/AliCFMuonResTask1.C

   1 // DEFINITION OF A FEW CONSTANTS
   2
   3 const Double_t nevtmin= 1;
   4 const Double_t nevtmax = 15000;
   5
   6 // Muons
   7 const Double_t ymin  = -4.0 ;
   8 const Double_t ymax  =  -2.5 ;
   9
  10 const Double_t phimin = -180;
  11 const Double_t phimax = 180;
  12
  13 // Resonance
  14 const Int_t    PDG = 443; //JPsi
  15
  16 const Double_t ptmin =  0.0 ;
  17 const Double_t ptmax =  30 ;
  18 const Double_t pmin =  0.0 ;
  19 const Double_t pmax =  700 ;
  20 const Int_t    charge  = 0 ;
  21 const Double_t mmin =  0.1 ;
  22 const Double_t mmax =  12 ;
  23 const Double_t mymin =  -4 ;
  24 const Double_t mymax =  -2.5 ;
  25 const Double_t costCSmin =  -1.;
  26 const Double_t costCSmax =  1.;
  27 const Double_t costHEmin =  -1.;
  28 const Double_t costHEmax =  1.;
  29 const Double_t phiCSmin =  -TMath::Pi();
  30 const Double_t phiCSmax =  TMath::Pi();
  31 const Double_t phiHEmin =  -TMath::Pi();
  32 const Double_t phiHEmax =  TMath::Pi();
  33
  34 //----------------------------------------------------
  35
  36 Bool_t AliCFMuonResTask1(
  37                             const Bool_t useGrid = 0,
  38                             const Bool_t readAOD = 0,
  39                             const char * kTagXMLFile="wn.xml" // XML file containing tags
  40                             )
  41 {
  42
  43   TBenchmark benchmark;
  44   benchmark.Start("AliMuonResTask1");
  45
  46   AliLog::SetGlobalDebugLevel(0);
  47
  48   Load() ; // load the required libraries
  49
  50   TChain * analysisChain ;
  51
  52 ///// INPUT
  53
  54   if (useGrid) { // data located on AliEn
  55     TGrid::Connect("alien://") ;    //  Create an AliRunTagCuts and an AliEventTagCuts Object
  56                                     //  and impose some selection criteria
  57     AliRunTagCuts      *runCuts   = new AliRunTagCuts();
  58     AliEventTagCuts    *eventCuts = new AliEventTagCuts();
  59     AliLHCTagCuts      *lhcCuts   = new AliLHCTagCuts();
  60     AliDetectorTagCuts *detCuts   = new AliDetectorTagCuts();
  61     eventCuts->SetMultiplicityRange(0,2000);
  62
  63     //  Create an AliTagAnalysis Object and chain the tags
  64     AliTagAnalysis   *tagAna = new AliTagAnalysis();
  65     if (readAOD) tagAna->SetType("AOD");  // for aliroot > v4-05
  66     else         tagAna->SetType("ESD");  // for aliroot > v4-05
  67     TAlienCollection *coll   = TAlienCollection::Open(kTagXMLFile);
  68     TGridResult      *tagResult = coll->GetGridResult("",0,0);
  69     tagResult->Print();
  70     tagAna->ChainGridTags(tagResult);
  71
  72     // Create a new esd chain and assign the chain that is returned by querying the tags
  73     analysisChain = tagAna->QueryTags(runCuts,lhcCuts,detCuts,eventCuts);
  74   }
  75
  76   else {// local data
  77     // here put your input data path
  78     printf("\n\nRunning on local file, please check the path\n\n");
  79
  80     if (readAOD) {
  81       analysisChain = new TChain("aodTree");
  82       analysisChain->Add("AliAOD.root");
  83     }
  84     else {
  85       analysisChain = new TChain("esdTree");
  86       analysisChain->Add("/alidata/alice/arnaldi/CORRFW/polar-1.1000.1000/AliESDs.root");
  87       analysisChain->Add("/alidata/alice/arnaldi/CORRFW/polar-1.1001.1000/AliESDs.root");
  88       analysisChain->Add("/alidata/alice/arnaldi/CORRFW/polar-1.1002.1000/AliESDs.root");
  89    }
  90   }
  91
  92 ///// END INPUT
  93
  94
  95   Info("AliCFMuonResTask1",Form("CHAIN HAS %d ENTRIES",(Int_t)analysisChain->GetEntries()));
  96
  97   // CONTAINER DEFINITION
  98   Info("AliCFMuonResTask1","SETUP CONTAINER");
  99
 100   // The sensitive variables (13 in this example), their indices
 101   UInt_t nevt  = 0;
 102   UInt_t y1  = 1;
 103   UInt_t phi1  = 2;
 104   UInt_t y2  = 3;
 105   UInt_t phi2  = 4;
 106   UInt_t imass  = 5;
 107   UInt_t y  = 6;
 108   UInt_t pt = 7;
 109   UInt_t p = 8;
 110   UInt_t costCS = 9;
 111   UInt_t phiCS = 10;
 112   UInt_t costHE = 11;
 113   UInt_t phiHE = 12;
 114
 115
 116   // Setting up the container grid
 117   UInt_t nstep = 2 ; //number of selection steps : MC and ESD
 118   const Int_t nvar   = 13 ;     //number of variables on the grid
 119   const Int_t nbin1  = nevtmax ;
 120   const Int_t nbin2  = 100 ;
 121   const Int_t nbin3  = 360 ;
 122   const Int_t nbin4  = 100 ;
 123   const Int_t nbin5  = 360 ;
 124   const Int_t nbin6  = 100 ;
 125   const Int_t nbin7  = 50 ;
 126   const Int_t nbin8  = 50 ;
 127   const Int_t nbin9  = 50 ;
 128   const Int_t nbin10  = 20 ;
 129   const Int_t nbin11  = 20 ;
 130   const Int_t nbin12  = 20 ;
 131   const Int_t nbin13  = 20 ;
 132
 133   // arrays for the number of bins in each dimension
 134   Int_t iBin[nvar];
 135   iBin[0]=nbin1;
 136   iBin[1]=nbin2;
 137   iBin[2]=nbin3;
 138   iBin[3]=nbin4;
 139   iBin[4]=nbin5;
 140   iBin[5]=nbin6;
 141   iBin[6]=nbin7;
 142   iBin[7]=nbin8;
 143   iBin[8]=nbin9;
 144   iBin[9]=nbin10;
 145   iBin[10]=nbin11;
 146   iBin[11]=nbin12;
 147   iBin[12]=nbin13;
 148
 149   // arrays for lower bounds :
 150   Double_t *binLim1=new Double_t[nbin1+1];
 151   Double_t *binLim2=new Double_t[nbin2+1];
 152   Double_t *binLim3=new Double_t[nbin3+1];
 153   Double_t *binLim4=new Double_t[nbin4+1];
 154   Double_t *binLim5=new Double_t[nbin5+1];
 155   Double_t *binLim6=new Double_t[nbin6+1];
 156   Double_t *binLim7=new Double_t[nbin7+1];
 157   Double_t *binLim8=new Double_t[nbin8+1];
 158   Double_t *binLim9=new Double_t[nbin9+1];
 159   Double_t *binLim10=new Double_t[nbin10+1];
 160   Double_t *binLim11=new Double_t[nbin11+1];
 161   Double_t *binLim12=new Double_t[nbin12+1];
 162   Double_t *binLim13=new Double_t[nbin13+1];
 163
 164   // values for bin lower bounds
 165   for(Int_t i=0; i<=nbin1; i++) binLim1[i]=(Double_t)nevtmin  + (nevtmax-nevtmin)  /nbin1*(Double_t)i ;
 166   for(Int_t i=0; i<=nbin2; i++) binLim2[i]=(Double_t)ymin  + (ymax-ymin)  /nbin2*(Double_t)i ;
 167   for(Int_t i=0; i<=nbin3; i++) binLim3[i]=(Double_t)phimin  + (phimax-phimin)  /nbin3*(Double_t)i ;
 168   for(Int_t i=0; i<=nbin4; i++) binLim4[i]=(Double_t)ymin  + (ymax-ymin)  /nbin4*(Double_t)i ;
 169   for(Int_t i=0; i<=nbin5; i++) binLim5[i]=(Double_t)phimin  + (phimax-phimin)  /nbin5*(Double_t)i ;
 170   for(Int_t i=0; i<=nbin6; i++) binLim6[i]=(Double_t)mmin  + (mmax-mmin)  /nbin6*(Double_t)i ;
 171   for(Int_t i=0; i<=nbin7; i++) binLim7[i]=(Double_t)mymin  + (mymax-mymin)  /nbin7*(Double_t)i ;
 172   for(Int_t i=0; i<=nbin8; i++) binLim8[i]=(Double_t)ptmin + (ptmax-ptmin)/nbin8*(Double_t)i ;
 173   for(Int_t i=0; i<=nbin9; i++) binLim9[i]=(Double_t)pmin + (pmax-pmin)/nbin9*(Double_t)i ;
 174   for(Int_t i=0; i<=nbin10; i++) binLim10[i]=(Double_t)costCSmin + (costCSmax-costCSmin)/nbin10*(Double_t)i ;
 175   for(Int_t i=0; i<=nbin11; i++) binLim11[i]=(Double_t)phiCSmin + (phiCSmax-phiCSmin)/nbin11*(Double_t)i ;
 176   for(Int_t i=0; i<=nbin12; i++) binLim12[i]=(Double_t)costHEmin + (costHEmax-costHEmin)/nbin12*(Double_t)i ;
 177   for(Int_t i=0; i<=nbin13; i++) binLim13[i]=(Double_t)phiHEmin + (phiHEmax-phiHEmin)/nbin13*(Double_t)i ;
 178
 179   // one container  of 2 steps (MC and ESD) with 12 variables
 180   AliCFContainer* container = new AliCFContainer("container","container for tracks",nstep,nvar,iBin);
 181   // setting the bin limits
 182   container -> SetBinLimits(nevt,binLim1);
 183   container -> SetBinLimits(y1,binLim2);
 184   container -> SetBinLimits(phi1,binLim3);
 185   container -> SetBinLimits(y2,binLim4);
 186   container -> SetBinLimits(phi2,binLim5);
 187   container -> SetBinLimits(imass,binLim6);
 188   container -> SetBinLimits(y,binLim7);
 189   container -> SetBinLimits(pt,binLim8);
 190   container -> SetBinLimits(p,binLim9);
 191   container -> SetBinLimits(costCS,binLim10);
 192   container -> SetBinLimits(phiCS,binLim11);
 193   container -> SetBinLimits(costHE,binLim12);
 194   container -> SetBinLimits(phiHE,binLim13);
 195
 196   // Set list
 197   TList* qaList = new TList();
 198
 199   //CREATE THE CUTS
 200   // Choice of the Resonance
 201   AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts");
 202   mcGenCuts->SetRequirePdgCode(PDG);
 203   mcGenCuts->SetQAOn(qaList);
 204
 205   // Set a pt range of the resonance
 206   AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts");
 207   mcKineCuts->SetChargeMC(charge);
 208   mcKineCuts->SetPtRange(ptmin,ptmax);
 209   mcKineCuts->SetQAOn(qaList);
 210
 211   // Create and fill the list associated
 212   TObjArray* mcList = new TObjArray(0) ;
 213   mcList->AddLast(mcKineCuts);
 214   mcList->AddLast(mcGenCuts);
 215
 216   // kinematic cuts on muons rapidity
 217   AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts","rec-level kine cuts");
 218 //  recKineCuts->SetRapidityRange(ymin,ymax);
 219   recKineCuts->SetRapidityRange(-4,-2.5);
 220   recKineCuts->SetQAOn(qaList);
 221   TObjArray* recList = new TObjArray(0) ;
 222   recList->AddLast(recKineCuts);
 223
 224   // CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
 225
 226   printf("CREATE INTERFACE AND CUTS\n");
 227   AliCFManager* man = new AliCFManager() ;
 228   man->SetParticleContainer     (container);
 229
 230   man->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList);
 231   man->SetParticleCutsList(AliCFManager::kPartAccCuts,recList);
 232
 233   //CREATE THE TASK
 234   printf("CREATE TASK\n");
 235   // create the task
 236   AliCFMuonResTask1 *task = new AliCFMuonResTask1("AliMuonResTask1");
 237   task->SetCFManager(man); //here is set the CF manager
 238   task->SetQAList(qaList);
 239   if (readAOD)       task->SetReadAODData() ;
 240
 241   //SETUP THE ANALYSIS MANAGER TO READ INPUT CHAIN AND WRITE DESIRED OUTPUTS
 242   printf("CREATE ANALYSIS MANAGER\n");
 243   // Make the analysis manager
 244   AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
 245
 246   if (useGrid) mgr->SetAnalysisType(AliAnalysisManager::kGridAnalysis);
 247   else mgr->SetAnalysisType(AliAnalysisManager::kLocalAnalysis);
 248
 249
 250   AliMCEventHandler*  mcHandler = new AliMCEventHandler();
 251   mgr->SetMCtruthEventHandler(mcHandler);
 252
 253   AliInputEventHandler* dataHandler ;
 254
 255   if   (readAOD) dataHandler = new AliAODInputHandler();
 256   else           dataHandler = new AliESDInputHandler();
 257   mgr->SetInputEventHandler(dataHandler);
 258
 259   // Create and connect containers for input/output
 260
 261   // input data
 262   AliAnalysisDataContainer *cinput0  = mgr->CreateContainer("cchain0",TChain::Class(),AliAnalysisManager::kInputContainer);
 263
 264   // output data
 265   Char_t file[256];
 266   sprintf(file,"CFMuonResTask1.root");
 267   printf("Analysis output in %s \n",file);
 268
 269   // output TH1I for event counting
 270   AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("chist0", TH1I::Class(),AliAnalysisManager::kOutputContainer,file);
 271   // output Correction Framework Container (for acceptance & efficiency calculations)
 272   AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("ccontainer0", AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,file);
 273
 274   cinput0->SetData(analysisChain);
 275   mgr->AddTask(task);
 276
 277   mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
 278
 279   mgr->ConnectOutput(task,1,coutput1);
 280   mgr->ConnectOutput(task,2,coutput2);
 281
 282   printf("READY TO RUN\n");
 283   //RUN !!!
 284   if (mgr->InitAnalysis()) {
 285     mgr->PrintStatus();
 286     mgr->StartAnalysis("local",analysisChain);
 287   }
 288
 289   benchmark.Stop("AliMuonResTask1");
 290   benchmark.Show("AliMuonResTask1");
 291
 292   return kTRUE ;
 293 }
 294
 295 void Load() {
 296
 297   //load the required aliroot libraries
 298   gSystem->Load("libANALYSIS") ;
 299   gSystem->Load("libANALYSISalice") ;
 300
 301 //  gSystem->Load("libCORRFW.so") ;
 302   gSystem->Load("$ALICE_ROOT/lib/tgt_linux/libCORRFW.so") ;
 303
 304   //compile online the task class
 305   gSystem->SetIncludePath("-I. -I$ALICE_ROOT/include -I$ROOTSYS/include");
 306   gROOT->LoadMacro("./AliCFMuonResTask1.cxx+");
 307 }