Split: fix refs to AddTaskPhysicsSelection.C
[u/mrichter/AliRoot.git] / PWGUD / multVScentPbPb / MyAnalysisMacroTrackletMulti.C
1 Bool_t needRecPoints = kFALSE;
2
3 void MyAnalysisMacroTrackletMulti
4 (TString dataset="/alice/sim/LHC10f8c_130844",
5  TString outFName = "trbg.root",
6  Int_t   nEvents   = -1,
7  Float_t etaMin     =-0.5,        // min eta range to fill in histos
8  Float_t etaMax     = 0.5,        // max eta range to fill in histos
9  Float_t zMin       = -7,         // process events with Z vertex min
10  Float_t zMax       =  7,         //                     max positions
11  Int_t   useCentVar = 0,          // centrality variable to use: 0=V0, 1=SPD2corr
12  //
13  Float_t cutSigNStd  = 1.5,        // cut on weighed distance used to extract signal
14  Float_t cutSigDPhiS = -1,        // cut on dPhi-phiBent used to extract signal (if negative -> dphi*sqrt(cutSigNStd)
15  Bool_t  useMC  = kTRUE,          // fill MC info (doRec=kTRUE)
16  Float_t scaleMCV0 = 0.7520,      // rescale MC V0 to match the data
17  //
18  Bool_t doRec  = kTRUE,           // fill data histos from new reco
19  Bool_t doInj  = kTRUE,           // create Inj. bg
20  Bool_t doRot  = kFALSE,          // create Rot. bg
21  Bool_t doMix  = kFALSE,//kTRUE,  // create Mix. bg
22  // 
23  // specific parameters for reconstruction
24  float  phiRot      = 3.14159e+00, // angle for bg. generation with rotation
25  float  injScale    = 1.,//0.7,    // inject injScale*Ncl(Lr1/Lr2) hits
26  Bool_t scaleDTheta = kTRUE,       // scale dTheta by 1/sin^2(theta) in trackleting
27  float  nStdDev     = 25.,         // number of st.dev. for tracklet cut to keep
28  float  dphi        = 0.06,        // dphi window (sigma of tracklet cut)
29  float  dtht        = 0.025,       // dtheta .... (if negative, abs will be used with additional cut on |dthetaX|, apart from w.distance
30  float  phishift    = 0.0045,      // bending shift
31  Bool_t remOvl      = kTRUE,       
32  float  ovlPhiCut   = 0.005, 
33  float  ovlZetaCut  = 0.05,
34  Int_t  nEventsSkip = 0,
35  //----------------------- Ntracklets selection parameters important for mixing, to be tuned
36  Float_t ntMin      =   1,         // process events with ESDmult 
37  Float_t ntMax      = 20000,       // within this range
38  Float_t ntMixBinSz = 20000,       // ESDMult bin size for mixing
39  //----------------------- Zv selection parameters important for mixing, to be tuned
40  Float_t zMixBinSz  =  14,       //0.1,  // Zv. bin for mixing
41  //---------------------------------------------------------------------------------
42  //
43  Bool_t checkReconstructables = kFALSE//kTRUE, // fill histos for reconstructable (needs useMC and doRec) 
44  //
45  )
46 {
47   //  
48   if (cutSigDPhiS<0) cutSigDPhiS = TMath::Sqrt(cutSigNStd)*dphi;
49   //
50   needRecPoints = doRec || doInj || doRot || doMix;
51     //
52   printf("Start Analysis for %s, max %d Events skipping %d, Event Cuts: %.1f<eta<%.1f, %.2f<Zv<%.2f\n",
53          dataset.Data(),nEvents,nEventsSkip,etaMin,etaMax,zMin,zMax);
54   printf("Centrality variable: %d\n",useCentVar);
55   printf("Tracklet cuts: dPhi:%.3f dTheta:%.3f phiShift:%.4f | Keep %.1f NstDev\n"
56          "Scale dTheta: %s | Signal Selection: NstDev:%.1f, dPhiS: %.3f\n", 
57          dphi,dtht,phishift,nStdDev,scaleDTheta ? "ON":"OFF",
58          cutSigNStd,cutSigDPhiS);
59   //
60   printf("UseMC: %s. V0 scale: %.4f\n",useMC ? "ON":"OFF",scaleMCV0);
61   printf("Operations:  \n"
62          "Reco:%s (RemOvl:%s phi:%.3f zeta:%.3f)\n"
63          "Inj:%s (scale: %.2f)\n"
64          "Rot:%s (phi: %.4f)\n"
65          "Mix:%s (Nt:%.1f:%.1f:%.1f Zv:%.1f:%.1f:%.2f)\n", 
66          doRec ? "ON":"OFF",remOvl? "ON":"OFF",ovlPhiCut,ovlZetaCut,
67          doInj ? "ON":"OFF",injScale,
68          doRot ? "ON":"OFF",phiRot,
69          doMix ? "ON":"OFF",ntMin,ntMax,ntMixBinSz, zMin,zMax,zMixBinSz);
70   //
71   if (nEvents<0) nEvents = int(1e9);
72   TString format = GetFormatFromDataSet(dataset);
73   //
74   // ALICE stuff
75   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
76   if (!mgr) mgr = new AliAnalysisManager("Test train");
77   //
78   InputHandlerSetup(format,useMC);
79   if (doMix) MixHandlerSetup(ntMin,ntMax,ntMixBinSz, zMin,zMax,zMixBinSz);
80   // compile our task
81   gProof->Load("AliITSMultRecBg.cxx++");
82   gProof->Load("AliTrackletTaskMulti.cxx++");
83   //
84   printf("Loading Centrality task\n");
85   gROOT->LoadMacro("$ALICE_ROOT/OADB/macros/AddTaskCentrality.C");
86   AliCentralitySelectionTask *taskCentrality = AddTaskCentrality();
87   //  taskCentrality->SetDebugLevel(2);
88   if (useMC) taskCentrality->SetMCInput();
89   //  taskCentrality->Dump();
90   //
91   // load and run AddTask macro
92   gROOT->LoadMacro("AddMultTaskTrackletMulti.C");
93   //
94   // create our task
95   AliTrackletTaskMulti *mltTask = AddMultTaskTrackletMulti(outFName.Data());
96   //
97   mltTask->SetUseCentralityVar(useCentVar);
98   mltTask->SetDoNormalReco(doRec);
99   mltTask->SetDoInjection(doInj);
100   mltTask->SetDoRotation(doRot);
101   mltTask->SetDoMixing(doMix);  
102   //
103   mltTask->SetUseMC(useMC);
104   mltTask->SetCheckReconstructables(checkReconstructables);
105   //
106   mltTask->SetEtaMin(etaMin);
107   mltTask->SetEtaMax(etaMax);
108   mltTask->SetZVertexMin(zMin);
109   mltTask->SetZVertexMax(zMax);
110   //
111   mltTask->SetDPhiSCut(cutSigDPhiS);
112   mltTask->SetNStdCut(cutSigNStd);
113   mltTask->SetScaleMCV0(scaleMCV0);
114   //
115   mltTask->SetScaleDThetaBySin2T(scaleDTheta);
116   mltTask->SetNStdDev(nStdDev);
117   mltTask->SetPhiWindow(dphi);
118   mltTask->SetThetaWindow(dtht);
119   mltTask->SetPhiShift(phishift);
120   mltTask->SetPhiOverlapCut(ovlPhiCut);
121   mltTask->SetZetaOverlapCut(ovlZetaCut);
122   mltTask->SetPhiRot(phiRot);
123   mltTask->SetInjScale(injScale);
124   mltTask->SetRemoveOverlaps(remOvl);
125   //
126   printf("new Task: %p\n",mltTask);
127   //
128   printf("Requesting physics selection in %s mode\n",useMC ? "MC":"Data");
129   gROOT->ProcessLine(".L $ALICE_ROOT/OADB/macros/AddTaskPhysicsSelection.C");
130   AliPhysicsSelectionTask* physicsSelectionTask = AddTaskPhysicsSelection(useMC,0);
131   mltTask->SelectCollisionCandidates();//AliVEvent::kMB);
132   //
133   // Run analysis
134   mgr->InitAnalysis();
135   // process dataset  
136   mgr->StartAnalysis("proof", dataset.Data(), nEvents, nEventsSkip); 
137   //
138   TString evstCmd = "if [ -e event_stat.root ]; then \nmv event_stat.root evstat_"; 
139   evstCmd += outFName;  evstCmd += " \nfi";
140   gSystem->Exec( evstCmd.Data() );
141   
142 }
143
144
145 TString GetFormatFromDataSet(TString dataset) {
146   
147 //   Info("runAAF.C","Detecting format from dataset (may take while, depends on network connection)...");
148   TString dsTreeName;
149   if (dataset.Contains("#")) {
150     Info("runAAF.C",Form("Detecting format from dataset name '%s' ...",dataset.Data()));
151     dsTreeName=dataset(dataset.Last('#'),dataset.Length());
152   } else {
153     Info("runAAF.C",Form("Detecting format from dataset '%s' (may take while, depends on network connection) ...",dataset.Data()));
154     TFileCollection *ds = gProof->GetDataSet(dataset.Data());
155     if (!ds) {
156       Error(Form("Dataset %s doesn't exist on proof cluster!!!!",dataset.Data()));
157       return "";
158     }
159     dsTreeName = ds->GetDefaultTreeName();
160   }
161
162   if (dsTreeName.Contains("esdTree")) {
163     Info("runAAF.C","ESD input format detected ...");
164     return "ESD";
165   } else if (dsTreeName.Contains("aodTree"))  {
166     Info("runAAF.C","AOD input format detected ...");
167     return "AOD";
168   } else {
169     Error("runAAF.C",Form("Tree %s is not supported !!!",dsTreeName.Data()));
170     Error("runAAF.C",Form("Maybe set your DS to %s#esdTree or %s#aodTree",dataset.Data(),dataset.Data()));
171   }
172   
173   return "";
174 }
175
176 Bool_t InputHandlerSetup(TString format = "esd", Bool_t useKine = kTRUE)
177 {
178   format.ToLower();
179
180   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
181
182   AliAnalysisDataContainer *cin = mgr->GetCommonInputContainer();
183
184   if (cin) return;
185
186   if (!format.CompareTo("esd"))
187   {
188     AliESDInputHandler *esdInputHandler = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
189
190     if (!esdInputHandler)
191     {
192       Info("CustomAnalysisTaskInputSetup", "Creating esdInputHandler ...");
193       if (needRecPoints)
194         esdInputHandler = new AliESDInputHandlerRP();
195       else 
196         esdInputHandler = new AliESDInputHandler();
197       //
198       mgr->SetInputEventHandler(esdInputHandler);
199     }
200     //
201     if (useKine)
202     {
203       AliMCEventHandler* mcInputHandler = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
204
205       if (!mcInputHandler)
206       {
207         Info("CustomAnalysisTaskInputSetup", "Creating mcInputHandler ...");
208         AliMCEventHandler* mcInputHandler = new AliMCEventHandler();
209         mgr->SetMCtruthEventHandler(mcInputHandler);
210       }
211     }
212
213   }
214   else if (!format.CompareTo("aod"))
215   {
216     AliAODInputHandler *aodInputHandler = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
217
218     if (!aodInputHandler)
219     {
220       Info("CustomAnalysisTaskInputSetup", "Creating aodInputHandler ...");
221       aodInputHandler = new AliAODInputHandler();
222       mgr->SetInputEventHandler(aodInputHandler);
223     }
224   }
225   else
226   {
227     Info("Wrong input format!!! Only ESD and AOD are supported. Skipping Task ...");
228     return kFALSE;
229   }
230
231   return kTRUE;
232 }
233
234 void MixHandlerSetup(float ntMin,float ntMax,float ntMixBinSz,
235                      float zMin, float zMax, float zMixBinSz)
236 {
237   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
238   if (!mgr) return;
239   int bufferSize = 1;
240   AliESDInputHandlerRP *esdH = dynamic_cast<AliESDInputHandlerRP*>(mgr->GetInputEventHandler());
241   if (!esdH) return;
242   //
243   AliMixEventInputHandler *esdMixH = new AliMixEventInputHandler(bufferSize);
244   esdMixH->SetInputHandlerForMixing(esdH);
245   AliMixEventPool *evPool = new AliMixEventPool("MyPool");
246   AliMixEventCutObj *tracklets = new AliMixEventCutObj(AliMixEventCutObj::kNumberTracklets, ntMin,ntMax,ntMixBinSz);
247   AliMixEventCutObj *zvertex = new AliMixEventCutObj(AliMixEventCutObj::kZVertex, zMin,zMax, zMixBinSz);
248   //  evPool->AddCut(tracklets);
249   evPool->AddCut(zvertex);
250   //evPool->Init();
251   evPool->Print();
252   esdMixH->SetEventPool(evPool);
253   esdH->SetMixingHandler(esdMixH);
254 }
255
256 void AddPhysicsSelection(Bool_t isMC)
257 {
258   // physics selection a la Michele
259   if(!isMC) {
260     //AliPhysicsSelection * physSel = physicsSelectionTask->GetPhysicsSelection();
261     //    physSel->AddCollisionTriggerClass("+CMBAC-B-NOPF-ALL");
262     /*
263     physSel->AddCollisionTriggerClass("+CMBS1C-B-NOPF-ALL");
264     physSel->AddCollisionTriggerClass("+CMBS1A-B-NOPF-ALL");
265     */
266     //
267     //    physSel->AddCollisionTriggerClass("+CMBS2C-B-NOPF-ALL");
268     //    physSel->AddCollisionTriggerClass("+CMBS2A-B-NOPF-ALL");
269     //
270     // This are needed only to fill the statistics tables
271     //    physSel->AddBGTriggerClass("+CMBAC-C-NOPF-ALL");
272     //    physSel->AddBGTriggerClass("+CMBAC-A-NOPF-ALL");
273     //    physSel->AddBGTriggerClass("+CMBAC-E-NOPF-ALL");
274     //
275     /*
276     physSel->AddBGTriggerClass("+CMBS1C-C-NOPF-ALL");
277     physSel->AddBGTriggerClass("+CMBS1C-A-NOPF-ALL");
278     physSel->AddBGTriggerClass("+CMBS1C-E-NOPF-ALL");
279     //
280     physSel->AddBGTriggerClass("+CMBS1A-C-NOPF-ALL");
281     physSel->AddBGTriggerClass("+CMBS1A-A-NOPF-ALL");
282     physSel->AddBGTriggerClass("+CMBS1A-E-NOPF-ALL");
283     //
284     */
285     /*
286     //
287     physSel->AddBGTriggerClass("+CMBS2C-C-NOPF-ALL");
288     physSel->AddBGTriggerClass("+CMBS2C-A-NOPF-ALL");
289     physSel->AddBGTriggerClass("+CMBS2C-E-NOPF-ALL");
290     //
291     physSel->AddBGTriggerClass("+CMBS2A-C-NOPF-ALL");
292     physSel->AddBGTriggerClass("+CMBS2A-A-NOPF-ALL");
293     physSel->AddBGTriggerClass("+CMBS2A-E-NOPF-ALL");
294     */
295   } 
296   // if you use the following line, your task only gets the selected events
297   //  task->SelectCollisionCandidates(AliVEvent::kUserDefined);
298   //  task->SelectCollisionCandidates();
299   //
300   //Alternatively, in the UserExec of your task:
301   //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kUserDefined);
302   //
303 }