87e7c1ffedfb6bf54695a6e864b05fe0303bae87
[u/mrichter/AliRoot.git] / ITS / macrosSDD / ITSQArecoparam.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TFile.h>
3 #include <TAlienFile.h>
4 #include <TSystem.h>
5 #include <TGrid.h>
6 #include <TGridResult.h>
7 #include <Riostream.h>
8 #include <TObjArray.h>
9 #include <TKey.h>
10 #include <TStyle.h>
11 #include <TTree.h>
12 #include <TROOT.h>
13 #include <TPluginManager.h>
14 #include <TClass.h>
15 #include <iostream>
16 #include <stdio.h>
17 #include <TGeoGlobalMagField.h>
18 #include <TMap.h>
19 #include "AliRawReader.h"
20 #include "AliRawReaderRoot.h"
21 #include "AliRawReaderDate.h"
22 #include "AliITSQADataMakerRec.h"
23 #include "AliITSQASDDDataMakerRec.h"
24 #include "AliITSQAChecker.h"
25 #include "AliQAChecker.h"
26 #include "AliITSQASDDChecker.h"
27 #include "AliReconstructor.h"
28 #include "AliCDBManager.h"
29 #include "AliQAv1.h"
30 #include "AliMagF.h"
31 #include "AliGeomManager.h"
32 #include "AliITSInitGeometry.h"
33 #include "AliITSgeom.h"
34 #include "AliRecoParam.h"
35 #include "AliCDBPath.h"
36 #include "AliCDBEntry.h"
37 #include "AliRecoParam.h"
38 #include "AliDetectorRecoParam.h"
39 #include "AliITSReconstructor.h"
40 #include "AliITSRecPointContainer.h"
41 #include "AliLog.h"
42 #include "AliGRPObject.h"
43 #include "AliRunInfo.h"
44 #endif
45 void ITSQArecoparam(char *iFile, Int_t runNb=150000, Int_t idet=2, Int_t FirstEvt=0, Int_t MaxEvts=1000000)
46 {
47   TString namefile(iFile);
48   if(namefile.Contains("alien"))
49     {
50       TGrid::Connect("alien://"); 
51       if(!gGrid) {
52         printf("gGrid not found! exit macro\n");
53         return;
54       }
55     }
56   gStyle->SetPalette(1);
57   Int_t ic=0;
58   AliRawReader *rd=NULL; 
59   if(strstr(iFile,".root")!=0){rd = new AliRawReaderRoot(iFile,FirstEvt);}
60   else{rd=new AliRawReaderDate(iFile,FirstEvt);}
61   //Int_t runNumber = rd->GetRunNumber();
62   cout << "ITS Quality Assurance Prototype for run "<< runNb << endl; 
63   //TStopwatch mytimer;
64   //TString namefile(iFile);
65   // Set OCDB if needed
66   AliCDBManager* man = AliCDBManager::Instance();
67   if (!man->IsDefaultStorageSet()) {
68     printf("Setting a local default storage and run number \n");
69     //  if(namefile.Contains("alien")){
70       man->SetDefaultStorage("alien://folder=/alice/data/2011/OCDB");
71       //}
72       //else{man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");}
73     man->SetRun(runNb);
74     AliQAv1::SetQARefStorage("local://$ALICE_ROOT/QARef") ;
75   }
76   AliCDBManager::Instance()->GetAll(Form("ITS/Calib/*"));
77   
78
79   AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
80
81   AliGRPObject *fGRPData=NULL;
82   
83   TMap* m = dynamic_cast<TMap*>(entry->GetObject());  // old GRP entry
84
85     if (m) {
86        printf("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject");
87        m->Print();
88        fGRPData = new AliGRPObject();
89        fGRPData->ReadValuesFromMap(m);
90     }
91
92     else {
93        printf("Found an AliGRPObject in GRP/GRP/Data, reading it");
94        fGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject());  // new GRP entry
95        entry->SetOwner(0);
96     }
97
98
99  if (!fGRPData) {
100      printf("Error   No GRP entry found in OCDB!");
101      return;
102   }
103
104
105   TString lhcState = fGRPData->GetLHCState();
106   if (lhcState==AliGRPObject::GetInvalidString()) {
107     printf("Error  GRP/GRP/Data entry:  missing value for the LHC state ! Using UNKNOWN");
108     lhcState = "UNKNOWN";
109   }
110
111   TString beamType = fGRPData->GetBeamType();
112   if (beamType==AliGRPObject::GetInvalidString()) {
113     printf("Error   GRP/GRP/Data entry:  missing value for the beam type ! Using UNKNOWN");
114     beamType = "UNKNOWN";
115   }
116
117   Float_t beamEnergy = fGRPData->GetBeamEnergy();
118   if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
119     printf("Error   GRP/GRP/Data entry:  missing value for the beam energy ! Using 0");
120     beamEnergy = 0;
121   }
122
123   TString runType = fGRPData->GetRunType();
124   if (runType==AliGRPObject::GetInvalidString()) {
125     printf("Error  GRP/GRP/Data entry:  missing value for the run type ! Using UNKNOWN");
126     runType = "UNKNOWN";
127   }
128
129   Int_t activeDetectors = fGRPData->GetDetectorMask();
130   if (activeDetectors==AliGRPObject::GetInvalidUInt()) {
131     printf("Error  GRP/GRP/Data entry:  missing value for the detector mask ! Using 1074790399");
132     activeDetectors = 1074790399;
133   }
134
135   AliRunInfo *fRunInfo = new AliRunInfo(lhcState, beamType, beamEnergy, runType, activeDetectors);
136   fRunInfo->Dump();
137
138  if ( TGeoGlobalMagField::Instance()->IsLocked() ) {
139     if (TGeoGlobalMagField::Instance()->GetField()->TestBit(AliMagF::kOverrideGRP)) {
140       printf("ExpertMode!!! GRP information will be ignored !");
141       printf("ExpertMode!!! Running with the externally locked B field !");
142     }
143     else {
144       printf("Destroying existing B field instance!");
145       delete TGeoGlobalMagField::Instance();
146     }    
147   }
148   if ( !TGeoGlobalMagField::Instance()->IsLocked() ) {
149     // Construct the field map out of the information retrieved from GRP.
150     Bool_t ok = kTRUE;
151     // L3
152     Float_t l3Current = fGRPData->GetL3Current((AliGRPObject::Stats)0);
153     if (l3Current == AliGRPObject::GetInvalidFloat()) {
154       printf("Error :  GRP/GRP/Data entry:  missing value for the L3 current !");
155       ok = kFALSE;
156     }
157     
158     Char_t l3Polarity = fGRPData->GetL3Polarity();
159     if (l3Polarity == AliGRPObject::GetInvalidChar()) {
160       printf("Error:   GRP/GRP/Data entry:  missing value for the L3 polarity !");
161       ok = kFALSE;
162     }
163
164    // Dipole
165     Float_t diCurrent = fGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
166     if (diCurrent == AliGRPObject::GetInvalidFloat()) {
167       printf("Error  GRP/GRP/Data entry:  missing value for the dipole current !");
168       ok = kFALSE;
169     }
170
171     Char_t diPolarity = fGRPData->GetDipolePolarity();
172     if (diPolarity == AliGRPObject::GetInvalidChar()) {
173       printf("Error GRP/GRP/Data entry:  missing value for the dipole polarity !");
174       ok = kFALSE;
175     }
176
177
178    Int_t  polConvention = fGRPData->IsPolarityConventionLHC() ? AliMagF::kConvLHC : AliMagF::kConvDCS2008;
179    Bool_t uniformB = fGRPData->IsUniformBMap();
180    
181    if (ok) { 
182      AliMagF* fld = AliMagF::CreateFieldMap(TMath::Abs(l3Current) * (l3Polarity ? -1:1), 
183                                             TMath::Abs(diCurrent) * (diPolarity ? -1:1), 
184                                             polConvention,uniformB,beamEnergy, beamType.Data());
185      if (fld) {
186        TGeoGlobalMagField::Instance()->SetField( fld );
187        TGeoGlobalMagField::Instance()->Lock();
188        printf("Running with the B field constructed out of GRP !");
189      }
190      else printf("Fatal Failed to create a B field map !");
191    }
192    else printf("Fatal  B field is neither set nor constructed from GRP ! Exitig...");
193   }
194   
195
196
197
198   AliITSQADataMakerRec *itsQAdm = new AliITSQADataMakerRec(kTRUE,idet,0); //online kTRUE
199   itsQAdm->SetWriteExpert() ;
200   itsQAdm->SetRunNumber(runNb);  
201   //________________________For the RecPoints____________________________________
202   /************************************************/
203   TPluginManager* pluginManager=NULL;
204   TPluginHandler* pluginHandler=NULL; 
205   AliReconstructor* reconstructor = NULL;
206   AliITSRecPointContainer* rpcont=NULL;
207   AliGeomManager::LoadGeometry("geometry.root");
208   AliGeomManager::ApplyAlignObjsFromCDB("ITS");
209   //   ITS initializations
210   
211   AliITSInitGeometry initgeom;
212   AliITSgeom *geom = initgeom.CreateAliITSgeom();
213   printf("Geometry name: %s\n",(initgeom.GetGeometryName()).Data());
214   
215   printf("Loading reconstruction parameter objects for detector ITS\n");
216   AliRecoParam fRecoParam; 
217   AliCDBPath path("ITS","Calib","RecoParam");
218   AliCDBEntry *entry2=AliCDBManager::Instance()->Get(path.GetPath());
219   Bool_t cacheStatus = AliCDBManager::Instance()->GetCacheFlag();
220   if(!entry2){printf("Couldn't find RecoParam entry in OCDB for detector ITS");entry2=NULL;}
221   else {
222     TObject *recoParamObj = entry2->GetObject();
223     if (dynamic_cast<TObjArray*>(recoParamObj)) {
224       printf("RecoParam TObjArray\n");
225       fRecoParam.AddDetRecoParamArray(0,dynamic_cast<TObjArray*>(recoParamObj));
226     }
227     else if (dynamic_cast<AliDetectorRecoParam*>(recoParamObj)) {
228       printf("RecoParam AliDetectorRecoParam\n");
229       printf("Single set of reconstruction parameters found for detector ITS");
230       (dynamic_cast<AliDetectorRecoParam*>(recoParamObj))->SetAsDefault();
231       fRecoParam.AddDetRecoParam(0,(dynamic_cast<AliDetectorRecoParam*>(recoParamObj)));
232     }
233     else {printf("Error: No valid RecoParam object found in the OCDB for detector ITS");}
234     entry2->SetOwner(0);
235   }
236   if(!cacheStatus)entry2->SetObject(NULL);
237   if(!cacheStatus){ delete entry2;}
238   
239   // load the reconstructor object
240  pluginManager = gROOT->GetPluginManager();
241   TString detName = "ITS";
242   TString recName = "Ali" + detName + "Reconstructor";  
243   
244 pluginHandler = pluginManager->FindHandler("AliReconstructor", "ITS");
245   // if not, add a plugin for it
246   if (!pluginHandler) {
247     printf("defining plugin for ITS\n");
248     TString libs = gSystem->GetLibraries();
249     if (libs.Contains("lib" + detName + "base.so") ||
250         (gSystem->Load("lib" + detName + "base") >= 0)) {pluginManager->AddHandler("AliReconstructor", detName,recName, detName + "rec", recName + "()");}
251     else {pluginManager->AddHandler("AliReconstructor", detName,recName, detName, recName + "()");}
252     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
253   }
254   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);}
255   if (fRecoParam.GetDetRecoParamArray(0) && !AliReconstructor::GetRecoParam(0)) {
256     const AliDetectorRecoParam *par = fRecoParam.GetDetRecoParam(0);
257     reconstructor->Init();
258     reconstructor->SetRecoParam(par);
259   }
260   
261   /*AliITSRecPointContainer**/ rpcont=AliITSRecPointContainer::Instance();
262   rpcont->PrepareToRead();
263
264   Int_t cycleLength = 5;  
265   //cout << "Processing Run " << runNumber << endl;
266   cout << "Init: " << AliQAv1::kRAWS << endl;
267   
268   TObjArray **objArray= itsQAdm->Init(AliQAv1::kRAWS, cycleLength);
269   cout<<"raw tobjarray :"<<objArray<<"\n"<<endl;
270   for(Int_t spec = 0 ; spec < 5 ; spec++){
271     if(spec==1){
272       itsQAdm->SetEventSpecie(AliRecoParam::kCosmic);
273       AliQAv1::Instance()->SetEventSpecie(AliRecoParam::kCosmic);
274       itsQAdm->InitRaws();
275     }
276     else{continue;}
277   }
278   itsQAdm->StartOfCycle(AliQAv1::kRAWS,runNb,kFALSE);
279   /*********************************************************************/
280
281     cout << "Init: " << AliQAv1::kRECPOINTS << endl;
282     TObjArray **objArray1=  itsQAdm->Init(AliQAv1::kRECPOINTS, cycleLength);
283     cout<<"recpoint tobjarray :"<<objArray1<<"\n"<<endl;
284     for(Int_t spec = 0 ; spec < 5 ; spec++){
285       if(spec==1){
286         itsQAdm->SetEventSpecie(AliRecoParam::kCosmic);
287         AliQAv1::Instance()->SetEventSpecie(AliRecoParam::kCosmic);
288         itsQAdm->InitRecPoints();
289       }
290       else{continue;}
291     }
292     itsQAdm->StartOfCycle(AliQAv1::kRECPOINTS,runNb,kTRUE);
293
294   /*********************************************************************/
295   Int_t iev = 0;
296   while(rd->NextEvent() && iev < MaxEvts ) {
297     cout<<">>>>>>>   Processing event number: "<<++iev<<endl;
298     /*******************************************************/
299     if(itsQAdm->IsCycleDone()) {
300       cout << "end of cycle" << endl;
301       AliQAChecker::Instance()->SetRunNumber(AliCDBManager::Instance()->GetRun());      
302       itsQAdm->EndOfCycle(AliQAv1::kRAWS);
303       itsQAdm->StartOfCycle(AliQAv1::kRAWS,ic++,kFALSE);
304     }
305     /******************************************************/
306     /*************************************************/
307
308       if(itsQAdm->IsCycleDone()) {
309         cout << "end of cycle" << endl;
310         AliQAChecker::Instance()->SetRunNumber(AliCDBManager::Instance()->GetRun());
311         itsQAdm->EndOfCycle(AliQAv1::kRECPOINTS);
312         itsQAdm->StartOfCycle(AliQAv1::kRECPOINTS,ic,kTRUE);
313       } 
314  
315     /*************************************************/
316     cout<<"Beginning Exec"<<endl;
317     cout<<"AliQAv1::kRAWS   "<<AliQAv1::kRAWS<<endl;
318     itsQAdm->SetEventSpecie(AliRecoParam::kCosmic);
319     AliQAv1::Instance()->SetEventSpecie(AliRecoParam::kCosmic);
320     itsQAdm->Exec(AliQAv1::kRAWS,rd);
321     /***************************************************/
322     //  if(kRecPoints){
323       cout<<"AliQAv1::kRECPOINTS   "<<AliQAv1::kRECPOINTS<<endl;
324       cout << "DigitsToRecPoints" << endl;
325       TTree* fTreeR = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
326       Char_t option[5];
327
328       if(idet==0)sprintf(option,"ALL");
329       else if(idet==1)sprintf(option,"SPD");
330       else if(idet==2)sprintf(option,"SDD");
331       else if(idet==3)sprintf(option,"SSD");
332       printf("\t\t===========>option is %s\n",option);
333
334       rpcont->PrepareToRead();
335       reconstructor->Reconstruct(rd,fTreeR);
336
337       itsQAdm->SetEventSpecie(AliRecoParam::kCosmic);
338       AliQAv1::Instance()->SetEventSpecie(AliRecoParam::kCosmic)  ; 
339       itsQAdm->Exec(AliQAv1::kRECPOINTS,fTreeR);
340       cout<<"Finishing Exec"<<endl;
341     
342       ((AliITSReconstructor*)reconstructor)->ResetRecPoints();
343       delete fTreeR; 
344       fTreeR=NULL; 
345       /*****************************************************/
346     }
347
348   cout << "end RAWS cycle: " << AliQAv1::kRAWS << endl;
349   cout << "refStorage: " << AliQAv1::GetQARefStorage() << endl;
350   cout << "end of cycle 2" << endl;
351   AliQAChecker::Instance()->SetRunNumber(AliCDBManager::Instance()->GetRun());
352   itsQAdm->EndOfCycle(AliQAv1::kRAWS);   
353   cout << "Raws QA completed for " << iev << " events" << endl;
354   /*******************************************************************/
355
356     AliQAChecker::Instance()->SetRunNumber(AliCDBManager::Instance()->GetRun());
357     itsQAdm->EndOfCycle(AliQAv1::kRECPOINTS);
358     cout << "RecPoints QA completed for " << iev << " events" << endl;
359
360   /*******************************************************************/
361   itsQAdm->Finish(); // write to the output File
362
363   cout << "Call AliITSQASDDDataMakerRec destructor" << endl;
364   delete itsQAdm;
365   itsQAdm=NULL;
366
367 }
368
369 void ITSQArecoparam(Int_t runNb=150000,Int_t year=2011,Char_t period[10]="LHC11a",Int_t idet=2,Int_t FirstEvt=0, Int_t MaxEvts=1000000)
370 {
371   //TString namefile(iFile);
372   //if(namefile.Contains("alien"))
373   //{
374       TGrid::Connect("alien://");
375       if(!gGrid) {
376         printf("gGrid not found! exit macro\n");
377         return;
378       }
379       //}
380
381   gSystem->AddIncludePath("-I. -I$ALICE_ROOT/include");
382
383   char path2[200];
384
385   sprintf(path2,"/alice/data/");
386   printf("path %s\n",path2);
387   char rawfilename[200];
388   sprintf(rawfilename,"%04i/%s/%09i/raw/%02i%09i*.*.root",year,period,runNb,year-2000,runNb);
389
390
391   Int_t number=0;
392
393   TGridResult *gr = gGrid->Query(path2,rawfilename);
394   if (gr->GetEntries() < 1) {
395     printf("In this run there are not raws files: Exit macro\n");
396     return;
397   }
398   //if(gr->GetEntries() > 10) number=3;
399   //else number=3;
400   Bool_t chunkok=kFALSE;
401   TString filenamedef; //= gr->GetKey(number,"turl");
402   //printf("-> FILE %s \n",filename);
403   while(chunkok==kFALSE || number<gr->GetEntries())
404     {
405       const char* filename = gr->GetKey(number,"turl");
406       printf("-> FILE %s \n",filename);
407       TString namefile(filename);
408       if(namefile.Contains(".10.root")==kFALSE){chunkok=kTRUE; filenamedef.Form("%s",filename); break;} 
409       else{
410         chunkok=kFALSE;
411         number++;
412       }
413       if(number==gr->GetEntries()-1) {chunkok=kTRUE; filenamedef.Form("%s",filename); break;}       
414     }
415
416   char *filetouse=(char*)filenamedef.Data();//=filenamedef.Data(); //=NULL;
417   //sprintf(filetouse,"%s",(char*)filenamedef.Data());
418
419   printf("File to use = %s\n\n",filetouse);
420
421   Int_t runn=runNb;
422
423   printf("Run: %d \n",runn);
424
425   ITSQArecoparam(filetouse,runn,idet,FirstEvt,MaxEvts);
426
427 }