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") ||
 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 }